#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;
}