Statistiques
| Révision :

root / BLAS / xGEMM / xGEMM.c @ 250

Historique | Voir | Annoter | Télécharger (21,4 ko)

1
/* 
2
   Performs matrix multiply on several BLAS implementations 
3
   Copyleft Emmanuel QUEMENER <emmanuel.quemener@gmail.com> under GPLv3
4

5
   2014-03-14 : Add clBLAS implementation
6

7
   Thanks for help from aurel32@debian.org
8
*/
9

    
10
#include <stdio.h>
11
#include <math.h>
12
#include <stdlib.h>
13
#include <sys/time.h>
14
#include <string.h>
15

    
16
#ifdef CLBLAS
17
#include <clBLAS.h>
18
/* Precise here to avoid new specific bench function */
19
int MyPlatform;
20
int MyDevice;
21
#elif CUBLAS
22
#include <cublas.h>
23
#define CUBLAS_WRAPPER_ERROR_NOERR      0
24
#define CUBLAS_WRAPPER_ERROR_ALLOC      1
25
#define CUBLAS_WRAPPER_ERROR_SET        2
26
#define CUBLAS_WRAPPER_ERROR_GET        3
27
#define CUBLAS_WRAPPER_ERROR_STUB       4
28
#elif THUNKING
29
#include <cublas.h>
30
#include "fortran_common.h"
31
#include "fortran_thunking.h"
32
#elif FBLAS
33
#include <f77blas.h>
34
#elif GSL
35
#include <gsl_cblas.h>
36
#elif ACML
37
#include <acml.h>
38
#else
39
#include <cblas.h>
40
//#include <blaswrap.h>
41
#endif
42

    
43
#ifdef CLBLAS
44

    
45
#ifdef FP64
46
#define LENGTH cl_double
47
#else
48
#define LENGTH cl_float
49
#endif
50

    
51
#else
52

    
53
#ifdef FP64
54
#define LENGTH double
55
#else
56
#define LENGTH float
57
#endif
58

    
59
#endif
60

    
61
/* #ifdef FBLAS */
62

    
63
/* #ifdef FP64 */
64

    
65
/* void F77_dgemm(FCHAR, FCHAR, FINT, FINT, FINT, const double *, const double *, FINT,  */
66
/*                const double *, FINT, const double *, double *, FINT); */
67

    
68
/* #else */
69

    
70
/* void F77_sgemm(FCHAR, FCHAR, FINT, FINT, FINT, const float *, const float *, FINT,  */
71
/*                const float *, FINT, const float *, float *, FINT); */
72

    
73
/* #endif */
74
/* #endif */
75

    
76
/* Matrix with only defined triangular terms */
77
/* Even if there are 0 in matrix, must be defined at all ! */
78

    
79
/* Get from fortran.c */
80

    
81
#ifdef CUBLAS
82
static char *errMsg[5] = 
83
{
84
    "no error",
85
    "allocation error",
86
    "setVector/setMatrix error",
87
    "getVector/getMatrix error",
88
    "not implemented"
89
};
90

    
91
static void wrapperError (const char *funcName, int error)
92
{
93
    printf ("cublas%s wrapper: %s\n", funcName, errMsg[error]);
94
    fflush (stdout);
95
}
96
#endif
97

    
98
int printVector(const int dimVector,const LENGTH *dataVector,
99
                char *nameVector,char *mesgVector)
100
{
101
#ifndef QUIET
102

    
103
  int i;
104
  printf("\n%s of %s, size %i:\n",mesgVector,nameVector,dimVector);
105
  for (i=0;i<dimVector;i++)
106
    {
107
      printf("%s[%i]=%2.10e\n",nameVector,i,dataVector[i]);
108
    }
109
#endif
110

    
111
  return 0;
112
}
113
  
114
int printResults(const int dimVector,const LENGTH *dataVector,
115
                 char *nameVector,char *mesgVector)
116
{
117
#ifdef RESULTS
118
  int i;
119

    
120
  printf("\n%s of %s, size %i:\n",mesgVector,nameVector,dimVector);
121
  for (i=0;i<dimVector;i++)
122
    {
123
      printf("%s[%i]=%2.10e\n",nameVector,i,dataVector[i]);
124
    }
125
#endif
126
  return 0;
127
}
128
  
129
#ifdef CUBLAS
130
int printVectorGPU(const int dimVector,const LENGTH *dataVector,
131
                   char *nameVector,char *mesgVector)
