#include "real.h"
#include <stdio.h>
#include <math.h>
#include "qr.h"
#include "mat.h"
//#include <gsl/gsl_math.h>
//#include <gsl/gsl_eigen.h>

void qrdiag(real**,int);

int main(void){
int i,n=4,j,k,nn;
real **r,**a,**q,**b,**c,*x,*y,s,eps=0.00001;

a = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) a[i] = (real *)malloc(n*sizeof(real)); 

r = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) r[i] = (real *)malloc(n*sizeof(real)); 

q = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) q[i] = (real *)malloc(n*sizeof(real)); 

b = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) b[i] = (real *)malloc(n*sizeof(real)); 

c = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) c[i] = (real *)malloc(n*sizeof(real)); 

x = (real *)malloc(n*sizeof(real)); 
y = (real *)malloc(n*sizeof(real)); 

//for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j]=1/(real)(i+j+1);
srand(2);
for(i=0;i<n;i++) for(j=i;j<n;j++)
	//{a[i][j]=(real)(i+1)*(j+1)*(j+1)+1;a[j][i]=a[i][j];}
	{a[i][j]=((real)rand())/RAND_MAX;a[j][i]=a[i][j];}

//	a[0][0]=1   ; a[1][0]=1   ; a[2][0]=0.3 ; // a[3][0]=4 ;
//	a[0][1]=1   ; a[1][1]=2   ; a[2][1]=1.2 ; // a[3][1]=3 ;
//	a[0][2]=0.3 ; a[1][2]=1.2 ; a[2][2]=2   ; // a[3][2]=-3 ;
//	a[0][3]=4   ; a[1][3]=3   ; a[2][3]=-3  ; // a[3][3]=3;

printf("matrix A:\n"); print_mat(a,n,n); lanczos(a,n,x,y,q);
printf("matrix Q:\n"); print_mat(q,n,n);

for(i=0;i<n;i++) for(j=0;j<n;j++) b[i][j]=vec_vec(q[i],q[j],n);
printf("matrix Q'.Q:\n"); print_mat(b,n,n);

for(i=0;i<n;i++) for(j=0;j<n;j++) b[i][j]=0;
for(i=0;i<n;i++) b[i][i]=x[i];
for(i=1;i<n;i++) {b[i][i-1]=y[i]; b[i-1][i]=y[i];};
mat_mat(a,n,n,q,n,r); mat_mat(q,n,n,b,n,c);

printf("matrix AQ:\n"); print_mat(r,n,n);
printf("matrix QT:\n"); print_mat(c,n,n);
printf("matrix T :\n"); print_mat(b,n,n);
qrdiag(b,n);
printf("matrix D :\n"); print_mat(b,n,n);

printf("eigenvalues: "); for(i=0;i<n;i++) printf("%g ",b[i][i]);
printf("\n");
}