Skip CBLAS DSP initializtion if OFFLOAD is disabled.
[dense-linear-algebra-libraries/linalg.git] / examples / dsponly / dgemm_test / dgemm_test.c
1 /******************************************************************************
2  * Copyright (c) 2015, Texas Instruments Incorporated - http://www.ti.com
3  *   All rights reserved.
4  *
5  *   Redistribution and use in source and binary forms, with or without
6  *   modification, are permitted provided that the following conditions are met:
7  *       * Redistributions of source code must retain the above copyright
8  *         notice, this list of conditions and the following disclaimer.
9  *       * Redistributions in binary form must reproduce the above copyright
10  *         notice, this list of conditions and the following disclaimer in the
11  *         documentation and/or other materials provided with the distribution.
12  *       * Neither the name of Texas Instruments Incorporated nor the
13  *         names of its contributors may be used to endorse or promote products
14  *         derived from this software without specific prior written permission.
15  *
16  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17  *   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18  *   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19  *   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
20  *   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21  *   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22  *   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23  *   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24  *   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25  *   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
26  *   THE POSSIBILITY OF SUCH DAMAGE.
27  *****************************************************************************/
28 /******************************************************************************
29 * FILE: dgemm_test.c 
30 * Purpose: to test BLAS function DGEMM.
31 *
32 * DGEMM  performs one of the matrix-matrix operations
33 *
34 *    C := alpha*op( A )*op( B ) + beta*C,
35 *
36 *    where  op( X ) is one of
37 *
38 *    op( X ) = X   or   op( X ) = X**T,
39 *
40 * alpha and beta are scalars, and A, B and C are matrices, with op(A)
41 * an m by k matrix,  op(B) a k by n matrix and C an m by n matrix.
42 ******************************************************************************/
43 #include <string.h>              
44 #include <stdio.h>
45 #include <omp.h>                 /* OpenMP header */
46 #include <ti/libarch/libarch.h>  /* Library Architecture header */
47 #include <ti/linalg/ticblas.h>   /* TI CBLAS API extension header */
48 #include <ti/linalg/cblas.h>     /* Standard CBLAS header */
50 extern void cleanup_after_ticblas();
51 extern void prepare_for_ticblas();
52 extern double omp_get_wtime(void);
54 void matrix_gen(double *A, double *B, double *C, int m, int k, int n);
55 void mat_mpy(const double * A, const double * B, double * C, int mat_N, 
56              int mat_K, int mat_M, double alpha, double beta);
57 double dotprod(const double * A, const double * B, int n);
58 void print_matrix(double *mat, int m, int n);
59 double diff_matrix(double *mat1, double * mat2, int m, int n);
61 int main (int argc, char *argv[]) 
62 {
63   double *A, *B, *C, *C_copy;  /* matrices */
64   int m, n, k;                 /* matrix dimensions */
65   double alpha, beta;          /* scalars for matrix multiplication */
66   double precision_diff;       /* precision difference between TI DGEMM implementation
67                                   and a reference implementation */
68   double time, time_diff;      /* to compute time spent on DGEMM */
69   double gflops;               /* convert time to GFLOPS */
70   int nthreads, tid;           /* number of OpenMp threads and thread id */
72   /* Verify OpenMP working properly */
73   #pragma omp parallel private(nthreads, tid)
74   {
75     tid = omp_get_thread_num(); /* Obtain thread number */
76     printf("Hello World from thread = %d\n", tid);
78     /* Only master thread does this */
79     if (tid == 0) {
80       nthreads = omp_get_num_threads();
81       printf("Number of threads = %d\n", nthreads);
82     }
84   }  /* All threads join master thread and disband */
86   /* hard code dgemm parameters */
87   m = k = n = 1000;
88   alpha = 0.7; 
89   beta  = 1.3; 
91   /* Allocate memory for matrices */
92   A = (double *)malloc( m*k*sizeof( double ) );
93   B = (double *)malloc( k*n*sizeof( double ) );
94   C = (double *)malloc( m*n*sizeof( double ) );
95   C_copy = (double *)malloc( m*n*sizeof( double ) );
96   if (A == NULL || B == NULL || C == NULL || C_copy == NULL) {
97     printf( "\nERROR: Can't allocate memory for matrices. Aborting... \n\n");
98     free(A);
99     free(B);
100     free(C);
101     return 1;
102   }   
104   /* Initialize random number generator */    
105   srand(123456789);
107   /* Configure memory and initialize TI CBLAS */
108   prepare_for_ticblas();
110   /* Generate matrices */
111   matrix_gen(A, B, C, m, k, n);
112   memcpy(C_copy, C, m*n*sizeof(double));  
114   /* Call standard CBLAS API to do matrix multiplication */
115   time = omp_get_wtime();
116   cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, beta, C, n);
117   time_diff = omp_get_wtime() - time;
118   gflops = ( 2.0 * m * n * k ) / time_diff / 1e9;
119   printf("DGEMM time for (m,n,k) = (%d,%d,%d) is %e, GFLOPS is %e.\n", m,n,k, time_diff, gflops);
121   /* Straightforward matrix multiplication as reference */
122   mat_mpy(A, B, C_copy, m, n, k, alpha, beta);
124   /* Find the difference between dgemm and reference */
125   precision_diff = diff_matrix(C, C_copy, m, k);
126   printf("Precision error is %e.\n", precision_diff);
127  
128   /* Finalize TI CBLAS and reconfigure memory */
129   cleanup_after_ticblas();
130   
131   printf("DGEMM test finished and passed.\n");
132   
133   return 0;
136 /*==============================================================================
137  * This function generates matrices of random data
138  *============================================================================*/
139 void matrix_gen(double *A, double *B, double *C, int m, int k, int n)
142     int i;
143     for (i = 0; i < (m*k); i++) {
144         A[i] = (double)rand()/RAND_MAX - 0.5;
145     }
147     for (i = 0; i < (k*n); i++) {
148         B[i] = (double)rand()/RAND_MAX - 0.5;
149     }
151     for (i = 0; i < (m*n); i++) {
152         C[i] = (double)rand()/RAND_MAX - 0.5;
153     }
154     
158 /******************************************************************************
159 * Straightforward implementation of matrix multiplication with row-major
160 ******************************************************************************/
161 void mat_mpy(const double * A, const double * B, double * C, int mat_N, 
162              int mat_K, int mat_M, double alpha, double beta)
164     int col, row;
165     double b_col[mat_K];
167     for (col = 0; col < mat_M; ++col)
168     {
169         for (row = 0; row < mat_K; ++row)
170             b_col[row] = B[row*mat_M+col];
172         for (row = 0; row < mat_N; ++row)
173             C[row*mat_M+col] =  alpha*dotprod(A + (row * mat_K), b_col, mat_K)
174                               + beta*C[row*mat_M+col];
175     }
178 /******************************************************************************
179 * dot product for matrix multiplication
180 ******************************************************************************/
181 double dotprod(const double * A, const double * B, int n)
183     int i;
184     
185     float result = 0;
186     for (i = 0; i < n; ++i) result += A[i] * B[i];
187     
188     return result;
191 /******************************************************************************
192 * Print a row-major matrix
193 ******************************************************************************/
194 void print_matrix(double *mat, int m, int n) 
196    int i, j;
197     
198    for(i=0; i<m; i++) {
199       for(j=0; j<n; j++) {
200          printf( " %10.5f ", mat[i*n+j]);
201       }
202       printf( "\n" );
203   }
206 /******************************************************************************
207 * Find the maximum absolute difference of two matrices
208 ******************************************************************************/
209 double diff_matrix(double *mat1, double * mat2, int m, int n)
211     int i, j;
212     double abs_max_err, err;
213     
214     abs_max_err = 0.0f;
215     for(i=0; i<m; i++) 
216     {
217        for(j=0; j<n; j++) 
218        {
219            err = fabs(mat1[i*n+j] - mat2[i*n+j]);
220            if(abs_max_err < err) {
221                abs_max_err = err;
222            }
223        }
224     }
225     
226     return (abs_max_err);
229 /* Nothing past this point */