132
{
133
#ifndef QUIET
134
  int i;
135
  cublasStatus stat;
136
  LENGTH *P=0;
137
  int incx=1;
138

    
139
  P=malloc(dimVector*sizeof(LENGTH));
140
  
141
  stat=cublasGetVector(dimVector,sizeof(P[0]),dataVector,incx,P,incx);
142

    
143
  if (stat != CUBLAS_STATUS_SUCCESS) {
144
    wrapperError ("ToGet", CUBLAS_WRAPPER_ERROR_GET);
145
  }  
146

    
147
  printf("\n%s of %s, size %i:\n",mesgVector,nameVector,dimVector);
148
  for (i=0;i<dimVector;i++)
149
    {
150
      printf("%s[%i]=%2.10e\n",nameVector,i,P[i]);
151
    }
152

    
153
  free(P);  
154
#endif
155

    
156
  return 0;
157
}
158
#endif
159

    
160
int bench(int dim,int RUNS)
161
{
162
  /*
163
  int dim=1000;
164
  int RUNS=100;
165
  int incx=1;
166
  */
167
#ifdef PRINT
168
  LENGTH factor=1.;
169
#endif
170

    
171
  LENGTH alpha=1.,beta=0.;
172
  LENGTH *A,*B,*C,*D;
173

    
174
  /* checkBefore checkAfter checks */
175
  LENGTH *checksA,*checksB;
176

    
177
  int i=0, j=0;
178

    
179
  double duration;
180

    
181
  struct timeval tv1,tv2;
182
  struct timezone tz;
183

    
184
  /* Create 4 Matrix of dimension dim by dim  */
185

    
186
  A=malloc(dim*dim*sizeof(LENGTH));
187
  B=malloc(dim*dim*sizeof(LENGTH));
188
  C=malloc(dim*dim*sizeof(LENGTH));
189
  D=malloc(dim*dim*sizeof(LENGTH));
190

    
191
  /* Create 2 vectors for checker Before and After */
192

    
193
  checksA=malloc(RUNS*sizeof(LENGTH));
194
  checksB=malloc(RUNS*sizeof(LENGTH));
195

    
196
  /* Initialize elements with random numbers */
197
  /* Initialize the seed for rand() */
198
  /* srand(time()); */
199

    
200
  for (i=0; i<dim; i++) {
201
    for (j=0; j<dim; j++) {
202
        A[i*dim+j]=(LENGTH)rand()/(RAND_MAX+1.)
203
          *(LENGTH)(i+1.)/(LENGTH)(j+1.); 
204
        B[i*dim+j]=(LENGTH)rand()/(RAND_MAX+1.)
205
          *(LENGTH)(i+1.)/(LENGTH)(j+1.);
206
        C[i*dim+j]=0.;
207
        D[i*dim+j]=0.;
208
    }
209
  }
210
  /*
211
  A[0]=1;
212
  A[1]=2;
213
  A[2]=3;
214
  A[3]=4;
215
  
216
  B[0]=5;
217
  B[1]=6;
218
  B[2]=7;
219
  B[3]=8;
220
  */
221

    
222
  /* Print the matrix */
223
 
224
#ifdef QUIET
225
#else
226
  for (i=0; i<dim; i++) {
227
    for (j=0; j<dim; j++) printf("A[%i,%i]=%1.5f ", i,j,A[i*dim+j]);
228
    putchar('\n');
229
  }
230
  putchar('\n');
231
  for (i=0; i<dim; i++) {
232
    for (j=0; j<dim; j++) printf("B[%i,%i]=%1.5f ", i,j,B[i*dim+j]);
233
    putchar('\n');
234
  }
235
  putchar('\n');
236
#endif
237

    
238
 /* Get first timer before launching */
239
  gettimeofday(&tv1, &tz);
240

    
241
  /* Compute with CLBLAS library  */
242
#ifdef CLBLAS
243

    
244
  cl_uint platformCount;
245
  cl_platform_id* platforms;
246
  cl_uint deviceCount;
247
  cl_device_id* devices;
248

    
249
  cl_int err,errA,errB,errC,errD;
250
  cl_context_properties props[3] = { CL_CONTEXT_PLATFORM, 0, 0 };
251
  cl_context ctx = 0;
252
  cl_command_queue queue = 0;
253
  cl_mem bufA, bufB, bufC, bufD;
254
  cl_event event = NULL;
255

    
256
  char* value;
257
  size_t valueSize;
258

    
259
  // tv3 Put on Device: Allocate & Write buffer
260
  // tv4 Compute
261
  struct timeval tv3,tv4;
262

    
263
  printf("Using CLBLAS: %i iterations for %ix%i matrix on (%d,%d)\n",
264
         RUNS,dim,dim,MyPlatform,MyDevice);
265

    
266
  /* Setup OpenCL environment. */
267
  /* - get all platforms and select MyPlatform */
268
  /* - get all devices from MyPlatform and select MyDevice */
269

    
270
  // Get all platforms
271
  err = clGetPlatformIDs(0, NULL, &platformCount);
272
  platforms = (cl_platform_id*) malloc(sizeof(cl_platform_id) * platformCount);
273
  err = clGetPlatformIDs(platformCount, platforms, NULL);
274

    
275
  // Get Device defined
276
  err = clGetDeviceIDs(platforms[MyPlatform], CL_DEVICE_TYPE_ALL, 0, NULL, &deviceCount);
277
  devices = (cl_device_id*) malloc(sizeof(cl_device_id) * deviceCount);
278
  err = clGetDeviceIDs(platforms[MyPlatform], CL_DEVICE_TYPE_ALL, deviceCount, devices, NULL);  
279

    
280
  // print device name
281
  err = clGetDeviceInfo(devices[MyDevice], CL_DEVICE_NAME, 0, NULL, &valueSize);
282
  value = (char*) malloc(valueSize);
283
  err = clGetDeviceInfo(devices[MyDevice], CL_DEVICE_NAME, valueSize, value, NULL);
284
  printf("Device (%d,%d): %s\n",MyPlatform,MyDevice, value);
285
  free(value);
286

    
287
  props[1] = (cl_context_properties)platforms[MyPlatform];
288

    
289
  /* Initialize Context */
290
  ctx = clCreateContext( props, 1, &devices[MyDevice], NULL, NULL, &err );
291
  queue = clCreateCommandQueue( ctx, devices[MyDevice], 0, &err );
292

    
293
  /* Setup clBLAS */
294
  err = clblasSetup( );
295

    
296
  /* Prepare OpenCL memory objects and place matrices inside them. */
297
  bufA = clCreateBuffer(ctx,CL_MEM_READ_ONLY,dim*dim*sizeof(*A),NULL,&errA );
298
  bufB = clCreateBuffer(ctx,CL_MEM_READ_ONLY,dim*dim*sizeof(*B),NULL,&errB );
299
  bufC = clCreateBuffer(ctx,CL_MEM_READ_WRITE,dim*dim*sizeof(*C),NULL,&errC );
300
  bufD = clCreateBuffer(ctx,CL_MEM_READ_WRITE,dim*dim*sizeof(*D),NULL,&errD );
301

    
302
  errA = clEnqueueWriteBuffer( queue,bufA,CL_TRUE,0,
303
                               dim*dim*sizeof(*A),A,0,NULL,NULL );
304
  errB = clEnqueueWriteBuffer( queue, bufB, CL_TRUE,0,
305
                               dim*dim*sizeof(*B),B,0,NULL,NULL );
306
  errC = clEnqueueWriteBuffer( queue, bufC, CL_TRUE,0,
307
                                 dim*dim*sizeof(*C),C,0,NULL,NULL );
308
  errD = clEnqueueWriteBuffer( queue, bufD, CL_TRUE,0,
309
                                 dim*dim*sizeof(*D),D,0,NULL,NULL );
310
  
311
  /* Get third timer after memory operation */
312
  gettimeofday(&tv3, &tz);
313

    
314
#ifdef FP64
315

    
316
  for (i=0;i<RUNS;i++)
317
    {
318
      err = clblasDgemm( clblasRowMajor,clblasNoTrans,clblasNoTrans, 
319
                         dim,dim,dim,alpha,bufA,0,dim,bufB,0,dim,beta,
320
                         bufC,0,dim,1,&queue,0,NULL,&event );
321

    
322
      err = clblasDgemm( clblasRowMajor,clblasTrans,clblasTrans, 
323
                         dim,dim,dim,alpha,bufB,0,dim,bufA,0,dim,beta,
324
                         bufD,0,dim,1,&queue,0,NULL,&event );
325

    
326
    }
327
  
328
  if (err != CL_SUCCESS) {
329
    printf("clblasDgemm() failed with %d\n", err);
330
  }
331

    
332
#else
333

    
334
  for (i=0;i<RUNS;i++)
335
    {
336

    
337
      err = clblasSgemm( clblasRowMajor,clblasNoTrans,clblasNoTrans, 
338
                         dim,dim,dim,alpha,bufA,0,dim,bufB,0,dim,beta,
339
                         bufC,0,dim,1,&queue,0,NULL,&event );
340

    
341
      err = clblasSgemm( clblasRowMajor,clblasTrans,clblasTrans, 
342
                         dim,dim,dim,alpha,bufB,0,dim,bufA,0,dim,beta,
343
                         bufD,0,dim,1,&queue,0,NULL,&event );
344
    }
345
  
346
  if (err != CL_SUCCESS) {
347
    printf("clblasSgemm() failed with %d\n", err);
348
  }
349

    
350
#endif
351

    
352
  /* Wait for calculations to be finished. */
353
  err = clWaitForEvents( 1, &event );
354
  
355
  /* Get fourth timer after memory free */
356
  gettimeofday(&tv4, &tz);
357

    
358
  /* Fetch results of calculations from GPU memory. */
359
  errC = clEnqueueReadBuffer( queue,bufC,CL_TRUE,0,dim*dim * sizeof(*C),
360
                             C,0,NULL,NULL );
361

    
362
  /* Fetch results of calculations from GPU memory. */
363
  errD = clEnqueueReadBuffer( queue,bufD,CL_TRUE,0,dim*dim*sizeof(*D),
364
                             D,0,NULL,NULL );
365

    
366
  /* Release OpenCL memory objects. */
367
  clReleaseMemObject( bufD );
368
  clReleaseMemObject( bufC );
369
  clReleaseMemObject( bufB );
370
  clReleaseMemObject( bufA );
371

    
372
  /* Finalize work with clBLAS */
373
  clblasTeardown( );
374

    
375
  /* Release OpenCL working objects. */
376
  clReleaseCommandQueue( queue );
377
  clReleaseContext( ctx );
378
  
379

    
380
  /* Compute with CuBLAS library  */
381
#elif CUBLAS
382
  LENGTH *devPtrA=0, *devPtrB=0, *devPtrC=0, *devPtrD=0;
383
  cublasStatus stat1, stat2, stat3, stat4;
384
  struct timeval tv3,tv4;
385

    
386
  /* Order is Row */
387
  /* Have to swap uplo and trans */
388
  char transa='N',transb='T';
389

    
390
  printf("Using CuBLAS: %i iterations for %ix%i matrix\n",
391
         RUNS,dim,dim);
392

    
393
  stat1=cublasAlloc(dim*dim,sizeof(devPtrA[0]),(void**)&devPtrA);
394
  stat2=cublasAlloc(dim*dim,sizeof(devPtrB[0]),(void**)&devPtrB);
395
  stat3=cublasAlloc(dim*dim,sizeof(devPtrC[0]),(void**)&devPtrC);
396
  stat4=cublasAlloc(dim*dim,sizeof(devPtrD[0]),(void**)&devPtrD);
397

    
398
  if ((stat1 != CUBLAS_STATUS_SUCCESS) || 
399
      (stat2 != CUBLAS_STATUS_SUCCESS) ||
400
      (stat3 != CUBLAS_STATUS_SUCCESS) ||
401
      (stat4 != CUBLAS_STATUS_SUCCESS) ) {
402
    wrapperError ("xGEMM", CUBLAS_WRAPPER_ERROR_ALLOC);
403
    cublasFree (devPtrA);
404
    cublasFree (devPtrB);
405
    cublasFree (devPtrC);
406
    cublasFree (devPtrD);
407
    return 1;
408
  }
409

    
410
  stat1=cublasSetMatrix(dim,dim,sizeof(A[0]),A,dim,devPtrA,dim);
411
  stat2=cublasSetMatrix(dim,dim,sizeof(B[0]),B,dim,devPtrB,dim);
412
  stat3=cublasSetMatrix(dim,dim,sizeof(C[0]),C,dim,devPtrC,dim);
413
  stat4=cublasSetMatrix(dim,dim,sizeof(D[0]),D,dim,devPtrD,dim);
414
  
415
  if ((stat1 != CUBLAS_STATUS_SUCCESS) ||
416
      (stat2 != CUBLAS_STATUS_SUCCESS) ||
417
      (stat3 != CUBLAS_STATUS_SUCCESS) ||
418
      (stat4 != CUBLAS_STATUS_SUCCESS) ) {
419
    wrapperError ("xGEMM", CUBLAS_WRAPPER_ERROR_SET);
420
    cublasFree (devPtrA);
421
    cublasFree (devPtrB);
422
    cublasFree (devPtrC);
423
    cublasFree (devPtrD);
424
    return 1;
425
  }
426

    
427
  /* Get third timer after memory operation */
428
  gettimeofday(&tv3, &tz);
429

    
430
#ifdef FP64
431

    
432
  for (i=0;i<RUNS;i++)
433
    {
434
      cublasDgemm(transa,transa,dim,dim,dim,alpha,devPtrB,dim,
435
                  devPtrA,dim,beta,devPtrC,dim);
436
      cublasDgemm(transb,transb,dim,dim,dim,alpha,devPtrA,dim,
437
                  devPtrB,dim,beta,devPtrD,dim);
438
    }
439
  
440
#else
441

    
442
  for (i=0;i<RUNS;i++)
443
    {
444
      cublasSgemm(transa,transa,dim,dim,dim,alpha,devPtrB,dim,
445
                  devPtrA,dim,beta,devPtrC,dim);
446
      cublasSgemm(transb,transb,dim,dim,dim,alpha,devPtrA,dim,
447
                  devPtrB,dim,beta,devPtrD,dim);
448
    }
449
  
450
#endif
451

    
452
  stat3=cublasGetMatrix(dim,dim,sizeof(C[0]),devPtrC,dim,C,dim);
453
  stat4=cublasGetMatrix(dim,dim,sizeof(D[0]),devPtrD,dim,D,dim);
454
  
455
  /* Get fourth timer before memory free */
456
  gettimeofday(&tv4, &tz);
457

    
458
  cublasFree (devPtrA);
459
  cublasFree (devPtrB);
460
  cublasFree (devPtrC);
461
  cublasFree (devPtrD);
462
  
463
  if ((stat1 != CUBLAS_STATUS_SUCCESS) ) {
464
    wrapperError ("xGEMM", CUBLAS_WRAPPER_ERROR_GET);
465
  }
466
  
467

    
468
#elif THUNKING
469
  
470
  /* Order is Row : Have to swap uplo='U' and trans='N' */
471
  char transa='N',transb='T';
472
  printf("Using CuBLAS/Thunking: %i iterations for %ix%i matrix\n",
473
         RUNS,dim,dim);
474

    
475
#ifdef FP64
476

    
477
  for (i=0;i<RUNS;i++)
478
    {      
479
      CUBLAS_DGEMM(&transa,&transa,
480
                         &dim,&dim,&dim,&alpha,B,&dim,A,&dim,&beta,C,&dim);
481
      CUBLAS_DGEMM(&transb,&transb,
482
                         &dim,&dim,&dim,&alpha,A,&dim,B,&dim,&beta,D,&dim);
483
    }
484

    
485
#else
486

    
487
  for (i=0;i<RUNS;i++)
488
    {      
489
      CUBLAS_SGEMM(&transa,&transa,
490
                         &dim,&dim,&dim,&alpha,B,&dim,A,&dim,&beta,C,&dim);
491
      CUBLAS_SGEMM(&transb,&transb,
492
                         &dim,&dim,&dim,&alpha,A,&dim,B,&dim,&beta,D,&dim);
493
    }
494
  
495
#endif
496

    
497
#elif FBLAS
498
  
499
  /* Order is Row : Have to swap uplo='U' and trans='N' */
500
      char transa='N',transb='T';
501
  
502
  printf("Using FBLAS: %i iterations for %ix%i matrix\n",
503
         RUNS,dim,dim);
504
  
505
#ifdef FP64
506

    
507
  for (i=0;i<RUNS;i++)
508
    {    
509
      dgemm_(&transa,&transa,&dim,&dim,&dim,&alpha,B,&dim,A,&dim,&beta,C,&dim);
510
      dgemm_(&transb,&transb,&dim,&dim,&dim,&alpha,A,&dim,B,&dim,&beta,D,&dim);
511
    }
512

    
513
#else
514

    
515
  for (i=0;i<RUNS;i++)
516
    {    
517
      sgemm_(&transa,&transa,&dim,&dim,&dim,&alpha,B,&dim,A,&dim,&beta,C,&dim);
518
      sgemm_(&transb,&transb,&dim,&dim,&dim,&alpha,A,&dim,B,&dim,&beta,D,&dim);
519
    }
520

    
521
#endif
522

    
523
#elif ACML
524
  
525
  /* Order is Row : Have to swap uplo='U' and trans='N' */
526
      char transa='N',transb='T';
527
  
528
  printf("Using ACML: %i iterations for %ix%i matrix\n",
529
         RUNS,dim,dim);
530
  
531
#ifdef FP64
532

    
533
  for (i=0;i<RUNS;i++)
534
    {    
535
      dgemm(transa,transa,dim,dim,dim,alpha,B,dim,A,dim,beta,C,dim);
536
      dgemm(transb,transb,dim,dim,dim,alpha,A,dim,B,dim,beta,D,dim);
537
    }
538

    
539
#else
540

    
541
  for (i=0;i<RUNS;i++)
542
    {    
543
      sgemm(transa,transa,dim,dim,dim,alpha,B,dim,A,dim,beta,C,dim);
544
      sgemm(transb,transb,dim,dim,dim,alpha,A,dim,B,dim,beta,D,dim);
545
    }
546

    
547
#endif
548

    
549
#elif GSL
550

    
551
  printf("Using GSL: %i iterations for %ix%i matrix\n",RUNS,dim,dim);
552

    
553
  /* 
554
     RowMajor : Matrix is read row by row
555
     Upper : the no null elements are on top
556
     NoTrans : no transposition before estimation
557
     NonUnit : Matrix is not unit
558
   */
559

    
560
#ifdef FP64
561

    
562
  for (i=0;i<RUNS;i++)
563
    {  
564
      cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
565
                  dim,dim,dim,alpha,A,dim,B,dim,beta,C,dim);
566
      cblas_dgemm(CblasRowMajor,CblasTrans,CblasTrans,
567
                  dim,dim,dim,alpha,B,dim,A,dim,beta,D,dim);
568
    }
