#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define TIME ((double)clock())/CLOCKS_PER_SEC
#define RND ((double)rand()/RAND_MAX)
#include "real.h"
#define max_print 20

void a_times_b(real** a,int n,int m,real** b,int l,real** c){
  for(int i=0;i<n;i++){
		for(int j=0;j<l;j++){
			c[j][i]=0;
			for(int k=0;k<m;k++) c[j][i]+=a[k][i]*b[j][k]; } } }

void aT_times_b(real** a,int n,int m,real** b,int l,real** c){
  for(int i=0;i<m;i++){
		for(int j=0;j<l;j++){
			c[j][i]=0;
			for(int k=0;k<n;k++) c[j][i]+=a[i][k]*b[j][k]; } } }

void print_matrix(real** a,int n,int m){
	for(int i=0;i<n;i++){
		for(int j=0;j<m;j++) printf("%.4f\t",a[j][i]); printf("\n"); } }

int qrdec(double** a,double** r,int n,int m);
int jacobi(real **A, int n, real** V);

int main(int argc, char** argv){
	int n=(argc>1? atoi(argv[1]):5);

	real** A=(real**)calloc(n,sizeof(real*));
	for(int i=0;i<n;i++) A[i]=(real*)calloc(n,sizeof(real));
	for(int i=0;i<n;i++) A[i][i]=(real)(i+1);

	real** V=(real**)calloc(n,sizeof(real*));
	for(int i=0;i<n;i++) V[i]=(real*)calloc(n,sizeof(real));

	real** R=(real**)calloc(n,sizeof(real*));
	for(int i=0;i<n;i++) R[i]=(real*)calloc(n,sizeof(real));

	for(int i=0;i<n;i++)for(int j=0;j<n;j++)V[j][i]=RND;
	qrdec(V,R,n,n);

	printf("\nmatrix A: \n");
	if(n<max_print) print_matrix(A,n,n);

	aT_times_b(V,n,n,A,n,R);
	a_times_b(R,n,n,V,n,A);

	printf("\nmatrix A in a random representation: \n");
	if(n<max_print) print_matrix(A,n,n);

	real time;
	time=TIME;
	int rotations=jacobi(A,n,V);
	time-=TIME;
	printf("\nthe result of Jacobi diagonalization: \n");
	printf("time\t = %g\n",time);
	printf("rotations\t = %d\n",rotations);
	if(n<max_print) print_matrix(A,n,n);
	else printf("Ann=%g\n",A[n-1][n-1]);
	return 0;
}