Changeset 90
- Timestamp:
- Jan 11, 2010, 4:41:13 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
proiecte/PPPP/eigenface/new/jacobi.c
r89 r90 1 #include <stdio.h> 1 2 #include <math.h> 2 #include "nrutil.h"3 #include <malloc.h> 3 4 #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ 4 5 a[k][l]=h+s*(g-h*tau); … … 13 14 { 14 15 int j,iq,ip,i; 15 float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; 16 float tresh,theta,tau,sm,s,h,g,c,*b,*z; 17 float t; 18 b = (float *) malloc(n * sizeof(float)); 19 z = (float *) malloc(n * sizeof(float)); 16 20 17 // TODO 18 // b=vector(1,n); 19 // z=vector(1,n); 20 21 22 for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. 23 for (iq=0;iq<n;iq++) v[ip][iq]=0.0; 24 v[ip][ip]=1.0; 21 printf("in jacobi\n"); 22 23 for (ip=0;ip<n;ip++) { // Initialize to the identity matrix. 24 for (iq=0;iq<n;iq++) 25 v[ip][iq]=0.0; 26 v[ip][ip]=1.0; 25 27 } 26 for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. 27 b[ip]=d[ip]=a[ip][ip]; 28 z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). 28 printf("in jacobi 2\n"); 29 for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a. 30 b[ip]=d[ip]=a[ip][ip]; 31 z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14). 29 32 } 30 33 34 printf("in jacobi 2\n"); 31 35 *nrot=0; 32 36 for (i=1;i<=50;i++) { 33 37 sm=0.0; 34 for (ip=0;ip<n-1;ip++) { 38 for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements. 35 39 for (iq=ip+1;iq<n;iq++) 36 40 sm += fabs(a[ip][iq]); 37 41 } 38 if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. 39 40 // TODO 41 //free_vector(z,1,n); 42 //free_vector(b,1,n); 42 if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow. 43 free(z); 44 free(b); 43 45 return; 44 46 } 45 47 if (i < 4) 46 tresh=0.2*sm/(n*n); 48 tresh=0.2*sm/(n*n); // ...on the first three sweeps. 47 49 else 48 tresh=0.0; 50 tresh=0.0; // ...thereafter. 49 51 for (ip=0;ip<n-1;ip++) { 50 52 for (iq=ip+1;iq<n;iq++) { 51 53 g=100.0*fabs(a[ip][iq]); 52 if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) 54 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. 53 55 && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) 54 56 a[ip][iq]=0.0; … … 56 58 h=d[iq]-d[ip]; 57 59 if ((float)(fabs(h)+g) == (float)fabs(h)) 58 t=(a[ip][iq])/h; 60 t=(a[ip][iq])/h; // t = 1/(2*theta) 59 61 else { 60 theta=0.5*h/(a[ip][iq]); Equation (11.1.10).62 theta=0.5*h/(a[ip][iq]); 61 63 t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); 62 64 if (theta < 0.0) t = -t; … … 73 75 a[ip][iq]=0.0; 74 76 75 for (j=0;j<ip-1;j++) { 77 for (j=0;j<ip-1;j++) { // Case of rotations 1 <= j < p. 76 78 ROTATE(a,j,ip,j,iq) 77 } 78 for (j=ip+1;j<iq-1;j++) { 79 } 80 for (j=ip+1;j<iq-1;j++) { // Case of rotations p < j < q. 79 81 ROTATE(a,ip,j,j,iq) 80 82 } 81 for (j=iq+1;j<n;j++) { 83 for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n. 82 84 ROTATE(a,ip,j,iq,j) 83 85 } … … 93 95 for (ip=0;ip<n;ip++) { 94 96 b[ip] += z[ip]; 95 d[ip]=b[ip]; 96 z[ip]=0.0; 97 d[ip]=b[ip]; // Update d with the sum of t*a[pq], 98 z[ip]=0.0; // and reinitialize z. 97 99 } 98 100 } … … 100 102 printf("Too many iterations in routine jacobi\n"); 101 103 } 104 105 int main(int argc, char * argv) 106 { 107 int n, i, j; 108 float *d; 109 float **a; 110 float **v; 111 int nrot = 0; 112 113 FILE * f = fopen("mat.in", "r"); 114 fscanf(f,"%d", &n); 115 116 a = (float **)malloc(n * sizeof(float*)); 117 for(i = 0; i < n; i++) 118 a[i] = (float *)malloc(n * sizeof(float)); 119 120 v = (float **)malloc(n * sizeof(float*)); 121 for(i = 0; i < n; i++) 122 v[i] = (float *)malloc(n * sizeof(float)); 123 124 d = (float *)malloc(n * sizeof(float)); 125 126 for(i = 0; i < n; i++) 127 for(j = 0; j < n; j++) 128 fscanf(f,"%f", &a[i][j]); 129 for(i = 0; i < n; i++) 130 { 131 for(j = 0; j < n; j++) 132 printf("%f ", a[i][j]); 133 printf("\n"); 134 } 135 136 jacobi(a, n, d, v, &nrot); 137 138 printf("v:\n"); 139 for(i = 0; i < n; i++) 140 { 141 for(j = 0; j < n; j++) 142 printf("%f ", v[i][j]); 143 printf("\n"); 144 } 145 146 printf("d:\n"); 147 for(i = 0; i < n; i++) 148 { 149 printf("%f ", d[i]); 150 } 151 152 return 0; 153 } 154
Note: See TracChangeset
for help on using the changeset viewer.