569
  
570
#else
571

    
572
  for (i=0;i<RUNS;i++)
573
    {  
574
      cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
575
                  dim,dim,dim,alpha,A,dim,B,dim,beta,C,dim);
576
      cblas_sgemm(CblasRowMajor,CblasTrans,CblasTrans,
577
                  dim,dim,dim,alpha,B,dim,A,dim,beta,D,dim);
578
    }
579
  
580
#endif
581
      
582
#else
583

    
584
  printf("Using CBLAS: %i iterations for %ix%i matrix\n",RUNS,dim,dim);
585

    
586
  /* 
587
     RowMajor : Matrix is read row bu row
588
     Upper : the no null elements are on top
589
     NoTrans : no transposition before estimation
590
     NonUnit : Matrix is not unit
591
   */
592

    
593
#ifdef FP64
594

    
595
  for (i=0;i<RUNS;i++)
596
    {  
597
      cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
598
                  dim,dim,dim,alpha,A,dim,B,dim,beta,C,dim);
599
      cblas_dgemm(CblasRowMajor,CblasTrans,CblasTrans,
600
                  dim,dim,dim,alpha,B,dim,A,dim,beta,D,dim);
601
    }
602
  
603
#else
604

    
605
  for (i=0;i<RUNS;i++)
606
    {  
607
      cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
608
                  dim,dim,dim,alpha,A,dim,B,dim,beta,C,dim);
