Olumide

#include <stdlib.h>
#include <string.h>
#include "f2c.h"
#include "clapack.h"

// header BLAS routines
double f2c_dnrm2( int *n, double *x, int *incx);
void f2c_dgemm( char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda,
           double *b, int *ldb, double *beta, double *c, int *ldc);


//Wrapper to ease call
double win_dnrm2( int n, double *x, int incx){
    return f2c_dnrm2( &n, x, &incx);
}

void win_dgemm( char transa, char transb, int m, int n, int k, double alpha, double *a, int lda,
               double *b, int ldb, double beta, double *c, int ldc){
    f2c_dgemm( &transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
}

void win_dgesv (int n, int nrhs, double *a, int lda, int *ipiv, double *b, int ldb, int *infoend)
      {
    dgesv_ ( &n, &nrhs, a, &lda, ipiv, b, &ldb, infoend);
      }
void win_dcopy (int n, double *dx, int incx, double *dy, int incy)
      {
        f2c_dcopy ( &n, dx, &incx, dy, &incy);
      }


// MAIN PROGRAM
int main(int argc, char **argv){

    int n, k;
    double *A, *b, *Acopy, *bcopy;
    int *ipiv, info;
    int i;
    double t1,t2,elapsed;
    //struct timeval tp;
    int rtn;
    double normr, normb;
    char notrans;

    // CALL WITHOUT WRAPPERS
    //int one;
    //double d_one, d_mone;
    //int n2;
    //char charN;

    //one =1;
    //d_one=1.0e+0; d_mone=-1.0e+0;
    //charN='N';

    n = 100; k = 1;
    for( i = 1; i < argc; i++ ) {
        if( strcmp( argv[i], "-n" ) == 0 ) {
            n  = atoi(argv[i+1]);
            i++;
        }
        if( strcmp( argv[i], "-k" ) == 0 ) {
            k  = atoi(argv[i+1]);
            i++;
        }
    }

    A = (double *)malloc(n*n*sizeof(double)) ;
    if (A==NULL){ printf("error of memory allocation\n"); exit(0); }
    Acopy = (double *)malloc(n*n*sizeof(double)) ;
    if (Acopy==NULL){ printf("error of memory allocation\n"); exit(0); }
    b = (double *)malloc(n*k*sizeof(double)) ;
    if (b==NULL){ printf("error of memory allocation\n"); exit(0); }
    bcopy = (double *)malloc(n*k*sizeof(double)) ;
    if (bcopy==NULL){ printf("error of memory allocation\n"); exit(0); }
    ipiv = (int *)malloc(n*sizeof(int)) ;
    if (ipiv==NULL){ printf("error of memory allocation\n"); exit(0); }

    for(i=0;i<n*n;i++)
        A[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
    for(i=0;i<n*k;i++)
        b[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;


    // CALL WITHOUT WRAPPERS
    //n2 = n*n;
    //f2c_dcopy(&n,b,&one,bcopy,&one);
    //f2c_dcopy(&n2,A,&one,Acopy,&one);
    win_dcopy(n*k,b,1,bcopy,1);
    win_dcopy(n*n,A,1,Acopy,1);

    //rtn=gettimeofday(&tp, NULL);
    //t1=(double)tp.tv_sec+(1.e-6)*tp.tv_usec;S
    t1=0;

    // CALL WITHOUT WRAPPERS
    //dgesv_( &n, &k, A, &n, ipiv, b, &n, &info);
    win_dgesv( n, k, A, n, ipiv, b, n, &info);
    printf("\tDGESV [info=%d]\n", info);

    //rtn=gettimeofday(&tp, NULL);
    //t2=(double)tp.tv_sec+(1.e-6)*tp.tv_usec;
    t2=1;
    elapsed=t2-t1;

    // CALL WITHOUT WRAPPERS
    //normb = f2c_dnrm2( &n ,bcopy, &one);
    normb = win_dnrm2( n ,bcopy, 1);

    // CALL WITHOUT WRAPPERS
    //f2c_dgemm(&charN,&charN, &n, &k, &n, &d_one, Acopy, &n, b, &n, &d_mone, bcopy, &n);
    win_dgemm('N','N', n, k, n, 1.0e+00, Acopy, n, b, n, -1.0e+00, bcopy, n);

    // CALL WITHOUT WRAPPERS
    //normr = f2c_dnrm2(&n,bcopy,&one);
    normr = win_dnrm2(n,bcopy,1);

    printf("\tDGESV [n=%d,k=%d] \t%fs\t%fMHz\t\t(check:%e)\n",
        n,k,elapsed,
        ((double) n*n*n/3)/elapsed/1.0e6,
        normr/normb);
    free(A);
    free(Acopy);
    free(b);
    free(bcopy);
    free(ipiv);

    return 0;

}