Index: /proiecte/PPPP/eigenface/new/jacobi_vector.c
===================================================================
--- /proiecte/PPPP/eigenface/new/jacobi_vector.c	(revision 93)
+++ /proiecte/PPPP/eigenface/new/jacobi_vector.c	(revision 93)
@@ -0,0 +1,123 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#define ROTATE(a,i,j,k,l) g=a[i*n + j];h=a[k*n + l];a[i*n + j]=g-s*(h+g*tau);\
+a[k*n + 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,t,sm,s,h,g,c,*b,*z;
+	
+	b = (float *) malloc(n * sizeof(float));
+	z = (float *) malloc(n * sizeof(float));
+	
+	for (ip=0;ip<n;ip++) { 															// Initialize to the identity matrix.
+		for (iq=0;iq<n;iq++) v[ip*n + iq]=0.0;
+			v[ip*n + ip]=1.0;
+	}
+	for (ip=0;ip<n;ip++) { 															// Initialize b and d to the diagonal of a.
+		b[ip]=d[ip]=a[ip*n + ip]; 														
+		z[ip]=0.0; 																// This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14).
+	}
+	
+	*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*n + 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*n + 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*n + iq]=0.0;
+				else if (fabs(a[ip*n + iq]) > tresh) {
+					h=d[iq]-d[ip];
+					if ((float)(fabs(h)+g) == (float)fabs(h))
+						t=(a[ip*n + iq])/h; 											// t = 1/(2*theta)
+					else {
+						theta=0.5*h/(a[ip*n + iq]); 									// Equation (11.1.10).
+						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*n + iq];
+					z[ip] -= h;
+					z[iq] += h;
+					d[ip] -= h;
+					d[iq] += h;
+					a[ip*n + 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]; 															// Updte d with the sum of t*a[pq],
+			z[ip]=0.0; 															// and reinitialize z.
+		}
+	}
+	
+	printf("Too many iterations in routine jacobi\n");
+}
+
+/* 
+Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi
+(x11.1) or tqli (x11.3), this routine sorts the eigenvalues into descending order, and rearranges
+the columns of v correspondingly. The method is straight insertion. 
+*/
+void eigsrt(float d[], float **v, int n)	
+{
+	int k,j,i;
+	float p;
+	for (i=0;i<n-1;i++) {
+		p=d[k=i];
+		for (j=i+1;j<n;j++)
+			if (d[j] >= p) p=d[k=j];
+		if (k != i) {
+			d[k]=d[i];
+			d[i]=p;
+			for (j=0;j<n;j++) {
+				p=v[j*n + i];
+				v[j*n + i]=v[j*n + k];
+				v[j*n + k]=p;
+			}
+		}
+	}
+}
