MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
testing_zblas.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  * @precisions normal z -> c d s
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  cuDoubleComplex c_neg_one = MAGMA_Z_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  cuDoubleComplex *A, *B, *C, *C2;
39  cuDoubleComplex *dA, *dB, *dC1, *dC2;
40  cuDoubleComplex alpha = MAGMA_Z_MAKE( 0.5, 0.1 );
41  cuDoubleComplex beta = MAGMA_Z_MAKE( 0.7, 0.2 );
42  double dalpha = 0.6;
43  double dbeta = 0.8;
44  double 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_zmalloc_host( &A , size ); assert( err == 0 );
63  err = magma_zmalloc_host( &B , size ); assert( err == 0 );
64  err = magma_zmalloc_host( &C , size ); assert( err == 0 );
65  err = magma_zmalloc_host( &C2, size ); assert( err == 0 );
66  err = magma_zmalloc( &dA, size ); assert( err == 0 );
67  err = magma_zmalloc( &dB, size ); assert( err == 0 );
68  err = magma_zmalloc( &dC1, size ); assert( err == 0 );
69  err = magma_zmalloc( &dC2, size ); assert( err == 0 );
70 
71  // initialize matrices
72  size = maxn*maxn;
73  lapackf77_zlarnv( &ione, ISEED, &size, A );
74  lapackf77_zlarnv( &ione, ISEED, &size, B );
75  lapackf77_zlarnv( &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 zswap\n" );
82  assert( k >= 4 );
83  magma_zsetmatrix( m, k, A, ld, dA, ld );
84  magma_zswap( m, dA(0,1), 1, dA(0,2), 1 );
85  magma_zgetmatrix( m, k, dA, ld, C2, ld );
86  blasf77_zaxpy( &m, &c_neg_one, A(0,0), &ione, C2(0,0), &ione );
87  blasf77_zaxpy( &m, &c_neg_one, A(0,1), &ione, C2(0,2), &ione ); // swapped
88  blasf77_zaxpy( &m, &c_neg_one, A(0,2), &ione, C2(0,1), &ione ); // swapped
89  blasf77_zaxpy( &m, &c_neg_one, A(0,3), &ione, C2(0,3), &ione );
90  size = 4;
91  error = lapackf77_zlange( "F", &m, &size, C2, &ld, work );
92  printf( "zswap diff %.2g\n", error );
93 
94  // ----- test IZAMAX
95  // get argmax of column of A
96  printf( "\ntesting izamax\n" );
97  magma_zsetmatrix( m, k, A, ld, dA, ld );
98  for( int j = 0; j < k; ++j ) {
99  int i1 = magma_izamax( m, dA(0,j), 1 );
100  int i2 = cublasIzamax( 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 zgemv\n" );
111  for( int ia = 0; ia < 3; ++ia ) {
112  magma_zsetmatrix( m, n, A, ld, dA, ld );
113  magma_zsetvector( maxn, B, 1, dB, 1 );
114  magma_zsetvector( maxn, C, 1, dC1, 1 );
115  magma_zsetvector( maxn, C, 1, dC2, 1 );
116  magma_zgemv( trans[ia], m, n, alpha, dA, ld, dB, 1, beta, dC1, 1 );
117  cublasZgemv( 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  cublasZaxpy( size, c_neg_one, dC1, 1, dC2, 1 );
122  magma_zgetvector( size, dC2, 1, C2, 1 );
123  error = lapackf77_zlange( "F", &size, &ione, C2, &ld, work );
124  printf( "zgemv( %c ) diff %.2g\n", trans[ia], error );
125  }
126 
127  // ----- test ZHEMV
128  // c = alpha*A*b + beta*c, with A m*m symmetric; b,c m-vectors
129  // try upper/lower
130  printf( "\ntesting zhemv\n" );
131  for( int iu = 0; iu < 2; ++iu ) {
132  magma_zsetmatrix( m, m, A, ld, dA, ld );
133  magma_zsetvector( m, B, 1, dB, 1 );
134  magma_zsetvector( m, C, 1, dC1, 1 );
135  magma_zsetvector( m, C, 1, dC2, 1 );
136  magma_zhemv( uplo[iu], m, alpha, dA, ld, dB, 1, beta, dC1, 1 );
137  cublasZhemv( 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  cublasZaxpy( m, c_neg_one, dC1, 1, dC2, 1 );
141  magma_zgetvector( m, dC2, 1, C2, 1 );
142  error = lapackf77_zlange( "F", &m, &ione, C2, &ld, work );
143  printf( "zhemv( %c ) diff %.2g\n", uplo[iu], error );
144  }
145 
146  // ----- test ZTRSV
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 ztrsv\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_zgetrf( &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_zsetmatrix( m, m, A, ld, dA, ld );
158  magma_zsetvector( m, C, 1, dC1, 1 );
159  magma_zsetvector( m, C, 1, dC2, 1 );
160  magma_ztrsv( uplo[iu], trans[it], diag[id], m, dA, ld, dC1, 1 );
161  cublasZtrsv( 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  cublasZaxpy( m, c_neg_one, dC1, 1, dC2, 1 );
165  magma_zgetvector( m, dC2, 1, C2, 1 );
166  error = lapackf77_zlange( "F", &m, &ione, C2, &ld, work );
167  printf( "ztrsv( %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 ZGEMM
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 zgemm\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_zsetmatrix( (nta ? m : k), (nta ? m : k), A, ld, dA, ld );
181  magma_zsetmatrix( (ntb ? k : n), (ntb ? n : k), B, ld, dB, ld );
182  magma_zsetmatrix( m, n, C, ld, dC1, ld );
183  magma_zsetmatrix( m, n, C, ld, dC2, ld );
184  magma_zgemm( trans[ia], trans[ib], m, n, k, alpha, dA, ld, dB, ld, beta, dC1, ld );
185  cublasZgemm( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
189  magma_zgetmatrix( m, n, dC2, ld, C2, ld );
190  error = lapackf77_zlange( "F", &m, &n, C2, &ld, work );
191  printf( "zgemm( %c, %c ) diff %.2g\n", trans[ia], trans[ib], error );
192  }}
193 
194  // ----- test ZHEMM
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 zhemm\n" );
199  for( int is = 0; is < 2; ++is ) {
200  for( int iu = 0; iu < 2; ++iu ) {
201  magma_zsetmatrix( m, m, A, ld, dA, ld );
202  magma_zsetmatrix( m, n, B, ld, dB, ld );
203  magma_zsetmatrix( m, n, C, ld, dC1, ld );
204  magma_zsetmatrix( m, n, C, ld, dC2, ld );
205  magma_zhemm( side[is], uplo[iu], m, n, alpha, dA, ld, dB, ld, beta, dC1, ld );
206  cublasZhemm( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
210  magma_zgetmatrix( m, n, dC2, ld, C2, ld );
211  error = lapackf77_zlange( "F", &m, &n, C2, &ld, work );
212  printf( "zhemm( %c, %c ) diff %.2g\n", side[is], uplo[iu], error );
213  }}
214 
215  // ----- test ZHERK
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 zherk\n" );
220  for( int iu = 0; iu < 2; ++iu ) {
221  for( int it = 0; it < 3; ++it ) {
222  magma_zsetmatrix( n, k, A, ld, dA, ld );
223  magma_zsetmatrix( n, n, C, ld, dC1, ld );
224  magma_zsetmatrix( n, n, C, ld, dC2, ld );
225  magma_zherk( uplo[iu], trans[it], n, k, dalpha, dA, ld, dbeta, dC1, ld );
226  cublasZherk( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
230  magma_zgetmatrix( n, n, dC2, ld, C2, ld );
231  error = lapackf77_zlange( "F", &n, &n, C2, &ld, work );
232  printf( "zherk( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
233  }}
234 
235  // ----- test ZHER2K
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 zher2k\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_zsetmatrix( (nt ? n : k), (nt ? n : k), A, ld, dA, ld );
244  magma_zsetmatrix( n, n, C, ld, dC1, ld );
245  magma_zsetmatrix( n, n, C, ld, dC2, ld );
246  magma_zher2k( uplo[iu], trans[it], n, k, alpha, dA, ld, dB, ld, dbeta, dC1, ld );
247  cublasZher2k( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
251  magma_zgetmatrix( n, n, dC2, ld, C2, ld );
252  error = lapackf77_zlange( "F", &n, &n, C2, &ld, work );
253  printf( "zher2k( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
254  }}
255 
256  // ----- test ZTRMM
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 ztrmm\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_zsetmatrix( (left ? m : n), (left ? m : n), A, ld, dA, ld );
267  magma_zsetmatrix( m, n, C, ld, dC1, ld );
268  magma_zsetmatrix( m, n, C, ld, dC2, ld );
269  magma_ztrmm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC1, ld );
270  cublasZtrmm( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
274  magma_zgetmatrix( m, n, dC2, ld, C2, ld );
275  error = lapackf77_zlange( "F", &n, &n, C2, &ld, work );
276  printf( "ztrmm( %c, %c ) diff %.2g\n", uplo[iu], trans[it], error );
277  }}}}
278 
279  // ----- test ZTRSM
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 ztrsm\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_zsetmatrix( (left ? m : n), (left ? m : n), A, ld, dA, ld );
290  magma_zsetmatrix( m, n, C, ld, dC1, ld );
291  magma_zsetmatrix( m, n, C, ld, dC2, ld );
292  magma_ztrsm( side[is], uplo[iu], trans[it], diag[id], m, n, alpha, dA, ld, dC1, ld );
293  cublasZtrsm( 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  cublasZaxpy( ld*n, c_neg_one, dC1, 1, dC2, 1 );
297  magma_zgetmatrix( m, n, dC2, ld, C2, ld );
298  error = lapackf77_zlange( "F", &n, &n, C2, &ld, work );
299  printf( "ztrsm( %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 }