Changeset 132 for proiecte

Timestamp:

Jan 14, 2010, 12:21:33 AM (14 years ago)

Author:

(none)

Message:

PPPP

File:

• proiecte/PPPP/eigenface2/eigenface_p.c (modified) (5 diffs)

proiecte/PPPP/eigenface2/eigenface_p.c

@@ -16 +16 @@
 a. nrot returns the number of Jacobi rotations that were required.
 */
-void jacobi(float **a, int n, float d[], float **v, int *nrot)
+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;        
+        float tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
+        int tid, nthreads;
+        
         b = (float *) malloc(n * sizeof(float));
         z = (float *) malloc(n * sizeof(float));
         
-        int nthreads, tid;
-
         /* Fork a team of threads giving them their own copies of variables */
         #pragma omp parallel private(nthreads, tid, i, j, ip, iq) shared(tresh, theta, tau, sm, s, h, g, c, b, z, a, n, d, v, nrot)
@@ -41 +40 @@
     }
 
-
-        #pragma omp parallel for
-        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;
-        }
-        
-        #pragma omp parallel for
-        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).
-        }
-        
-        }  /* All threads join master thread and disband */
-        
+        
+        #pragma omp parallel for
+        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;
+        }
+        #pragma omp parallel for
+        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;
-
-                #pragma omp parallel for reduction(+:sm)
-                for (ip=0;ip<n-1;ip++) {                        // Sum off-diagonal elements.
-                        for (iq=ip+1;iq<n;iq++) {
-                                sm += fabs(a[ip][iq]);
-                        //printf("Hello World from thread = %d\n", omp_get_thread_num());
-                        }
-                }
-                if (sm == 0.0) {                                // The normal return, which relies on quadratic convergence to machine underflow.
+                sm=0.0;
+        #pragma omp parallel for reduction(+:sm)
+                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);
@@ -75 +66 @@
                 }
                 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.
+                                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][iq]=0.0;
-                                else if (fabs(a[ip][iq]) > tresh) {
+                                        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][iq])/h;        // t = 1/(2*theta)
+                                                t=(a[ip*n + iq])/h;                                                                                     // t = 1/(2*theta)
                                         else {
-                                                theta=0.5*h/(a[ip][iq]); 
+                                                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;
@@ -98 +88 @@
                                         s=t*c;
                                         tau=s/(1.0+c);
-                                        h=t*a[ip][iq];
+                                        h=t*a[ip*n + iq];
                                         z[ip] -= h;
                                         z[iq] += h;
                                         d[ip] -= h;
                                         d[iq] += h;
-                                        a[ip][iq]=0.0;
-                                
+                                        a[ip*n + iq]=0.0;
+                                        
                                         #pragma omp parallel for
-                                        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)
-                                        } 
+                                        }
                                         #pragma omp parallel for
-                                        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)
                                         }
                                         #pragma omp parallel for
-                                        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)
                                         }
@@ -130 +120 @@
                 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");
-}
-
+                        d[ip]=b[ip];                                                                                                                    // Updte d with the sum of t*a[pq],
+                        z[ip]=0.0;                                                                                                                      // and reinitialize z.
+                }
+        }
+}
 /* 
 Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi