daily pastebin goal
50%
SHARE
TWEET

Olumide

a guest Aug 10th, 2009 279 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdlib.h>
  2. #include <string.h>
  3. #include "f2c.h"
  4. #include "clapack.h"
  5.  
  6. // header BLAS routines
  7. double f2c_dnrm2( int *n, double *x, int *incx);
  8. void f2c_dgemm( char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda,
  9.                    double *b, int *ldb, double *beta, double *c, int *ldc);
  10.  
  11.  
  12. //Wrapper to ease call
  13. double win_dnrm2( int n, double *x, int incx){
  14.         return f2c_dnrm2( &n, x, &incx);
  15. }
  16.  
  17. void win_dgemm( char transa, char transb, int m, int n, int k, double alpha, double *a, int lda,
  18.                            double *b, int ldb, double beta, double *c, int ldc){
  19.         f2c_dgemm( &transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  20. }
  21.  
  22. void win_dgesv (int n, int nrhs, double *a, int lda, int *ipiv, double *b, int ldb, int *infoend)
  23.       {
  24.         dgesv_ ( &n, &nrhs, a, &lda, ipiv, b, &ldb, infoend);
  25.       }
  26. void win_dcopy (int n, double *dx, int incx, double *dy, int incy)
  27.       {
  28.                 f2c_dcopy ( &n, dx, &incx, dy, &incy);
  29.       }
  30.  
  31.  
  32. // MAIN PROGRAM
  33. int main(int argc, char **argv){
  34.  
  35.         int n, k;
  36.         double *A, *b, *Acopy, *bcopy;
  37.         int *ipiv, info;
  38.         int i;
  39.         double t1,t2,elapsed;
  40.         //struct timeval tp;
  41.         int rtn;
  42.         double normr, normb;
  43.     char notrans;
  44.  
  45.         // CALL WITHOUT WRAPPERS
  46.         //int one;
  47.         //double d_one, d_mone;
  48.         //int n2;
  49.     //char charN;
  50.  
  51.         //one =1;
  52.         //d_one=1.0e+0; d_mone=-1.0e+0;
  53.     //charN='N';
  54.  
  55.         n = 100; k = 1;
  56.         for( i = 1; i < argc; i++ ) {
  57.                 if( strcmp( argv[i], "-n" ) == 0 ) {
  58.                         n  = atoi(argv[i+1]);
  59.                         i++;
  60.                 }
  61.                 if( strcmp( argv[i], "-k" ) == 0 ) {
  62.                         k  = atoi(argv[i+1]);
  63.                         i++;
  64.                 }
  65.         }
  66.  
  67.         A = (double *)malloc(n*n*sizeof(double)) ;
  68.         if (A==NULL){ printf("error of memory allocation\n"); exit(0); }
  69.         Acopy = (double *)malloc(n*n*sizeof(double)) ;
  70.         if (Acopy==NULL){ printf("error of memory allocation\n"); exit(0); }
  71.         b = (double *)malloc(n*k*sizeof(double)) ;
  72.         if (b==NULL){ printf("error of memory allocation\n"); exit(0); }
  73.         bcopy = (double *)malloc(n*k*sizeof(double)) ;
  74.         if (bcopy==NULL){ printf("error of memory allocation\n"); exit(0); }
  75.         ipiv = (int *)malloc(n*sizeof(int)) ;
  76.         if (ipiv==NULL){ printf("error of memory allocation\n"); exit(0); }
  77.  
  78.         for(i=0;i<n*n;i++)
  79.                 A[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
  80.         for(i=0;i<n*k;i++)
  81.                 b[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
  82.        
  83.  
  84.         // CALL WITHOUT WRAPPERS
  85.         //n2 = n*n;
  86.         //f2c_dcopy(&n,b,&one,bcopy,&one);
  87.         //f2c_dcopy(&n2,A,&one,Acopy,&one);
  88.         win_dcopy(n*k,b,1,bcopy,1);
  89.         win_dcopy(n*n,A,1,Acopy,1);
  90.  
  91.         //rtn=gettimeofday(&tp, NULL);
  92.         //t1=(double)tp.tv_sec+(1.e-6)*tp.tv_usec;S
  93.         t1=0;
  94.  
  95.         // CALL WITHOUT WRAPPERS
  96.         //dgesv_( &n, &k, A, &n, ipiv, b, &n, &info);
  97.     win_dgesv( n, k, A, n, ipiv, b, n, &info);
  98.         printf("\tDGESV [info=%d]\n", info);
  99.  
  100.         //rtn=gettimeofday(&tp, NULL);
  101.         //t2=(double)tp.tv_sec+(1.e-6)*tp.tv_usec;
  102.         t2=1;
  103.         elapsed=t2-t1;
  104.  
  105.         // CALL WITHOUT WRAPPERS
  106.         //normb = f2c_dnrm2( &n ,bcopy, &one);
  107.         normb = win_dnrm2( n ,bcopy, 1);
  108.  
  109.         // CALL WITHOUT WRAPPERS
  110.         //f2c_dgemm(&charN,&charN, &n, &k, &n, &d_one, Acopy, &n, b, &n, &d_mone, bcopy, &n);
  111.         win_dgemm('N','N', n, k, n, 1.0e+00, Acopy, n, b, n, -1.0e+00, bcopy, n);
  112.  
  113.     // CALL WITHOUT WRAPPERS
  114.         //normr = f2c_dnrm2(&n,bcopy,&one);
  115.         normr = win_dnrm2(n,bcopy,1);
  116.                
  117.         printf("\tDGESV [n=%d,k=%d] \t%fs\t%fMHz\t\t(check:%e)\n",
  118.                 n,k,elapsed,
  119.                 ((double) n*n*n/3)/elapsed/1.0e6,
  120.                 normr/normb);
  121.         free(A);
  122.         free(Acopy);
  123.         free(b);
  124.         free(bcopy);
  125.         free(ipiv);
  126.  
  127.         return 0;
  128.  
  129. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top