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