609
      cblas_sgemm(CblasRowMajor,CblasTrans,CblasTrans,
610
                  dim,dim,dim,alpha,B,dim,A,dim,beta,D,dim);
611
    }
612
  
613
#endif
614

    
615
#endif
616

    
617
  /* Get second timer after launching */
618
  gettimeofday(&tv2, &tz);
619

    
620
  /* Store the checker of errors */
621
  checksA[0]=0.;
622
  
623
  for (i=0; i<dim; i++) {
624
    for (j=0; j<dim; j++) {
625
      checksA[0]=checksA[0]+fabs(D[i*dim+j]-C[j*dim+i]);
626
    }
627
  }
628

    
629
  /* Print the matrix */
630
 
631
#ifdef QUIET
632
#else
633
  for (i=0; i<dim; i++) {
634
    for (j=0; j<dim; j++) printf("C[%i,%i]=%1.5f ", i,j,C[i*dim+j]);
635
    putchar('\n');
636
  }
637
  putchar('\n');
638
  for (i=0; i<dim; i++) {
639
    for (j=0; j<dim; j++) printf("D[%i,%i]=%1.5f ", i,j,D[i*dim+j]);
640
    putchar('\n');
641
  }
642
  putchar('\n');
643
#endif
644

    
645
  /* Free 1 Matrix and 2 Vectors of dimension dim  */
