SVD (singular value decomposition) algorithm _ compute arbitrary n*m matrix _c language code __ algorithm

Help students solve a problem.

SVD (singular value decomposition algorithm) of the C language code can be found on the Internet, such as the "numerical Recipes in C," a book gives the code: http://cacs.usc.edu/education/phys516/src/TB/ SVDCMP.C, but also gives the demo code HTTP://CACS.USC.EDU/EDUCATION/PHYS516/SRC/TB/SINGULAR.C. However, the demo code can only compute the n*n matrix, the calculation of m*n (that is, M, N unequal) matrix will be error.

I have modified the SINGULAR.C, svdcmp.c unmodified, so that I can compute any m*n matrix. The main finding is a workaround, which is commented in the following code. Because the study of this algorithm is not in-depth, can not be guaranteed to be correct, but tested 6 sets of data, the operation results of this program and the results of the operation of MATLAB is consistent.

In the file singular.c, 2 sets of test data are given, the test data is the 3x3 matrix, which can change the M and N values to become 2x3 matrix at runtime. And the program at the end of the test results, whether the correct look at the know.

The two files are placed in the same directory and can be compiled via the GCC command in Linux:

gcc./singular.c./svdcmp.c-o./singular

When the compilation is complete, execution./singular is ready to run.

File singular.c:

    
    
     
     #include <stdio.h> #include <stdlib.h>   #define M
     
     2 #define N 3//input matrix size #if (m>n)//Take the larger number of M, n #define MN m #else #define MN N #endif   double **dmatrix (int, int, int, in
     
     T);
     
     Double *dvector (int, int);
     
     void svdcmp (double * *, int, int, double *, double * *);
     
         void Print_r (double **a, int m, int n) {int i, J; for (i = 1; I <= m. i++) {for (j = 1; J <= N; j) {printf ("%le", a[i][
     
             J]);
     
         printf ("\ n");
     
          }   int Main () {* * * is declared using * *, and then uses Dmatrix to request space
     
         * So when the function references the matrix, you do not have to declare the matrix size/double **a;
     
 Double *w;        Double **u,**v;
     
         int i,j,k;
     
         Double T;
     
     Double T1[MN],T2[MN];
     
          /* matrices apply for space with the maximum value in M,n to avoid cross-border/a = Dmatrix (1, MN, 1, MN);
     
         U = Dmatrix (1, MN, 1, MN);
     
         W = dvector (1, MN);
     
     v = Dmatrix (1, MN, 1, MN);
     
           A[1][1] = 3.0;
     
         A[1][2] = 2.0;
     
         A[1][3] = 1.0;
     
         A[2][1] = 1.0;
     
         A[2][2] = 3.0;
     
         A[2][3] = 2.0;
     
         A[3][1] = 1.0;
     
         A[3][2] = 2.0;
     
     A[3][3] = 3.0;
     
            /* Alternative Authentication data *//A[1][1] = 1.0;
     
         A[1][2] = 2.0;
     
         A[1][3] = 3.0;
     
         A[2][1] = 4.0;
     
         A[2][2] = 5.0;
     
         A[2][3] = 6.0;
     
         A[3][1] = 7.0;
     
         A[3][2] = 8.0;
     
     A[3][3] = 9.0;  
     
           printf ("= = = a ===\n");
     
      Print_r (A, M, N);