• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# diagonalize

a guest Aug 24th, 2011 136 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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. }
RAW Paste Data
Top