646

    
647
  free(A);
648
  free(B);
649
  free(C);
650
  free(D);
651

    
652
  putchar('\n');
653

    
654
#ifdef CLBLAS
655
  double memoryIn,memoryOut;
656

    
657
  memoryIn=(double)((tv3.tv_sec-tv1.tv_sec) * 1000000L +        \
658
                    (tv3.tv_usec-tv1.tv_usec))/1000000.;  
659

    
660
  memoryOut=(double)((tv2.tv_sec-tv4.tv_sec) * 1000000L +        \
661
                    (tv2.tv_usec-tv4.tv_usec))/1000000.;  
662

    
663
  duration=(double)((tv4.tv_sec-tv3.tv_sec) * 1000000L +        \
664
                    (tv4.tv_usec-tv3.tv_usec))/1000000./RUNS;  
665

    
666
  printf("Duration of memory allocation : %2.10f s\n",memoryIn);
667
  printf("Duration of memory free : %2.10f s\n",memoryOut);
668
#elif CUBLAS
669
  double memoryIn,memoryOut;
670

    
671
  memoryIn=(double)((tv3.tv_sec-tv1.tv_sec) * 1000000L +        \
672
                    (tv3.tv_usec-tv1.tv_usec))/1000000.;  
673

    
674
  memoryOut=(double)((tv2.tv_sec-tv4.tv_sec) * 1000000L +        \
675
                    (tv2.tv_usec-tv4.tv_usec))/1000000.;  
