/*
dtrsm_test.c

compilation:
        gcc -o dtrsm_test -Wall -g -O1 -L<path to GotoBLAS> -lgoto -lpthread   dtrsm_test.c
*/

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

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

/*
 * lower triangular random values
 */
void randL(double * vec, int size)
{
   int l,r;
   for( l=0; l<size; l++ )
      for ( r=0; r<l+1; r++)
	 {
	    printf( "idx: %d\n", l*size+r );
	    vec[l*size+r] = drand48()*10;
	 }
}

/*
 * identity matrix
 */
void eye(double * vec, int size)
{
   int l;
   for( l=0; l<size; l++ )
      vec[l*size+l] = 1.0;
}

/*
 * print out matrix (debug purposes)
 */
void matprint(double * vec, int size)
{
        int l,r;
        printf("matrix of size %d x %d:\n", size, size);
	printf("M=[ \n");
        for( l=0; l<size; l++ )
        {
                for ( r=0; r<size; r++)
                {
                        printf("%.12g\t", vec[l*size+r]);
                }
                printf(";\n");
        }
	printf("]\n");
}


int main(int argc, char *argv[] )
{
        double *M;
        double *Mi;
        int n;

        // for BLAS call
        double alpha = 1.0;
        char side = 'L', uplo = 'L', trans = 'N', diag = 'N';

        if (argc < 2)
        {
                printf("please give matrix size as first argument!\n");
                return -1;
        }

        // initialize RNG
        srand48( (unsigned)time((time_t *)NULL) );

        n = atoi(argv[1]);

        // allocate and set memory
        M = (double *) malloc( n*n * sizeof(double) );
        memset(M, 0, n*n*sizeof(double));
        Mi = (double *) malloc( n*n * sizeof(double) );
        memset(Mi, 0, n*n*sizeof(double));

        // create lower triangular and identity matrix
        randL(M, n);
        eye(Mi, n);

        // print out
        matprint(M, n);

	int i,j;

        matprint(Mi, n);
        
        printf("now calling blas.\n");
        dtrsm_( &side , &uplo, &trans, &diag, &n, &n, &alpha, M, &n, Mi, &n );

        printf("done.\nresults:\n");
        // print out
        matprint(Mi, n);

        return 0;
}













/*
//randL
//printf(" idx: %d\n", l*size+r );


	// my print 
	for( i=0; i<n; i++ ){
	   for( j=0; j<n; j++ )
	      printf("M[%d,%d]=%5.3g\t", i, j, M[i+j*n] );	 
	   printf("\n");
	};
*/