#include "real.h"
#include "mat.h"
#include <stdlib.h>
#include <math.h>

void lanczos(real **A,int n,real *a,real *b,real **q){
int i,j; real *w, norm; w=(real*)malloc(n*sizeof(real));
for(i=0;i<n;i++)for(j=0;j<n;j++)
	if(A[i][j]!=A[j][i]) printf("lanczos: nonsymmetric matrix");

for(i=0;i<n;i++) q[0][i]=((real)rand())/RAND_MAX;

norm=sqrt(vec_vec(q[0],q[0],n)); for(i=0;i<n;i++) q[0][i]/=norm;
mat_vec(A,n,n,q[0],w); a[0]=vec_vec(q[0],w,n);

for(i=0;i<n-1;i++){

if(i==0) for(j=0;j<n;j++) q[i+1][j]=w[j]-a[i]*q[i][j];
    else for(j=0;j<n;j++) q[i+1][j]=w[j]-a[i]*q[i][j]-b[i]*q[i-1][j];

b[i+1]=sqrt(vec_vec(q[i+1],q[i+1],n)); for(j=0;j<n;j++) q[i+1][j]/=b[i+1];

mat_vec(A,n,n,q[i+1],w); a[i+1]=vec_vec(q[i+1],w,n);
		}
free(w);
}