1 /* testMMM.c */
2
3 #include "../spoolesMT.h"
4 #include "../../Drand.h"
5 #include "../../timings.h"
6
7 /*--------------------------------------------------------------------*/
8 int
main(int argc,char * argv[])9 main ( int argc, char *argv[] )
10 /*
11 ------------------------------------------------------------------
12 generate a random matrix and test a matrix-matrix multiply method.
13 the output is a matlab file to test correctness.
14
15 created -- 98jan29, cca
16 --------------------------------------------------------------------
17 */
18 {
19 DenseMtx *X, *Y, *Y2 ;
20 double alpha[2] ;
21 double alphaImag, alphaReal, t1, t2 ;
22 double *zvec ;
23 Drand *drand ;
24 int col, dataType, ii, msglvl, ncolA, nitem, nops, nrhs,
25 nrowA, nrowX, nrowY, nthread, row, seed,
26 storageMode, symflag, transposeflag ;
27 int *colids, *rowids ;
28 InpMtx *A ;
29 FILE *msgFile ;
30
31 if ( argc != 15 ) {
32 fprintf(stdout,
33 "\n\n %% usage : %s msglvl msgFile symflag storageMode "
34 "\n %% nrow ncol nent nrhs seed alphaReal alphaImag nthread"
35 "\n %% msglvl -- message level"
36 "\n %% msgFile -- message file"
37 "\n %% dataType -- type of matrix entries"
38 "\n %% 1 -- real"
39 "\n %% 2 -- complex"
40 "\n %% symflag -- symmetry flag"
41 "\n %% 0 -- symmetric"
42 "\n %% 1 -- hermitian"
43 "\n %% 2 -- nonsymmetric"
44 "\n %% storageMode -- storage mode"
45 "\n %% 1 -- by rows"
46 "\n %% 2 -- by columns"
47 "\n %% 3 -- by chevrons, (requires nrow = ncol)"
48 "\n %% transpose -- transpose flag"
49 "\n %% 0 -- Y := Y + alpha * A * X"
50 "\n %% 1 -- Y := Y + alpha * A^H * X, nonsymmetric only"
51 "\n %% 2 -- Y := Y + alpha * A^T * X, nonsymmetric only"
52 "\n %% nrowA -- number of rows in A"
53 "\n %% ncolA -- number of columns in A"
54 "\n %% nitem -- number of items"
55 "\n %% nrhs -- number of right hand sides"
56 "\n %% seed -- random number seed"
57 "\n %% alphaReal -- y := y + alpha*A*x"
58 "\n %% alphaImag -- y := y + alpha*A*x"
59 "\n %% nthread -- # of threads"
60 "\n", argv[0]) ;
61 return(0) ;
62 }
63 msglvl = atoi(argv[1]) ;
64 if ( strcmp(argv[2], "stdout") == 0 ) {
65 msgFile = stdout ;
66 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
67 fprintf(stderr, "\n fatal error in %s"
68 "\n unable to open file %s\n",
69 argv[0], argv[2]) ;
70 return(-1) ;
71 }
72 dataType = atoi(argv[3]) ;
73 symflag = atoi(argv[4]) ;
74 storageMode = atoi(argv[5]) ;
75 transposeflag = atoi(argv[6]) ;
76 nrowA = atoi(argv[7]) ;
77 ncolA = atoi(argv[8]) ;
78 nitem = atoi(argv[9]) ;
79 nrhs = atoi(argv[10]) ;
80 seed = atoi(argv[11]) ;
81 alphaReal = atof(argv[12]) ;
82 alphaImag = atof(argv[13]) ;
83 nthread = atoi(argv[14]) ;
84 fprintf(msgFile,
85 "\n %% %s "
86 "\n %% msglvl -- %d"
87 "\n %% msgFile -- %s"
88 "\n %% dataType -- %d"
89 "\n %% symflag -- %d"
90 "\n %% storageMode -- %d"
91 "\n %% transposeflag -- %d"
92 "\n %% nrowA -- %d"
93 "\n %% ncolA -- %d"
94 "\n %% nitem -- %d"
95 "\n %% nrhs -- %d"
96 "\n %% seed -- %d"
97 "\n %% alphaReal -- %e"
98 "\n %% alphaImag -- %e"
99 "\n %% nthread -- %d"
100 "\n",
101 argv[0], msglvl, argv[2], dataType, symflag, storageMode,
102 transposeflag, nrowA, ncolA, nitem, nrhs, seed,
103 alphaReal, alphaImag, nthread) ;
104 fflush(msgFile) ;
105 if ( dataType != 1 && dataType != 2 ) {
106 fprintf(stderr, "\n invalid value %d for dataType\n", dataType) ;
107 exit(-1) ;
108 }
109 if ( symflag != 0 && symflag != 1 && symflag != 2 ) {
110 fprintf(stderr, "\n invalid value %d for symflag\n", symflag) ;
111 exit(-1) ;
112 }
113 if ( storageMode != 1 && storageMode != 2 && storageMode != 3 ) {
114 fprintf(stderr,
115 "\n invalid value %d for storageMode\n", storageMode) ;
116 exit(-1) ;
117 }
118 if ( transposeflag < 0
119 || transposeflag > 2 ) {
120 fprintf(stderr, "\n error, transposeflag = %d, must be 0, 1 or 2",
121 transposeflag) ;
122 exit(-1) ;
123 }
124 if ( (transposeflag == 1 && symflag != 2)
125 || (transposeflag == 2 && symflag != 2) ) {
126 fprintf(stderr, "\n error, transposeflag = %d, symflag = %d",
127 transposeflag, symflag) ;
128 exit(-1) ;
129 }
130 if ( transposeflag == 1 && dataType != 2 ) {
131 fprintf(stderr, "\n error, transposeflag = %d, dataType = %d",
132 transposeflag, dataType) ;
133 exit(-1) ;
134 }
135 if ( symflag == 1 && dataType != 2 ) {
136 fprintf(stderr,
137 "\n symflag = 1 (hermitian), dataType != 2 (complex)") ;
138 exit(-1) ;
139 }
140 if ( nrowA <= 0 || ncolA <= 0 || nitem <= 0 ) {
141 fprintf(stderr,
142 "\n invalid value: nrow = %d, ncol = %d, nitem = %d",
143 nrowA, ncolA, nitem) ;
144 exit(-1) ;
145 }
146 if ( symflag < 2 && nrowA != ncolA ) {
147 fprintf(stderr,
148 "\n invalid data: symflag = %d, nrow = %d, ncol = %d",
149 symflag, nrowA, ncolA) ;
150 exit(-1) ;
151 }
152 alpha[0] = alphaReal ;
153 alpha[1] = alphaImag ;
154 /*
155 ----------------------------
156 initialize the matrix object
157 ----------------------------
158 */
159 A = InpMtx_new() ;
160 InpMtx_init(A, storageMode, dataType, 0, 0) ;
161 drand = Drand_new() ;
162 /*
163 ----------------------------------
164 generate a vector of nitem triples
165 ----------------------------------
166 */
167 rowids = IVinit(nitem, -1) ;
168 Drand_setUniform(drand, 0, nrowA) ;
169 Drand_fillIvector(drand, nitem, rowids) ;
170 colids = IVinit(nitem, -1) ;
171 Drand_setUniform(drand, 0, ncolA) ;
172 Drand_fillIvector(drand, nitem, colids) ;
173 Drand_setUniform(drand, 0.0, 1.0) ;
174 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
175 zvec = DVinit(nitem, 0.0) ;
176 Drand_fillDvector(drand, nitem, zvec) ;
177 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
178 zvec = ZVinit(nitem, 0.0, 0.0) ;
179 Drand_fillDvector(drand, 2*nitem, zvec) ;
180 }
181 /*
182 -----------------------------------
183 assemble the entries entry by entry
184 -----------------------------------
185 */
186 if ( msglvl > 1 ) {
187 fprintf(msgFile, "\n\n A = zeros(%d,%d) ;", nrowA, ncolA) ;
188 }
189 if ( symflag == 1 ) {
190 /*
191 ----------------
192 hermitian matrix
193 ----------------
194 */
195 for ( ii = 0 ; ii < nitem ; ii++ ) {
196 if ( rowids[ii] == colids[ii] ) {
197 zvec[2*ii+1] = 0.0 ;
198 }
199 if ( rowids[ii] <= colids[ii] ) {
200 row = rowids[ii] ; col = colids[ii] ;
201 } else {
202 row = colids[ii] ; col = rowids[ii] ;
203 }
204 InpMtx_inputComplexEntry(A, row, col, zvec[2*ii], zvec[2*ii+1]) ;
205 }
206 } else if ( symflag == 0 ) {
207 /*
208 ----------------
209 symmetric matrix
210 ----------------
211 */
212 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
213 for ( ii = 0 ; ii < nitem ; ii++ ) {
214 if ( rowids[ii] <= colids[ii] ) {
215 row = rowids[ii] ; col = colids[ii] ;
216 } else {
217 row = colids[ii] ; col = rowids[ii] ;
218 }
219 InpMtx_inputRealEntry(A, row, col, zvec[ii]) ;
220 }
221 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
222 for ( ii = 0 ; ii < nitem ; ii++ ) {
223 if ( rowids[ii] <= colids[ii] ) {
224 row = rowids[ii] ; col = colids[ii] ;
225 } else {
226 row = colids[ii] ; col = rowids[ii] ;
227 }
228 InpMtx_inputComplexEntry(A, row, col,
229 zvec[2*ii], zvec[2*ii+1]) ;
230 }
231 }
232 } else {
233 /*
234 -------------------
235 nonsymmetric matrix
236 -------------------
237 */
238 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
239 for ( ii = 0 ; ii < nitem ; ii++ ) {
240 InpMtx_inputRealEntry(A, rowids[ii], colids[ii], zvec[ii]) ;
241 }
242 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
243 for ( ii = 0 ; ii < nitem ; ii++ ) {
244 InpMtx_inputComplexEntry(A, rowids[ii], colids[ii],
245 zvec[2*ii], zvec[2*ii+1]) ;
246 }
247 }
248 }
249 InpMtx_changeStorageMode(A, INPMTX_BY_VECTORS) ;
250 DVfree(zvec) ;
251 if ( symflag == 0 || symflag == 1 ) {
252 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
253 nops = 4*A->nent*nrhs ;
254 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
255 nops = 16*A->nent*nrhs ;
256 }
257 } else {
258 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
259 nops = 2*A->nent*nrhs ;
260 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
261 nops = 8*A->nent*nrhs ;
262 }
263 }
264 if ( msglvl > 1 ) {
265 /*
266 -------------------------------------------
267 write the assembled matrix to a matlab file
268 -------------------------------------------
269 */
270 InpMtx_writeForMatlab(A, "A", msgFile) ;
271 if ( symflag == 0 ) {
272 fprintf(msgFile,
273 "\n for k = 1:%d"
274 "\n for j = k+1:%d"
275 "\n A(j,k) = A(k,j) ;"
276 "\n end"
277 "\n end", nrowA, ncolA) ;
278 } else if ( symflag == 1 ) {
279 fprintf(msgFile,
280 "\n for k = 1:%d"
281 "\n for j = k+1:%d"
282 "\n A(j,k) = ctranspose(A(k,j)) ;"
283 "\n end"
284 "\n end", nrowA, ncolA) ;
285 }
286 }
287 /*
288 -------------------------------
289 generate dense matrices X and Y
290 -------------------------------
291 */
292 if ( transposeflag == 0 ) {
293 nrowX = ncolA ;
294 nrowY = nrowA ;
295 } else {
296 nrowX = nrowA ;
297 nrowY = ncolA ;
298 }
299 X = DenseMtx_new() ;
300 Y = DenseMtx_new() ;
301 Y2 = DenseMtx_new() ;
302 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
303 DenseMtx_init(X, SPOOLES_REAL, 0, 0, nrowX, nrhs, 1, nrowX) ;
304 Drand_fillDvector(drand, nrowX*nrhs, DenseMtx_entries(X)) ;
305 DenseMtx_init(Y, SPOOLES_REAL, 0, 0, nrowY, nrhs, 1, nrowY) ;
306 Drand_fillDvector(drand, nrowY*nrhs, DenseMtx_entries(Y)) ;
307 DenseMtx_init(Y2, SPOOLES_REAL, 0, 0, nrowY, nrhs, 1, nrowY) ;
308 DVcopy(nrowY*nrhs, DenseMtx_entries(Y2), DenseMtx_entries(Y)) ;
309 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
310 DenseMtx_init(X, SPOOLES_COMPLEX, 0, 0, nrowX, nrhs, 1, nrowX) ;
311 Drand_fillDvector(drand, 2*nrowX*nrhs, DenseMtx_entries(X)) ;
312 DenseMtx_init(Y, SPOOLES_COMPLEX, 0, 0, nrowY, nrhs, 1, nrowY) ;
313 Drand_fillDvector(drand, 2*nrowY*nrhs, DenseMtx_entries(Y)) ;
314 DenseMtx_init(Y2, SPOOLES_COMPLEX, 0, 0, nrowY, nrhs, 1, nrowY) ;
315 DVcopy(2*nrowY*nrhs, DenseMtx_entries(Y2), DenseMtx_entries(Y)) ;
316 }
317 if ( msglvl > 1 ) {
318 fprintf(msgFile, "\n X = zeros(%d,%d) ;", nrowX, nrhs) ;
319 DenseMtx_writeForMatlab(X, "X", msgFile) ;
320 fprintf(msgFile, "\n Y = zeros(%d,%d) ;", nrowY, nrhs) ;
321 DenseMtx_writeForMatlab(Y, "Y", msgFile) ;
322 }
323 /*
324 --------------------------------------------
325 perform the matrix-matrix multiply in serial
326 --------------------------------------------
327 */
328 if ( msglvl > 1 ) {
329 fprintf(msgFile, "\n alpha = %20.12e + %20.2e*i;",
330 alpha[0], alpha[1]);
331 fprintf(msgFile, "\n Z = zeros(%d,1) ;", nrowY) ;
332 }
333 if ( transposeflag == 0 ) {
334 MARKTIME(t1) ;
335 if ( symflag == 0 ) {
336 InpMtx_sym_mmm(A, Y, alpha, X) ;
337 } else if ( symflag == 1 ) {
338 InpMtx_herm_mmm(A, Y, alpha, X) ;
339 } else if ( symflag == 2 ) {
340 InpMtx_nonsym_mmm(A, Y, alpha, X) ;
341 }
342 MARKTIME(t2) ;
343 if ( msglvl > 1 ) {
344 DenseMtx_writeForMatlab(Y, "Z", msgFile) ;
345 fprintf(msgFile, "\n maxerr = max(Z - Y - alpha*A*X) ") ;
346 fprintf(msgFile, "\n") ;
347 }
348 } else if ( transposeflag == 1 ) {
349 MARKTIME(t1) ;
350 InpMtx_nonsym_mmm_H(A, Y, alpha, X) ;
351 MARKTIME(t2) ;
352 if ( msglvl > 1 ) {
353 DenseMtx_writeForMatlab(Y, "Z", msgFile) ;
354 fprintf(msgFile,
355 "\n maxerr = max(Z - Y - alpha*ctranspose(A)*X) ") ;
356 fprintf(msgFile, "\n") ;
357 }
358 } else if ( transposeflag == 2 ) {
359 MARKTIME(t1) ;
360 InpMtx_nonsym_mmm_T(A, Y, alpha, X) ;
361 MARKTIME(t2) ;
362 if ( msglvl > 1 ) {
363 DenseMtx_writeForMatlab(Y, "Z", msgFile) ;
364 fprintf(msgFile,
365 "\n maxerr = max(Z - Y - alpha*transpose(A)*X) ") ;
366 fprintf(msgFile, "\n") ;
367 }
368 }
369 fprintf(msgFile, "\n %% %d ops, %.3f time, %.3f serial mflops",
370 nops, t2 - t1, 1.e-6*nops/(t2 - t1)) ;
371 /*
372 --------------------------------------------------------
373 perform the matrix-matrix multiply in multithreaded mode
374 --------------------------------------------------------
375 */
376 if ( msglvl > 1 ) {
377 fprintf(msgFile,
378 "\n alpha = %20.12e + %20.2e*i;", alpha[0], alpha[1]);
379 fprintf(msgFile, "\n Z = zeros(%d,1) ;", nrowY) ;
380 }
381 if ( transposeflag == 0 ) {
382 MARKTIME(t1) ;
383 if ( symflag == 0 ) {
384 InpMtx_MT_sym_mmm(A, Y2, alpha, X, nthread, msglvl, msgFile) ;
385 } else if ( symflag == 1 ) {
386 InpMtx_MT_herm_mmm(A, Y2, alpha, X, nthread, msglvl, msgFile) ;
387 } else if ( symflag == 2 ) {
388 InpMtx_MT_nonsym_mmm(A, Y2, alpha, X, nthread, msglvl, msgFile) ;
389 }
390 MARKTIME(t2) ;
391 if ( msglvl > 1 ) {
392 DenseMtx_writeForMatlab(Y2, "Z2", msgFile) ;
393 fprintf(msgFile, "\n maxerr2 = max(Z2 - Y - alpha*A*X) ") ;
394 fprintf(msgFile, "\n") ;
395 }
396 } else if ( transposeflag == 1 ) {
397 MARKTIME(t1) ;
398 InpMtx_MT_nonsym_mmm_H(A, Y2, alpha, X, nthread, msglvl, msgFile) ;
399 MARKTIME(t2) ;
400 if ( msglvl > 1 ) {
401 DenseMtx_writeForMatlab(Y2, "Z2", msgFile) ;
402 fprintf(msgFile,
403 "\n maxerr2 = max(Z2 - Y - alpha*ctranspose(A)*X) ") ;
404 fprintf(msgFile, "\n") ;
405 }
406 } else if ( transposeflag == 2 ) {
407 MARKTIME(t1) ;
408 InpMtx_MT_nonsym_mmm_T(A, Y2, alpha, X, nthread, msglvl, msgFile) ;
409 MARKTIME(t2) ;
410 if ( msglvl > 1 ) {
411 DenseMtx_writeForMatlab(Y2, "Z2", msgFile) ;
412 fprintf(msgFile,
413 "\n maxerr2 = max(Z2 - Y - alpha*transpose(A)*X) ") ;
414 fprintf(msgFile, "\n") ;
415 }
416 }
417 fprintf(msgFile, "\n %% %d ops, %.3f time, %.3f MT mflops",
418 nops, t2 - t1, 1.e-6*nops/(t2 - t1)) ;
419 /*
420 ------------------------
421 free the working storage
422 ------------------------
423 */
424 InpMtx_free(A) ;
425 DenseMtx_free(X) ;
426 DenseMtx_free(Y) ;
427 DenseMtx_free(Y2) ;
428 IVfree(rowids) ;
429 IVfree(colids) ;
430 Drand_free(drand) ;
431
432 fclose(msgFile) ;
433
434 return(1) ; }
435
436 /*--------------------------------------------------------------------*/
437