#include "matutil.c"
#include "house.c"
void tridiagonalize(int n,matrix a,vector d,vector e)
{
	int i,j,k;
	double s,t,p,q;
	vector v,w;
	
	for(k=0;k<n-2;k++){
		v=a[k];d[k]=v[k];
		e[k]=house(n-k-1,&v[k+1]);
		if(e[k]==0) continue;
		for(i=k+1;i<n;i++){
			s=0;
			for(j=k+1;j<i;j++) s+=a[j][i]*v[j];
			for(j=i  ;j<n;j++) s+=a[i][j]*v[j];
			d[i]=s;
		}
		t=innerproduct(n-k-1,&v[k+1],&d[k+1])/2;
		for(i=n-1;i>k;i--){
			p=v[i];q=d[i]-t*p;d[i]=q;
			for(j=i;j<n;j++)
					a[i][j]-=p*d[j]+q*v[j];
		}
	}
	if(n>=2){d[n-2]=a[n-2][n-2];
			 e[n-2]=a[n-2][n-1];}
	if(n>=1) d[n-1]=a[n-1][n-1];
	for(k=n-1;k>=0;k--){
		v=a[k];
		if(k<n-2){
			for(i=k+1;i<n;i++){
				v=a[i];
				t=innerproduct(n-k-1,&v[k+1],&w[k+1]);
				for(j=k+1;j<n;j++)
					w[j]-=t*v[j];
			}
		}
		for(i=0;i<n;i++) v[i]=0;
		v[k]=1;
	}
}
void main(void){
	matrix mx=new_matrix(4,4);
	vector v1=new_vector(4);
	vector v2=new_vector(4);
	int i,j;
	mx[0][0]=1;
	mx[0][1]=1;
	mx[0][2]=1;
	mx[0][3]=1;
	
	mx[1][0]=1;
	mx[1][1]=2;
	mx[1][2]=2;
	mx[1][3]=2;
	
	mx[2][0]=1;
	mx[2][1]=2;
	mx[2][2]=3;
	mx[2][3]=3;
	
	mx[3][0]=1;
	mx[3][1]=2;
	mx[3][2]=3;
	mx[3][3]=4;
	tridiagonalize(4,mx,v1,v2);
	for(i=0;i<4;i++){
		for(j=0;j<4;j++){
			printf("%f ",mx[i][j]);
		}
		printf("\n");
	}
}