MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
testing_cblas.cpp
Go to the documentation of this file.
1 /*
2  * -- MAGMA (version 1.2.0) --
3  * Univ. of Tennessee, Knoxville
4  * Univ. of California, Berkeley
5  * Univ. of Colorado, Denver
6  * May 2012
7  *
8  * @author Mark Gates
9  * @generated c Thu May 10 22:27:35 2012
10  *
11  **/
12 #include <stdio.h>
13 
14 // make sure that asserts are enabled
15 #ifdef NDEBUG
16 #undef NDEBUG
17 #endif
18 #include <assert.h>
19 
20 #include "common_magma.h"
21 #include "testings.h"
22 
23 #define A(i,j) &A[ (i) + (j)*ld ]
24 #define dA(i,j) &dA[ (i) + (j)*ld ]
25 #define C2(i,j) &C2[ (i) + (j)*ld ]
26 
27 int main( int argc, char** argv )
28 {
29  TESTING_CUDA_INIT();
30 
31  cuFloatComplex c_neg_one = MAGMA_C_NEG_ONE;
32  magma_int_t ione = 1;
33  const char trans[] = { 'N', 'C', 'T' };
34  const char uplo[] = { 'L', 'U' };
35  const char diag[] = { 'U', 'N' };
36  const char side[] = { 'L', 'R' };
37 
38  cuFloatComplex *A, *B, *C, *C2;
39  cuFloatComplex *dA, *dB, *dC1, *dC2;
40  cuFloatComplex alpha = MAGMA_C_MAKE( 0.5, 0.1 );
41  cuFloatComplex beta = MAGMA_C_MAKE( 0.7, 0.2 );
42  float dalpha = 0.6;
43  float dbeta = 0.8;
44  float work[1], error;
45  magma_int_t m = 50;
46  magma_int_t n = 35;
47  magma_int_t k = 40;
48  magma_int_t ISEED[4] = { 0, 1, 2, 3 };
49  magma_int_t size, maxn, ld, info;
50  magma_int_t *piv;
51  magma_err_t err;
52 
53  printf( "\n" );
54 
55  // allocate matrices
56  // over-allocate so they can be any combination of {m,n,k} x {m,n,k}.
57  maxn = max( max( m, n ), k );
58  ld = maxn;
59  size = maxn*maxn;
60  piv = (magma_int_t*) malloc( maxn * sizeof(magma_int_t) );
61  assert( piv != NULL );
62  err = magma_cmalloc_host( &A , size ); assert( err == 0 );
63  err = magma_cmalloc_host( &B , size ); assert( err == 0 );
64  err = magma_cmalloc_host( &C , size ); assert( err == 0 );
65  err = magma_cmalloc_host( &C2, size ); assert( err == 0 );
66  err = magma_cmalloc( &dA, size ); assert( err == 0 );
67  err = magma_cmalloc( &dB, size ); assert( err == 0 );
68  err = magma_cmalloc( &dC1, size ); assert( err == 0 );
69  err = magma_cmalloc( &dC2, size ); assert( err == 0 );
70 
71  // initialize matrices
72  size = maxn*maxn;
73  lapackf77_clarnv( &ione, ISEED, &size, A );
74  lapackf77_clarnv( &ione, ISEED, &size, B );
75  lapackf77_clarnv( &ione, ISEED, &size, C );
76 
77  printf( "========== Level 1 BLAS ==========\n" );
78 
79  // ----- test ZSWAP
80  // swap 2nd and 3rd columns of dA, then copy to C2 and compare with A
81  printf( "\ntesting cswap\n" );
82  assert( k >= 4 );
83  magma_csetmatrix( m, k, A, ld, dA, ld );
84  magma_cswap( m, dA(0,1), 1, dA(0,2), 1 );
85  magma_cgetmatrix( m, k, dA, ld, C2, ld );
86  blasf77_caxpy( &m, &c_neg_one, A(0,0), &ione, C2(0,0), &ione );
87  blasf77_caxpy( &m, &c_neg_one, A(0,1), &ione, C2(0,2), &ione ); // swapped
88  blasf77_caxpy( &m, &c_neg_one, A(0,2), &ione, C2(0,1), &ione ); // swapped
89  blasf77_caxpy( &m, &c_neg_one, A(0,3), &ione, C2(0,3), &ione );
90  size = 4;
91  error = lapackf77_clange( "F", &m, &size, C2, &ld, work );
92  printf( "cswap diff %.2g\n", error );
93 
94  // ----- test IZAMAX
95  // get argmax of column of A
96  printf( "\ntesting icamax\n" );
97  magma_csetmatrix( m, k, A, ld, dA, ld );
98  for( int j = 0; j < k; ++j ) {
99  int i1 = magma_icamax( m, dA(0,j), 1 );
100  int i2 = cublasIcamax( m, dA(0,j), 1 );
101  assert( i1 == i2 );
102  printf( "i1 %4d, i2 %4d, diff %d\n", i1, i2, i1-i2 );
103  }
104 
105  printf( "\n========== Level 2 BLAS ==========\n" );
106 
107  // ----- test ZGEMV
108  // c = alpha*A*b + beta*c, with A m*n; b,c m or n-vectors
109  // try no-trans/trans
110  printf( "\ntesting cgemv\n" );
111  for( int ia = 0; ia < 3; ++ia ) {
112  magma_csetmatrix( m, n, A, ld, dA, ld );
113  magma_csetvector( maxn, B, 1, dB, 1 );
114  magma_csetvector( maxn, C, 1, dC1, 1 );
115  magma_csetvector( maxn, C, 1, dC2, 1 );
116  magma_cgemv( trans[ia], m, n, alpha, dA, ld, dB, 1, beta, dC1, 1 );
117  cublasCgemv( trans[ia], m, n, alpha, dA, ld, dB, 1, beta, dC2, 1 );
118 
119  // check results, storing diff between magma and cuda call in C2
120  size = (trans[ia] == 'N' ? m : n);
121  cublasCaxpy( size, c_neg_one, dC1, 1, dC2, 1 );
122  magma_cgetvector( size, dC2, 1, C2, 1 );
123  error = lapackf77_clange( "F", &size, &ione, C2, &ld, work );
124  printf( "cgemv( %c ) diff %.2g\n", trans[ia], error );
125  }
126 
127  // ----- test CHEMV
128  // c = alpha*A*b + beta*c, with A m*m symmetric; b,c m-vectors
129  // try upper/lower
130  printf( "\ntesting chemv\n" );
131  for( int iu = 0; iu < 2; ++iu ) {
132  magma_csetmatrix( m, m, A, ld, dA, ld );
133  magma_csetvector( m, B, 1, dB, 1 );
134  magma_csetvector( m, C, 1, dC1, 1 );
135  magma_csetvector( m, C, 1, dC2, 1 );
136  magma_chemv( uplo[iu], m, alpha, dA, ld, dB, 1, beta, dC1, 1 );
137  cublasChemv( uplo[iu], m, alpha, dA, ld, dB, 1, beta, dC2, 1 );
138 
139  // check results, storing diff between magma and cuda call in C2
140  cublasCaxpy( m, c_neg_one, dC1, 1, dC2, 1 );
141  magma_cgetvector( m, dC2, 1, C2, 1 );
142  error = lapackf77_clange( "F", &m, &ione, C2, &ld, work );
143  printf( "chemv( %c ) diff %.2g\n", uplo[iu], error );
144  }
145 
146  // ----- test CTRSV
147  // solve A*c = c, with A m*m triangular; c m-vector
148  // try upper/lower, no-trans/trans, unit/non-unit diag
149  printf( "\ntesting ctrsv\n" );
150  // Factor A into LU to get well-conditioned triangles, else solve yields garbage.
151  // Still can give garbage if solves aren't consistent with LU factors,
152  // e.g., using unit diag for U.
153  lapackf77_cgetrf( &m, &m, A, &ld, piv, &info );
154  for( int iu = 0; iu < 2; ++iu ) {
155  for( int it = 0; it < 3; ++it ) {
156  for( int id = 0; id < 2; ++id ) {
157  magma_csetmatrix( m, m, A, ld, dA, ld );
158  magma_csetvector( m, C, 1, dC1, 1 );
159  magma_csetvector( m, C, 1, dC2, 1 );
160  magma_ctrsv( uplo[iu], trans[it], diag[id], m, dA, ld, dC1, 1 );
161  cublasCtrsv( uplo[iu], trans[it], diag[id], m, dA, ld, dC2, 1 );
162 
163  // check results, storing diff between magma and cuda call in C2
164  cublasCaxpy( m, c_neg_one, dC1, 1, dC2, 1 );
165  magma_cgetvector( m, dC2, 1, C2, 1 );
166  error = lapackf77_clange( "F", &m, &ione, C2, &ld, work );
167  printf( "ctrsv( %c, %c, %c ) diff %.2g\n", uplo[iu], trans[it], diag[id], error );
168  }}}
169 
170  printf( "\n========== Level 3 BLAS ==========\n" );
171 
172  // ----- test CGEMM
173  // C = alpha*A*B + beta*C, with A m*k or k*m; B k*n or n*k; C m*n
174  // try combinations of no-trans/trans
175  printf( "\ntesting cgemm\n" );
176  for( int ia = 0; ia < 3; ++ia ) {
177  for( int ib = 0; ib < 3; ++ib ) {
178  bool nta = (trans[ia] == 'N');
179  bool ntb = (trans[ib] == 'N');
180  magma_csetmatrix( (nta ? m : k), (nta ? m : k), A, ld, dA, ld );
181  magma_csetmatrix( (ntb ? k : n), (ntb ? n : k), B, ld, dB, ld );
182  magma_csetmatrix( m, n, C, ld, dC1, ld );
183  magma_csetmatrix( m, n, C, ld, dC2, ld );
184  magma_cgemm( trans[ia], trans[ib], m, n, k, alpha, dA, ld, dB, ld, beta, dC1, ld );
185  cublasCgemm( trans[ia], trans[ib], m, n, k, alpha, dA, ld, dB, ld, beta, dC2, ld );
186 
187  // check results, storing diff between magma and cuda call in C2
188  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
189  magma_cgetmatrix( m, n, dC2, ld, C2, ld );
190  error = lapackf77_clange( "F", &m, &n, C2, &ld, work );
191  printf( "cgemm( %c, %c ) diff %.2g\n", trans[ia], trans[ib], error );
192  }}
193 
194  // ----- test CHEMM
195  // C = alpha*A*B + beta*C (left) with A m*m symmetric; B,C m*n; or
196  // C = alpha*B*A + beta*C (right) with A n*n symmetric; B,C m*n
197  // try left/right, upper/lower
198  printf( "\ntesting chemm\n" );
199  for( int is = 0; is < 2; ++is ) {
200  for( int iu = 0; iu < 2; ++iu ) {
201  magma_csetmatrix( m, m, A, ld, dA, ld );
202  magma_csetmatrix( m, n, B, ld, dB, ld );
203  magma_csetmatrix( m, n, C, ld, dC1, ld );
204  magma_csetmatrix( m, n, C, ld, dC2, ld );
205  magma_chemm( side[is], uplo[iu], m, n, alpha, dA, ld, dB, ld, beta, dC1, ld );
206  cublasChemm( side[is], uplo[iu], m, n, alpha, dA, ld, dB, ld, beta, dC2, ld );
207 
208  // check results, storing diff between magma and cuda call in C2
209  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
210  magma_cgetmatrix( m, n, dC2, ld, C2, ld );
211  error = lapackf77_clange( "F", &m, &n, C2, &ld, work );
212  printf( "chemm( %c, %c ) diff %.2g\n", side[is], uplo[iu], error );
213  }}
214 
215  // ----- test CHERK
216  // C = alpha*A*A^H + beta*C (no-trans) with A m*k and C m*m symmetric; or
217  // C = alpha*A^H*A + beta*C (trans) with A k*m and C m*m symmetric
218  // try upper/lower, no-trans/trans
219  printf( "\ntesting cherk\n" );
220  for( int iu = 0; iu < 2; ++iu ) {
221  for( int it = 0; it < 3; ++it ) {
222  magma_csetmatrix( n, k, A, ld, dA, ld );
223  magma_csetmatrix( n, n, C, ld, dC1, ld );
224  magma_csetmatrix( n, n, C, ld, dC2, ld );
225  magma_cherk( uplo[iu], trans[it], n, k, dalpha, dA, ld, dbeta, dC1, ld );
226  cublasCherk( uplo[iu], trans[it], n, k, dalpha, dA, ld, dbeta, dC2, ld );
227 
228  // check results, storing diff between magma and cuda call in C2
229  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
230  magma_cgetmatrix( n, n, dC2, ld, C2, ld );
231  error = lapackf77_clange( "F", &n, &n, C2, &ld, work );
232  printf( "cherk( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
233  }}
234 
235  // ----- test CHER2K
236  // C = alpha*A*B^H + ^alpha*B*A^H + beta*C (no-trans) with A,B n*k; C n*n symmetric; or
237  // C = alpha*A^H*B + ^alpha*B^H*A + beta*C (trans) with A,B k*n; C n*n symmetric
238  // try upper/lower, no-trans/trans
239  printf( "\ntesting cher2k\n" );
240  for( int iu = 0; iu < 2; ++iu ) {
241  for( int it = 0; it < 3; ++it ) {
242  bool nt = (trans[it] == 'N');
243  magma_csetmatrix( (nt ? n : k), (nt ? n : k), A, ld, dA, ld );
244  magma_csetmatrix( n, n, C, ld, dC1, ld );
245  magma_csetmatrix( n, n, C, ld, dC2, ld );
246  magma_cher2k( uplo[iu], trans[it], n, k, alpha, dA, ld, dB, ld, dbeta, dC1, ld );
247  cublasCher2k( uplo[iu], trans[it], n, k, alpha, dA, ld, dB, ld, dbeta, dC2, ld );
248 
249  // check results, storing diff between magma and cuda call in C2
250  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
251  magma_cgetmatrix( n, n, dC2, ld, C2, ld );
252  error = lapackf77_clange( "F", &n, &n, C2, &ld, work );
253  printf( "cher2k( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
254  }}
255 
256  // ----- test CTRMM
257  // C = alpha*A*C (left) with A m*m triangular; C m*n; or
258  // C = alpha*C*A (right) with A n*n triangular; C m*n
259  // try left/right, upper/lower, no-trans/trans, unit/non-unit
260  printf( "\ntesting ctrmm\n" );
261  for( int is = 0; is < 2; ++is ) {
262  for( int iu = 0; iu < 2; ++iu ) {
263  for( int it = 0; it < 3; ++it ) {
264  for( int id = 0; id < 2; ++id ) {
265  bool left = (side[is] == 'L');
266  magma_csetmatrix( (left ? m : n), (left ? m : n), A, ld, dA, ld );
267  magma_csetmatrix( m, n, C, ld, dC1, ld );
268  magma_csetmatrix( m, n, C, ld, dC2, ld );
269  magma_ctrmm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC1, ld );
270  cublasCtrmm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC2, ld );
271 
272  // check results, storing diff between magma and cuda call in C2
273  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
274  magma_cgetmatrix( m, n, dC2, ld, C2, ld );
275  error = lapackf77_clange( "F", &n, &n, C2, &ld, work );
276  printf( "ctrmm( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
277  }}}}
278 
279  // ----- test CTRSM
280  // solve A*X = alpha*B (left) with A m*m triangular; B m*n; or
281  // solve X*A = alpha*B (right) with A n*n triangular; B m*n
282  // try left/right, upper/lower, no-trans/trans, unit/non-unit
283  printf( "\ntesting ctrsm\n" );
284  for( int is = 0; is < 2; ++is ) {
285  for( int iu = 0; iu < 2; ++iu ) {
286  for( int it = 0; it < 3; ++it ) {
287  for( int id = 0; id < 2; ++id ) {
288  bool left = (side[is] == 'L');
289  magma_csetmatrix( (left ? m : n), (left ? m : n), A, ld, dA, ld );
290  magma_csetmatrix( m, n, C, ld, dC1, ld );
291  magma_csetmatrix( m, n, C, ld, dC2, ld );
292  magma_ctrsm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC1, ld );
293  cublasCtrsm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC2, ld );
294 
295  // check results, storing diff between magma and cuda call in C2
296  cublasCaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
297  magma_cgetmatrix( m, n, dC2, ld, C2, ld );
298  error = lapackf77_clange( "F", &n, &n, C2, &ld, work );
299  printf( "ctrsm( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
300  }}}}
301 
302  // cleanup
303  magma_free_host( A );
304  magma_free_host( B );
305  magma_free_host( C );
306  magma_free_host( C2 );
307  magma_free( dA );
308  magma_free( dB );
309  magma_free( dC1 );
310  magma_free( dC2 );
311 
312  TESTING_CUDA_FINALIZE();
313  return 0;
314 }