Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

diagonalize

By: a guest on Aug 24th, 2011  |  syntax: C  |  size: 6.67 KB  |  views: 102  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6. //definition for lapack / blas routines
  7. int dgemm_( char *transa, char *transb, const int *m, const int *n, const int *k,
  8.             double *alpha, const double *a, const int *lda, const double *b,
  9.             const int *ldb, double *beta, double *c, const int *ldc);
  10. int dsyev_( char *jobz, char *uplo, int *n, double *a, int *lda, double *w,
  11.             double *work, int *lwork, int *info);
  12.  
  13. /**  [mk_diag diag mat n m] ~ makes a diagonal matrix from an array
  14.  * Places the [m] elements of [diag] along diagonal of [mat] with [n] cols */
  15. void mk_diag(double* diag, double* mat, int n, const int m)
  16. {
  17.     int i;
  18.     for(i =0; i < m; ++i)
  19.         mat[i*n+i] = diag[i];
  20. }
  21.  
  22. int main()
  23. {
  24.     int i, j;
  25.     double m[16] = { //rowmajor
  26.         0, 1, 2, 3,
  27.         1, 2, 6, 4,
  28.         2, 6, 7, 5,
  29.         3, 4, 5, 6 };
  30.     int n = 4;
  31.     int info = 0;
  32.     char jobz = 'V';
  33.     char uplo = 'U';
  34.     int lwork = n * n;
  35.  
  36. //print initial matrix
  37.     printf ("\nInitial Matrix (A):\n\t");
  38.     for (i=0; i<n; ++i) {
  39.         for (j=0; j<n; ++j){
  40.             //printf("[%5.1f] ", m[i][j] );
  41.             printf("[%6.1f] ", m[i*n+j]);
  42.         }
  43.         printf("\n\t");
  44.     }
  45. //find eigenvalues and eigenvectors of a symmetric matrix
  46.     double* x = (double*) malloc( n*sizeof(double));
  47.     double* u = (double*) malloc( n*n*sizeof(double));
  48.     double* work = (double*) malloc( n*n*sizeof(double));
  49.     memcpy (u, m, n*n*sizeof(double));
  50.  
  51.     dsyev_(&jobz,   //include eigenvectors? (V=yes)
  52.            &uplo,   //Upper triangle of A is stored
  53.            &n,      //order of the matrix
  54.            u,       //** input matrix, returns orthonormal eigenvectors
  55.            &n,      //leading dimension of array A
  56.            x,       //** returns eigenvalues (as ROWS)
  57.            work,    //work space
  58.            &lwork,  //length of work; LWORK >= max(1,3*n-1)
  59.            &info);  //**returns info on results
  60.                     //== 0, success;
  61.                     // < 0, -ith argument had issue;
  62.                     // > 0, failed in convergence
  63.     if ( info != 0 ){
  64.         printf("error: %d", info);
  65.         return 0;
  66.     }
  67. //create all the vectors/matrices we need
  68. //  diagonal matrix of eigenvalues...
  69.     double* d = (double*) malloc( n*n*sizeof(double));
  70.     mk_diag( x, d, n, n);
  71.  
  72.     printf("\nD :\n\t");
  73.     for (i=0; i<n; ++i) {
  74.         for (j=0; j<n; ++j)
  75.             printf("[%6.5f] ", d[i*n+j]);
  76.         printf("\n\t");
  77.     }
  78.  
  79.     printf("\nUt :\n\t");
  80.     for (i=0; i<n; ++i) {
  81.         for (j=0; j<n; ++j)
  82.             printf("[%6.5f] ", u[i*n+j]);
  83.         printf("\n\t");
  84.     }
  85.     char _tran = 'T', ntran = 'N';
  86.     double one = 1.0;
  87.     double* I = (double*) malloc( n*n*sizeof(double));
  88.     dgemm_(   &ntran, //const char BLAS_TRANSPOSE  ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
  89.               &_tran, //const char BLAS_TRANSPOSE  ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
  90.                &n,    //const int  M                ~number of rows in A
  91.                &n,    //const int  N                ~number of cols in B
  92.                &n,    //const int  K                ~number of cols in A
  93.                &one,  //const double alpha          ~scalar to multiply after computation
  94.                 u,    //const double *A             ~matrix A
  95.                &n,    //const int lda               ~stride of A (inc for inner loop)
  96.                 u,    //const double *B             ~matrix B
  97.                &n,    //const int ldb               ~stride of B (inc for inner loop)
  98.                &one,  //const double beta           ~scalar to multiply C by
  99.                 I,    //      double *C             ~matrix C -- and output
  100.                &n );  //const int ldc               ~stride of C (inc for inner loop)
  101.  
  102.     printf("\nU * Ut == I :\n\t");
  103.     for (i=0; i<n; ++i) {
  104.         for (j=0; j<n; ++j)
  105.             printf("[%6.5f] ", I[i*n+j]);
  106.         printf("\n\t");
  107.     }
  108.  
  109. //---------------------------------------------------------------------------------------
  110.     double* ud = (double*) malloc( n*n*sizeof(double));
  111.     dgemm_(   &ntran, //const char BLAS_TRANSPOSE  ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
  112.               &ntran, //const char BLAS_TRANSPOSE  ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
  113.                &n,    //const int  M                ~number of rows in A
  114.                &n,    //const int  N                ~number of cols in B
  115.                &n,    //const int  K                ~number of cols in A
  116.                &one,  //const double alpha          ~scalar to multiply after computation
  117.                 u,    //const double *A             ~matrix A
  118.                &n,    //const int lda               ~stride of A (inc for inner loop)
  119.                 d,    //const double *B             ~matrix B
  120.                &n,    //const int ldb               ~stride of B (inc for inner loop)
  121.                &one,  //const double beta           ~scalar to multiply C by
  122.                 ud,   //      double *C             ~matrix C -- and output
  123.                &n );  //const int ldc               ~stride of C (inc for inner loop)
  124.  
  125.     printf("\n U * D :\n\t");
  126.     for (i=0; i<n; ++i) {
  127.         for (j=0; j<n; ++j)
  128.             printf("[%6.5f] ", ud[i*n+j]);
  129.         printf("\n\t");
  130.     }
  131.  
  132.     double* udut = (double*) malloc( n*n*sizeof(double));
  133.     //PERFORMS: C = b*(C)+a*(A*B)
  134.     dgemm_(   &ntran, //const char BLAS_TRANSPOSE  ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
  135.               &_tran, //const char BLAS_TRANSPOSE  ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
  136.                &n,    //const int  M                ~number of rows in A
  137.                &n,    //const int  N                ~number of cols in B
  138.                &n,    //const int  K                ~number of cols in A
  139.                &one,  //const double alpha          ~scalar to multiply after computation
  140.                 ud,    //const double *A             ~matrix A
  141.                &n,    //const int lda               ~stride of A (inc for inner loop)
  142.                 u,    //const double *B             ~matrix B
  143.                &n,    //const int ldb               ~stride of B (inc for inner loop)
  144.                &one,  //const double beta           ~scalar to multiply C by
  145.                udut,  //      double *C             ~matrix C -- and output
  146.                &n );  //const int ldc               ~stride of C (inc for inner loop)
  147.  
  148.     printf("\n(U * D) * U^t == A :\n\t");
  149.     for (i=0; i<n; ++i) {
  150.         for (j=0; j<n; ++j)
  151.             printf("[%6.5f] ", udut[i*n+j]);
  152.         printf("\n\t");
  153.     }
  154.     putchar ('\n');
  155.     return 0;
  156. }