676

    
677
  duration=(double)((tv4.tv_sec-tv3.tv_sec) * 1000000L +        \
678
                    (tv4.tv_usec-tv3.tv_usec))/1000000./RUNS;  
679

    
680
  printf("Duration of memory allocation : %2.10f s\n",memoryIn);
681
  printf("Duration of memory free : %2.10f s\n",memoryOut);
682
#else
683
  duration=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +        \
684
                    (tv2.tv_usec-tv1.tv_usec))/1000000./RUNS;  
685

    
686
#endif
687

    
688
  printf("Duration of each cycle : %2.10f s\n",duration);
689

    
690
  printf("Number of GFlops : %2.3f \n",
691
         dim*dim*2.*(2.*dim-1)/duration/1000000000.);
692

    
693
  printf("Error %1.10f\n",checksA[0]);
694
  printResults(RUNS,checksA,"C","Errors cumulated");
695

    
696
  putchar('\n');
697

    
698
  /* Free 2 vectors for checker Before and After */
699

    
700
  free(checksA);
701
  free(checksB);
702

    
703
  return 0;
704
}
705

    
706
#ifdef CLBLAS
707

    
708
int DelectOpenCLDevices() 
709
{
710
  /* */
711
  /* Not needed to import CL.h, already done in CLBLAS.h */
712

    
713
  int i, j;
714
  char* value;
715
  size_t valueSize;
716
  cl_uint platformCount;
717
  cl_platform_id* platforms;
718
  cl_uint deviceCount;
719
  cl_device_id* devices;
720
  cl_uint maxComputeUnits;
721
  cl_int maxWorkGroupSize;
722
  cl_int maxWorkItemSizes;
723
  cl_device_type dev_type;
724

    
725
  // get all platforms
726
  clGetPlatformIDs(0, NULL, &platformCount);
727
  platforms = (cl_platform_id*) malloc(sizeof(cl_platform_id) * platformCount);
728
  clGetPlatformIDs(platformCount, platforms, NULL);
729

    
730

    
731
  printf("OpenCL statistics: %d platform(s) detected\n\n",platformCount);
732

    
733
  for (i = 0; i < platformCount; i++) {
734

    
735
    // get all devices
736
    clGetDeviceIDs(platforms[i], CL_DEVICE_TYPE_ALL, 0, NULL, &deviceCount);
737
    devices = (cl_device_id*) malloc(sizeof(cl_device_id) * deviceCount);
738
    clGetDeviceIDs(platforms[i], CL_DEVICE_TYPE_ALL, deviceCount, devices, NULL);
739

    
740
    // for each device print critical attributes
741
    for (j = 0; j < deviceCount; j++) {
742
      
743
      // print device name
744
      clGetDeviceInfo(devices[j], CL_DEVICE_NAME, 0, NULL, &valueSize);
745
      value = (char*) malloc(valueSize);
746
      clGetDeviceInfo(devices[j], CL_DEVICE_NAME, valueSize, value, NULL);
747
      printf("Device (%d,%d): %s\n",i, j, value);
748
      free(value);
749

    
750
      // print type device CPU/GPU/ACCELERATOR
751
      clGetDeviceInfo(devices[j], CL_DEVICE_TYPE, sizeof(dev_type), &dev_type, NULL);
752
      printf("\tDevice Type: ");
753
      if(dev_type & CL_DEVICE_TYPE_GPU)
754
        printf("CL_DEVICE_TYPE_GPU ");
755
      if(dev_type & CL_DEVICE_TYPE_CPU)
756
        printf("CL_DEVICE_TYPE_CPU ");
757
      if(dev_type & CL_DEVICE_TYPE_ACCELERATOR)
758
        printf("CL_DEVICE_TYPE_ACCELERATOR ");
759
      if(dev_type & CL_DEVICE_TYPE_DEFAULT)
760
        printf("CL_DEVICE_TYPE_DEFAULT ");
761
      printf("\n");
762

    
763
      // print device vendor
764
      clGetDeviceInfo(devices[j], CL_DEVICE_VENDOR, 0, NULL, &valueSize);
765
      value = (char*) malloc(valueSize);
766
      clGetDeviceInfo(devices[j], CL_DEVICE_VENDOR, valueSize, value, NULL);
767
      printf("\tDevice vendor: %s\n", value);
768
      free(value);
769

    
770
      // print hardware device version
771
      clGetDeviceInfo(devices[j], CL_DEVICE_VERSION, 0, NULL, &valueSize);
772
      value = (char*) malloc(valueSize);
773
      clGetDeviceInfo(devices[j], CL_DEVICE_VERSION, valueSize, value, NULL);
774
      printf("\tHardware version: %s\n", value);
775
      free(value);
776

    
777
      // print software driver version
778
      clGetDeviceInfo(devices[j], CL_DRIVER_VERSION, 0, NULL, &valueSize);
779
      value = (char*) malloc(valueSize);
780
      clGetDeviceInfo(devices[j], CL_DRIVER_VERSION, valueSize, value, NULL);
781
      printf("\tSoftware version: %s\n", value);
782
      free(value);
783
      
784
      // print c version supported by compiler for device
785
      clGetDeviceInfo(devices[j], CL_DEVICE_OPENCL_C_VERSION, 0, NULL, &valueSize);
786
      value = (char*) malloc(valueSize);
787
      clGetDeviceInfo(devices[j], CL_DEVICE_OPENCL_C_VERSION, valueSize, value, NULL);
788
      printf("\tOpenCL C version: %s\n", value);
789
      free(value);
790

    
791
      // print parallel compute units
792
      clGetDeviceInfo(devices[j], CL_DEVICE_MAX_COMPUTE_UNITS,
793
                      sizeof(maxComputeUnits), &maxComputeUnits, NULL);
794
      printf("\tParallel compute units: %d\n", maxComputeUnits);
795
      
796
      // print max work group size
797
      clGetDeviceInfo(devices[j], CL_DEVICE_MAX_WORK_GROUP_SIZE,
798
                      sizeof(maxWorkGroupSize), &maxWorkGroupSize, NULL);
799
      printf("\tMaximum Work Group Size: %d\n", maxWorkGroupSize);
800
      
801
      // print max work items size
802
      clGetDeviceInfo(devices[j], CL_DEVICE_MAX_WORK_ITEM_SIZES,
803
                      sizeof(maxWorkItemSizes), &maxWorkItemSizes, NULL);
804
      printf("\tMaximum Work Item Sizes: %d\n", maxWorkItemSizes);
805
      
806
    }
807
    printf("\n");
808
    free(devices);
809
  }
810

    
811
  free(platforms);
812
  return 0;
813

    
814
}
815
#endif
816

    
817
int main(int argc,char **argv)
818
{
819
  if ((argc==1)||
820
      (strcmp(argv[1],"-h")==0)||
821
      (strcmp(argv[1],"--help")==0))
822
    {
823
#ifdef CLBLAS
824
      printf("\nPerforms a bench using BLAS library implementation:\n\n"
825
             "\t#1 Size of square matrices \n"
826
             "\t#2 Number of iterations \n"
827
             "\t#3 OpenCL Plateform ID\n"
828
             "\t#4 OpenCL Device ID\n\n");
829
      DelectOpenCLDevices();
830
#else
831
      printf("\nPerforms a bench using BLAS library implementation:\n\n"
832
             "\t#1 Size of square matrices \n"
833
             "\t#2 Number of iterations\n\n");
834
#endif
835
    }
836
  else if ((atoi(argv[1])>=2)&&
837
           (atoi(argv[2])>=1))
838
    {
839
#ifdef CLBLAS
840
      MyPlatform=atoi(argv[3]);
841
      MyDevice=atoi(argv[4]);
842
#endif
843
      bench(atoi(argv[1]),atoi(argv[2]));
844
    }
845

    
846
  return 0;
847
}