Changeset 90


Ignore:
Timestamp:
Jan 11, 2010, 4:41:13 PM (14 years ago)
Author:
(none)
Message:

PPPP - tested jacobi.c

File:
1 edited

Legend:

Unmodified
Added
Removed
  • proiecte/PPPP/eigenface/new/jacobi.c

    r89 r90  
     1#include <stdio.h>
    12#include <math.h>
    2 #include "nrutil.h"
     3#include <malloc.h>
    34#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
    45a[k][l]=h+s*(g-h*tau);
     
    1314{
    1415        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));
    1620       
    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;
    2527        }
    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).
    2932        }
    3033       
     34        printf("in jacobi 2\n");
    3135        *nrot=0;
    3236        for (i=1;i<=50;i++) {
    3337                sm=0.0;
    34                 for (ip=0;ip<n-1;ip++) {                                                                                                        // Sum off-diagonal elements.
     38                for (ip=0;ip<n-1;ip++) {                        // Sum off-diagonal elements.
    3539                        for (iq=ip+1;iq<n;iq++)
    3640                        sm += fabs(a[ip][iq]);
    3741                }
    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);
    4345                        return;
    4446                }
    4547                if (i < 4)
    46                         tresh=0.2*sm/(n*n);                                                                                                             // ...on the first three sweeps.
     48                        tresh=0.2*sm/(n*n);                             // ...on the first three sweeps.
    4749                else
    48                         tresh=0.0;                                                                                                                      // ...thereafter.
     50                        tresh=0.0;                                      // ...thereafter.
    4951                for (ip=0;ip<n-1;ip++) {
    5052                        for (iq=ip+1;iq<n;iq++) {
    5153                                g=100.0*fabs(a[ip][iq]);
    52                                 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.
     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.
    5355                                        && (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
    5456                                        a[ip][iq]=0.0;
     
    5658                                        h=d[iq]-d[ip];
    5759                                        if ((float)(fabs(h)+g) == (float)fabs(h))
    58                                                 t=(a[ip][iq])/h;                                                                                        // t = 1/(2*theta)
     60                                                t=(a[ip][iq])/h;        // t = 1/(2*theta)
    5961                                        else {
    60                                                 theta=0.5*h/(a[ip][iq]); Equation (11.1.10).
     62                                                theta=0.5*h/(a[ip][iq]);
    6163                                                t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
    6264                                                if (theta < 0.0) t = -t;
     
    7375                                        a[ip][iq]=0.0;
    7476                                       
    75                                         for (j=0;j<ip-1;j++) {                                                                                  // Case of rotations 1 <= j < p.
     77                                        for (j=0;j<ip-1;j++) { // Case of rotations 1 <= j < p.
    7678                                                ROTATE(a,j,ip,j,iq)
    77                                         }
    78                                         for (j=ip+1;j<iq-1;j++) {                                                                               // Case of rotations p < j < q.
     79                                        } 
     80                                        for (j=ip+1;j<iq-1;j++) { // Case of rotations p < j < q.
    7981                                                ROTATE(a,ip,j,j,iq)
    8082                                        }
    81                                         for (j=iq+1;j<n;j++) {                                                                                  // Case of rotations q < j <= n.
     83                                        for (j=iq+1;j<n;j++) {  // Case of rotations q < j <= n.
    8284                                                ROTATE(a,ip,j,iq,j)
    8385                                        }
     
    9395                for (ip=0;ip<n;ip++) {
    9496                        b[ip] += z[ip];
    95                         d[ip]=b[ip];                                                                                                                    // Update d with the sum of t*a[pq],
    96                         z[ip]=0.0;                                                                                                                      // and reinitialize z.
     97                        d[ip]=b[ip];            // Update d with the sum of t*a[pq],
     98                        z[ip]=0.0;                      // and reinitialize z.
    9799                }
    98100        }
     
    100102        printf("Too many iterations in routine jacobi\n");
    101103}
     104
     105int 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.