#include "real.h"
#include "mat.h"
real abs(real x){if(x>0) return x; else return -x;}

void qrdiag(real **a,int n){ int i,j; real s, **r, **q, eps=0.0001;

q = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) q[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));

for(i=0;i<n;i++) for(j=0;j<n;j++) r[i][j]=0;

while(n>1){
s=a[n-1][n-1];
for(i=0;i<n;i++) for(j=0;j<n;j++) q[i][j]=a[i][j];
for(i=0;i<n;i++) q[i][i]-=s;
QR_dec(q,n,n,r); mat_mat(r,n,n,q,n,a);
for(i=0;i<n;i++) a[i][i]+=s;
printf("matrix T:\n"); print_mat(a,n,n);
s=0;for(i=0;i<n-1;i++)s+=abs(a[i][n-1]);
if(s<eps) n--;
		}
free(r); free(q);
}