// gcc my-trsv.c -o trsm.x -Wall -L/Users/pauldj/works/libs/OpenBLAS-develop/ -lopenblas 
// 
// https://github.com/elmar-peise/ELAPS/tree/master/resources/headers

#include <stdlib.h>
#include <stdio.h> 
#include <time.h>

int dgeqrf_( int *, int *, double *, int *, double *, double *, int *, int * );

int dormqr_( char *, char *, int *, int *, int *, double *, int *,
	     double *, double *, int *, double *, int *, int * );

int dpotrf_(char *, int *, double *, int *, int *);

void dtrmm_(char *, char *, char *, char *, int *, int *, double *,
	   double *, int *, double *, int *);


int main( int argc, char *argv[] ){
   int n, i, j, info; 
   double *M, *A, *tau, *work;
   int nwork;
   
   if(argc == 1){ printf(" argument needed: pass me one argument\n"); return(-1);}
   if(argc > 1)  n = atoi(argv[1]);
   
   M = (double *) malloc( n * n * sizeof(double) );
   
   // M is a square random matrix  
   srand48( (unsigned)time((time_t *)NULL) );
   for( i=0; i<n; i++ )     
      for( j=0; j<n; j++ )  
	 M[ j+i*n ] = drand48();  

   
   // QR factorization of M
   tau   = (double *) malloc( n * sizeof(double) );
   nwork = 50*n;
   work  = (double *) malloc( nwork * sizeof(double) );

   dgeqrf_( &n, &n, M, &n, tau, work, &nwork, &info );
   printf(" info QRfact: %d\n", info );

   
   // intial matrix definition
   A = (double *) malloc( n * n * sizeof(double) );   
   for( j=0; j<n; j++ )
      for( i=0; i<n; i++ )
   	 if( j == i )  A[ j+i*n ] = i+1;
   	 else A[ j+i*n ] = 0;

   
   //  A := Q*A*Q' 
   dormqr_( "L", "N", &n, &n, &n, M, &n, tau, A, &n, work, &nwork, &info );
   printf(" info QRgemm: %d\n", info );

   dormqr_( "R", "T", &n, &n, &n, M, &n, tau, A, &n, work, &nwork, &info );
   printf(" info QRgemm: %d\n", info );

   free(M);
   free(work);
   free(tau);

   // DONE:   Ax = lx   

   
   // B is SPD
   double *B, *L;
   B  = (double *) malloc( n*n * sizeof(double) );
   L  = (double *) malloc( n*n * sizeof(double) );
   double one=1;

   for( i=0; i<n; i++ )     
      for( j=0; j<n; j++ ){
	 if(i==j) B[j+i*n] = drand48() + n;
	 if(i<j)  B[j+i*n] = drand48();
	 if(i>j)  B[j+i*n] = 0;
	 L[j+i*n] = B[j+i*n];
      }	    

   
   // Cholesky:   L*L' = B
   dpotrf_( "L", &n, L, &n, &info );
   printf(" info Chol: %d\n", info );


   // A := L*A*L'
   dtrmm_( "L", "L", "N", "N", &n, &n, &one, L, &n, A, &n );
   dtrmm_( "R", "L", "T", "N", &n, &n, &one, L, &n, A, &n ); 
   free(L);

   
   for( i=0; i<n; i++ )     
      for( j=0; j<n; j++ )
	 if(i>j)  B[j+i*n] = B[i+j*n];

   // DONE:   Ax = lBx   
   
   free(A);
   free(B);   
   return( 0 );
}