Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include <math.h>
- //definition for lapack / blas routines
- int dgemm_( char *transa, char *transb, const int *m, const int *n, const int *k,
- double *alpha, const double *a, const int *lda, const double *b,
- const int *ldb, double *beta, double *c, const int *ldc);
- int dsyev_( char *jobz, char *uplo, int *n, double *a, int *lda, double *w,
- double *work, int *lwork, int *info);
- /** [mk_diag diag mat n m] ~ makes a diagonal matrix from an array
- * Places the [m] elements of [diag] along diagonal of [mat] with [n] cols */
- void mk_diag(double* diag, double* mat, int n, const int m)
- {
- int i;
- for(i =0; i < m; ++i)
- mat[i*n+i] = diag[i];
- }
- int main()
- {
- int i, j;
- double m[16] = { //rowmajor
- 0, 1, 2, 3,
- 1, 2, 6, 4,
- 2, 6, 7, 5,
- 3, 4, 5, 6 };
- int n = 4;
- int info = 0;
- char jobz = 'V';
- char uplo = 'U';
- int lwork = n * n;
- //print initial matrix
- printf ("\nInitial Matrix (A):\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j){
- //printf("[%5.1f] ", m[i][j] );
- printf("[%6.1f] ", m[i*n+j]);
- }
- printf("\n\t");
- }
- //find eigenvalues and eigenvectors of a symmetric matrix
- double* x = (double*) malloc( n*sizeof(double));
- double* u = (double*) malloc( n*n*sizeof(double));
- double* work = (double*) malloc( n*n*sizeof(double));
- memcpy (u, m, n*n*sizeof(double));
- dsyev_(&jobz, //include eigenvectors? (V=yes)
- &uplo, //Upper triangle of A is stored
- &n, //order of the matrix
- u, //** input matrix, returns orthonormal eigenvectors
- &n, //leading dimension of array A
- x, //** returns eigenvalues (as ROWS)
- work, //work space
- &lwork, //length of work; LWORK >= max(1,3*n-1)
- &info); //**returns info on results
- //== 0, success;
- // < 0, -ith argument had issue;
- // > 0, failed in convergence
- if ( info != 0 ){
- printf("error: %d", info);
- return 0;
- }
- //create all the vectors/matrices we need
- // diagonal matrix of eigenvalues...
- double* d = (double*) malloc( n*n*sizeof(double));
- mk_diag( x, d, n, n);
- printf("\nD :\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j)
- printf("[%6.5f] ", d[i*n+j]);
- printf("\n\t");
- }
- printf("\nUt :\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j)
- printf("[%6.5f] ", u[i*n+j]);
- printf("\n\t");
- }
- char _tran = 'T', ntran = 'N';
- double one = 1.0;
- double* I = (double*) malloc( n*n*sizeof(double));
- dgemm_( &ntran, //const char BLAS_TRANSPOSE ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &_tran, //const char BLAS_TRANSPOSE ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &n, //const int M ~number of rows in A
- &n, //const int N ~number of cols in B
- &n, //const int K ~number of cols in A
- &one, //const double alpha ~scalar to multiply after computation
- u, //const double *A ~matrix A
- &n, //const int lda ~stride of A (inc for inner loop)
- u, //const double *B ~matrix B
- &n, //const int ldb ~stride of B (inc for inner loop)
- &one, //const double beta ~scalar to multiply C by
- I, // double *C ~matrix C -- and output
- &n ); //const int ldc ~stride of C (inc for inner loop)
- printf("\nU * Ut == I :\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j)
- printf("[%6.5f] ", I[i*n+j]);
- printf("\n\t");
- }
- //---------------------------------------------------------------------------------------
- double* ud = (double*) malloc( n*n*sizeof(double));
- dgemm_( &ntran, //const char BLAS_TRANSPOSE ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &ntran, //const char BLAS_TRANSPOSE ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &n, //const int M ~number of rows in A
- &n, //const int N ~number of cols in B
- &n, //const int K ~number of cols in A
- &one, //const double alpha ~scalar to multiply after computation
- u, //const double *A ~matrix A
- &n, //const int lda ~stride of A (inc for inner loop)
- d, //const double *B ~matrix B
- &n, //const int ldb ~stride of B (inc for inner loop)
- &one, //const double beta ~scalar to multiply C by
- ud, // double *C ~matrix C -- and output
- &n ); //const int ldc ~stride of C (inc for inner loop)
- printf("\n U * D :\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j)
- printf("[%6.5f] ", ud[i*n+j]);
- printf("\n\t");
- }
- double* udut = (double*) malloc( n*n*sizeof(double));
- //PERFORMS: C = b*(C)+a*(A*B)
- dgemm_( &ntran, //const char BLAS_TRANSPOSE ~(op) on A (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &_tran, //const char BLAS_TRANSPOSE ~(op) on B (Cblas[Trans/NoTrans/Conj/ConjTrans])
- &n, //const int M ~number of rows in A
- &n, //const int N ~number of cols in B
- &n, //const int K ~number of cols in A
- &one, //const double alpha ~scalar to multiply after computation
- ud, //const double *A ~matrix A
- &n, //const int lda ~stride of A (inc for inner loop)
- u, //const double *B ~matrix B
- &n, //const int ldb ~stride of B (inc for inner loop)
- &one, //const double beta ~scalar to multiply C by
- udut, // double *C ~matrix C -- and output
- &n ); //const int ldc ~stride of C (inc for inner loop)
- printf("\n(U * D) * U^t == A :\n\t");
- for (i=0; i<n; ++i) {
- for (j=0; j<n; ++j)
- printf("[%6.5f] ", udut[i*n+j]);
- printf("\n\t");
- }
- putchar ('\n');
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement