Index: /proiecte/PPPP/eigenface/new/jacobi.c
===================================================================
--- /proiecte/PPPP/eigenface/new/jacobi.c	(revision 89)
+++ /proiecte/PPPP/eigenface/new/jacobi.c	(revision 90)
@@ -1,4 +1,5 @@
+#include <stdio.h>
 #include <math.h>
-#include "nrutil.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);
@@ -13,42 +14,43 @@
 {
 	int j,iq,ip,i;
-	float tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
+	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));
 	
-	// TODO
-	// b=vector(1,n);
-	// z=vector(1,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\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;
 	}
-	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");
+	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 (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.
-			
-			// TODO
-			//free_vector(z,1,n);
-			//free_vector(b,1,n);
+		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.
+			tresh=0.2*sm/(n*n);				// ...on the first three sweeps.
 		else
-			tresh=0.0;															// ...thereafter.
+			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.
+				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;
@@ -56,7 +58,7 @@
 					h=d[iq]-d[ip];
 					if ((float)(fabs(h)+g) == (float)fabs(h))
-						t=(a[ip][iq])/h; 											// t = 1/(2*theta)
+						t=(a[ip][iq])/h; 	// t = 1/(2*theta)
 					else {
-						theta=0.5*h/(a[ip][iq]); Equation (11.1.10).
+						theta=0.5*h/(a[ip][iq]); 
 						t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
 						if (theta < 0.0) t = -t;
@@ -73,11 +75,11 @@
 					a[ip][iq]=0.0;
 					
-					for (j=0;j<ip-1;j++) { 											// Case of rotations 1 <= j < p.
+					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.
+					} 
+					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.
+					for (j=iq+1;j<n;j++) { 	// Case of rotations q < j <= n.
 						ROTATE(a,ip,j,iq,j)
 					}
@@ -93,6 +95,6 @@
 		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.
+			d[ip]=b[ip]; 		// Update d with the sum of t*a[pq],
+			z[ip]=0.0; 			// and reinitialize z.
 		}
 	}
@@ -100,2 +102,53 @@
 	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;
+}
+
