#include <stdio.h>
#include <math.h>
#include <malloc.h>
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
a[k][l]=h+s*(g-h*tau);

/*
Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On
output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a.
v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of
a. nrot returns the number of Jacobi rotations that were required.
*/
void jacobi(float **a, int n, float d[], float **v, int *nrot)
{
	int j,iq,ip,i;
	float tresh,theta,tau,sm,s,h,g,c,*b,*z;
	float t;	
	b = (float *) malloc(n * sizeof(float));
	z = (float *) malloc(n * sizeof(float));
	
	printf("in jacobi\n");
		
	for (ip=0;ip<n;ip++) { 		// Initialize to the identity matrix.
		for (iq=0;iq<n;iq++) 
			v[ip][iq]=0.0;
		v[ip][ip]=1.0;
	}
	printf("in jacobi 2\n");
	for (ip=0;ip<n;ip++) {     // Initialize b and d to the diagonal of a.
		b[ip]=d[ip]=a[ip][ip]; 					
		z[ip]=0.0; 			   // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14).
	}
	
	printf("in jacobi 2\n");
	*nrot=0;
	for (i=1;i<=50;i++) {
		sm=0.0;
		for (ip=0;ip<n-1;ip++) { 			// Sum off-diagonal elements.
			for (iq=ip+1;iq<n;iq++)
			sm += fabs(a[ip][iq]);
		}
		if (sm == 0.0) { 		     		// The normal return, which relies on quadratic convergence to machine underflow.
			free(z);
			free(b);
			return;
		}
		if (i < 4)
			tresh=0.2*sm/(n*n);				// ...on the first three sweeps.
		else
			tresh=0.0;		     			// ...thereafter.
		for (ip=0;ip<n-1;ip++) {
			for (iq=ip+1;iq<n;iq++) {
				g=100.0*fabs(a[ip][iq]);
				if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip])	// After four sweeps, skip the rotation if the off-diagonal element is small.
					&& (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
					a[ip][iq]=0.0;
				else if (fabs(a[ip][iq]) > tresh) {
					h=d[iq]-d[ip];
					if ((float)(fabs(h)+g) == (float)fabs(h))
						t=(a[ip][iq])/h; 	// t = 1/(2*theta)
					else {
						theta=0.5*h/(a[ip][iq]); 
						t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
						if (theta < 0.0) t = -t;
					}
					
					c=1.0/sqrt(1+t*t);
					s=t*c;
					tau=s/(1.0+c);
					h=t*a[ip][iq];
					z[ip] -= h;
					z[iq] += h;
					d[ip] -= h;
					d[iq] += h;
					a[ip][iq]=0.0;
					
					for (j=0;j<ip-1;j++) { // Case of rotations 1 <= j < p.
						ROTATE(a,j,ip,j,iq)
					} 
					for (j=ip+1;j<iq-1;j++) { // Case of rotations p < j < q.
						ROTATE(a,ip,j,j,iq)
					}
					for (j=iq+1;j<n;j++) { 	// Case of rotations q < j <= n.
						ROTATE(a,ip,j,iq,j)
					}
					for (j=0;j<n;j++) {
						ROTATE(v,j,ip,j,iq)
					}
					
					++(*nrot);
				}
			}
		}
		
		for (ip=0;ip<n;ip++) {
			b[ip] += z[ip];
			d[ip]=b[ip]; 		// Update d with the sum of t*a[pq],
			z[ip]=0.0; 			// and reinitialize z.
		}
	}
	
	printf("Too many iterations in routine jacobi\n");
}

int main(int argc, char * argv)
{
	int n, i, j;
	float *d;
	float **a;
	float **v;
	int nrot = 0;
	
	FILE * f = fopen("mat.in", "r");
	fscanf(f,"%d", &n);
	
	a = (float **)malloc(n * sizeof(float*));
	for(i = 0; i < n; i++)
		a[i] = (float *)malloc(n * sizeof(float));
	
	v  = (float **)malloc(n * sizeof(float*));
	for(i = 0; i < n; i++)
		v[i] = (float *)malloc(n * sizeof(float));
	
	d = (float *)malloc(n * sizeof(float));
	
	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
			fscanf(f,"%f", &a[i][j]);
	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
			printf("%f ", a[i][j]);
		printf("\n");
	}
	
	jacobi(a, n, d, v, &nrot); 
	
	printf("v:\n");
	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
			printf("%f ", v[i][j]);
		printf("\n");
	}

	printf("d:\n");
	for(i = 0; i < n; i++)
	{
		printf("%f ", d[i]);
	}

	return 0;
}