Changeset 132
- Timestamp:
- Jan 14, 2010, 12:21:33 AM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
proiecte/PPPP/eigenface2/eigenface_p.c
r115 r132 16 16 a. nrot returns the number of Jacobi rotations that were required. 17 17 */ 18 void jacobi(float * *a, int n, float d[], float **v, int *nrot)18 void jacobi(float *a, int n, float d[], float *v, int *nrot) 19 19 { 20 20 int j,iq,ip,i; 21 float tresh,theta,tau,sm,s,h,g,c,*b,*z; 22 float t; 21 float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; 22 int tid, nthreads; 23 23 24 b = (float *) malloc(n * sizeof(float)); 24 25 z = (float *) malloc(n * sizeof(float)); 25 26 26 int nthreads, tid;27 28 27 /* Fork a team of threads giving them their own copies of variables */ 29 28 #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 } 42 41 43 44 #pragma omp parallel for 45 for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. 46 for (iq=0;iq<n;iq++) 47 v[ip][iq]=0.0; 48 v[ip][ip]=1.0; 49 } 50 51 #pragma omp parallel for 52 for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. 53 b[ip]=d[ip]=a[ip][ip]; 54 z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). 55 } 56 57 } /* All threads join master thread and disband */ 58 42 43 #pragma omp parallel for 44 for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. 45 for (iq=0;iq<n;iq++) v[ip*n + iq]=0.0; 46 v[ip*n + ip]=1.0; 47 } 48 #pragma omp parallel for 49 for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. 50 b[ip]=d[ip]=a[ip*n + ip]; 51 z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). 52 } 53 } 59 54 *nrot=0; 60 55 for (i=1;i<=50;i++) { 61 62 sm = 0.0; 63 64 #pragma omp parallel for reduction(+:sm) 65 for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements. 66 for (iq=ip+1;iq<n;iq++) { 67 sm += fabs(a[ip][iq]); 68 //printf("Hello World from thread = %d\n", omp_get_thread_num()); 69 } 70 } 71 if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. 56 sm=0.0; 57 #pragma omp parallel for reduction(+:sm) 58 for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements. 59 for (iq=ip+1;iq<n;iq++) 60 sm += fabs(a[ip*n + iq]); 61 } 62 if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. 72 63 free(z); 73 64 free(b); … … 75 66 } 76 67 if (i < 4) 77 tresh=0.2*sm/(n*n); // ...on the first three sweeps.68 tresh=0.2*sm/(n*n); // ...on the first three sweeps. 78 69 else 79 tresh=0.0; // ...thereafter. 80 70 tresh=0.0; // ...thereafter. 81 71 for (ip=0;ip<n-1;ip++) { 82 72 for (iq=ip+1;iq<n;iq++) { 83 g=100.0*fabs(a[ip ][iq]);84 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.73 g=100.0*fabs(a[ip*n + iq]); 74 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. 85 75 && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) 86 a[ip ][iq]=0.0;87 else if (fabs(a[ip ][iq]) > tresh) {76 a[ip*n + iq]=0.0; 77 else if (fabs(a[ip*n + iq]) > tresh) { 88 78 h=d[iq]-d[ip]; 89 79 if ((float)(fabs(h)+g) == (float)fabs(h)) 90 t=(a[ip ][iq])/h;// t = 1/(2*theta)80 t=(a[ip*n + iq])/h; // t = 1/(2*theta) 91 81 else { 92 theta=0.5*h/(a[ip ][iq]);82 theta=0.5*h/(a[ip*n + iq]); // Equation (11.1.10). 93 83 t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); 94 84 if (theta < 0.0) t = -t; … … 98 88 s=t*c; 99 89 tau=s/(1.0+c); 100 h=t*a[ip ][iq];90 h=t*a[ip*n + iq]; 101 91 z[ip] -= h; 102 92 z[iq] += h; 103 93 d[ip] -= h; 104 94 d[iq] += h; 105 a[ip ][iq]=0.0;106 95 a[ip*n + iq]=0.0; 96 107 97 #pragma omp parallel for 108 for (j=0;j<=ip-1;j++) { // Case of rotations 1 <= j < p.98 for (j=0;j<=ip-1;j++) { // Case of rotations 1 <= j < p. 109 99 ROTATE(a,j,ip,j,iq) 110 } 100 } 111 101 #pragma omp parallel for 112 for (j=ip+1;j<=iq-1;j++) { // Case of rotations p < j < q.102 for (j=ip+1;j<=iq-1;j++) { // Case of rotations p < j < q. 113 103 ROTATE(a,ip,j,j,iq) 114 104 } 115 105 #pragma omp parallel for 116 for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n.106 for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n. 117 107 ROTATE(a,ip,j,iq,j) 118 108 } … … 130 120 for (ip=0;ip<n;ip++) { 131 121 b[ip] += z[ip]; 132 d[ip]=b[ip]; // Update d with the sum of t*a[pq], 133 z[ip]=0.0; // and reinitialize z. 134 } 135 } 136 137 printf("Too many iterations in routine jacobi\n"); 138 } 139 122 d[ip]=b[ip]; // Updte d with the sum of t*a[pq], 123 z[ip]=0.0; // and reinitialize z. 124 } 125 } 126 } 140 127 /* 141 128 Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi
Note: See TracChangeset
for help on using the changeset viewer.