1 /* test_solveupd.c */
2
3 #include "../SubMtx.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 test the SubMtx_solveupd() method.
13
14 created -- 98may02, cca
15 ----------------------------------
16 */
17 {
18 SubMtx *mtxA, *mtxX, *mtxY ;
19 double ops, t1, t2 ;
20 double *entX, *entY ;
21 Drand *drand ;
22 FILE *msgFile ;
23 int inc1, inc2, irowA, jcolX, mode, msglvl, ncolA, nentA, nrowA,
24 ncolX, nrowX, ncolY, nrowY, seed, type ;
25 int *colindA, *ivec, *rowindA ;
26
27 if ( argc != 12 ) {
28 fprintf(stdout,
29 "\n\n usage : %s msglvl msgFile type mode"
30 "\n nrowY ncolY nrowA ncolA nentA nrowX seed"
31 "\n msglvl -- message level"
32 "\n msgFile -- message file"
33 "\n type -- type of matrix A"
34 "\n 1 -- real"
35 "\n 2 -- complex"
36 "\n mode -- mode of matrix A"
37 "\n 0 -- dense stored by rows"
38 "\n 1 -- dense stored by columns"
39 "\n 2 -- sparse stored by rows"
40 "\n 3 -- sparse stored by columns"
41 "\n nrowY -- # of rows in vector Y"
42 "\n ncolY -- # of columns in vector Y"
43 "\n nrowA -- # of rows in matrix A"
44 "\n ncolA -- # of columns in matrix A"
45 "\n nentA -- # of entries in matrix A"
46 "\n nrowX -- # of rows in matrix X"
47 "\n seed -- random number seed"
48 "\n", argv[0]) ;
49 return(0) ;
50 }
51 if ( (msglvl = atoi(argv[1])) < 0 ) {
52 fprintf(stderr, "\n message level must be positive\n") ;
53 exit(-1) ;
54 }
55 if ( strcmp(argv[2], "stdout") == 0 ) {
56 msgFile = stdout ;
57 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
58 fprintf(stderr, "\n unable to open file %s\n", argv[2]) ;
59 return(-1) ;
60 }
61 type = atoi(argv[3]) ;
62 mode = atoi(argv[4]) ;
63 nrowY = atoi(argv[5]) ;
64 ncolY = atoi(argv[6]) ;
65 nrowA = atoi(argv[7]) ;
66 ncolA = atoi(argv[8]) ;
67 nentA = atoi(argv[9]) ;
68 nrowX = atoi(argv[10]) ;
69 seed = atoi(argv[11]) ;
70 fprintf(msgFile, "\n %% %s:"
71 "\n %% msglvl = %d"
72 "\n %% msgFile = %s"
73 "\n %% type = %d"
74 "\n %% mode = %d"
75 "\n %% nrowY = %d"
76 "\n %% ncolY = %d"
77 "\n %% nrowA = %d"
78 "\n %% ncolA = %d"
79 "\n %% nentA = %d"
80 "\n %% nrowX = %d"
81 "\n %% seed = %d",
82 argv[0], msglvl, argv[2], type, mode, nrowY, ncolY,
83 nrowA, ncolA, nentA, nrowX, seed) ;
84 ncolX = ncolY ;
85 /*
86 -----------------------------
87 check for errors in the input
88 -----------------------------
89 */
90 if ( nrowA <= 0 || nrowA > nrowY
91 || ncolA <= 0 || ncolA > nrowX
92 || nentA > nrowA*ncolA
93 || nrowX <= 0 ) {
94 fprintf(stderr, "\n invalid input\n") ;
95 exit(-1) ;
96 }
97 switch ( type ) {
98 case SPOOLES_REAL :
99 case SPOOLES_COMPLEX :
100 break ;
101 default :
102 fprintf(stderr, "\n invalid type %d\n", type) ;
103 exit(-1) ;
104 }
105 switch ( mode ) {
106 case SUBMTX_DENSE_ROWS :
107 case SUBMTX_DENSE_COLUMNS :
108 case SUBMTX_SPARSE_ROWS :
109 case SUBMTX_SPARSE_COLUMNS :
110 break ;
111 default :
112 fprintf(stderr, "\n invalid mode %d\n", mode) ;
113 exit(-1) ;
114 }
115 /*
116 --------------------------------------
117 initialize the random number generator
118 --------------------------------------
119 */
120 drand = Drand_new() ;
121 Drand_init(drand) ;
122 Drand_setSeed(drand, seed) ;
123 Drand_setNormal(drand, 0.0, 1.0) ;
124 /*
125 ------------------------------
126 initialize the X SubMtx object
127 ------------------------------
128 */
129 MARKTIME(t1) ;
130 mtxX = SubMtx_new() ;
131 SubMtx_init(mtxX, type, SUBMTX_DENSE_COLUMNS, 0, 0,
132 nrowX, ncolX, nrowX*ncolX) ;
133 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
134 if ( SUBMTX_IS_REAL(mtxX) ) {
135 Drand_fillDvector(drand, nrowX*ncolX, entX) ;
136 } else if ( SUBMTX_IS_COMPLEX(mtxX) ) {
137 Drand_fillDvector(drand, 2*nrowX*ncolX, entX) ;
138 }
139 MARKTIME(t2) ;
140 fprintf(msgFile, "\n %% CPU : %.3f to initialize X SubMtx object",
141 t2 - t1) ;
142 /*
143 ------------------------------
144 initialize the Y SubMtx object
145 ------------------------------
146 */
147 MARKTIME(t1) ;
148 mtxY = SubMtx_new() ;
149 SubMtx_init(mtxY, type, SUBMTX_DENSE_COLUMNS, 0, 0,
150 nrowY, ncolY, nrowY*ncolY) ;
151 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
152 if ( SUBMTX_IS_REAL(mtxX) ) {
153 Drand_fillDvector(drand, nrowY*ncolY, entY) ;
154 DVzero(nrowY*ncolY, entY) ;
155 } else if ( SUBMTX_IS_COMPLEX(mtxX) ) {
156 Drand_fillDvector(drand, 2*nrowY*ncolY, entY) ;
157 DVzero(2*nrowY*ncolY, entY) ;
158 }
159 MARKTIME(t2) ;
160 fprintf(msgFile, "\n %% CPU : %.3f to initialize Y SubMtx object",
161 t2 - t1) ;
162 /*
163 -------------------------------------
164 initialize the A matrix SubMtx object
165 -------------------------------------
166 */
167 mtxA = SubMtx_new() ;
168 SubMtx_initRandom(mtxA, type, mode, 0, 0, nrowA, ncolA, nentA, seed) ;
169 /*
170 -------------------------
171 load the row indices of A
172 -------------------------
173 */
174 SubMtx_rowIndices(mtxA, &nrowA, &rowindA) ;
175 ivec = IVinit(nrowY, -1) ;
176 IVramp(nrowY, ivec, 0, 1) ;
177 IVshuffle(nrowY, ivec, seed+1) ;
178 IVcopy(nrowA, rowindA, ivec) ;
179 IVqsortUp(nrowA, rowindA) ;
180 IVfree(ivec) ;
181 if ( msglvl > 3 ) {
182 fprintf(msgFile, "\n %% row indices of A") ;
183 IVfprintf(msgFile, nrowA, rowindA) ;
184 fflush(msgFile) ;
185 }
186 /*
187 ----------------------------
188 load the column indices of A
189 ----------------------------
190 */
191 SubMtx_columnIndices(mtxA, &ncolA, &colindA) ;
192 ivec = IVinit(nrowX, -1) ;
193 IVramp(nrowX, ivec, 0, 1) ;
194 IVshuffle(nrowX, ivec, seed+2) ;
195 IVcopy(ncolA, colindA, ivec) ;
196 IVqsortUp(ncolA, colindA) ;
197 IVfree(ivec) ;
198 if ( msglvl > 3 ) {
199 fprintf(msgFile, "\n %% column indices of A") ;
200 IVfprintf(msgFile, ncolA, colindA) ;
201 fflush(msgFile) ;
202 }
203 /*
204 ----------------------------------
205 compute the matrix-matrix multiply
206 ----------------------------------
207 */
208 if ( type == SPOOLES_REAL ) {
209 double *colX, *pYij, *rowA ;
210 double sum ;
211 DV *colDV, *rowDV ;
212 int ii ;
213
214 ops = 2*nrowA*ncolA*ncolX ;
215 colDV = DV_new() ;
216 DV_init(colDV, nrowX, NULL) ;
217 colX = DV_entries(colDV) ;
218 rowDV = DV_new() ;
219 DV_init(rowDV, ncolA, NULL) ;
220 rowA = DV_entries(rowDV) ;
221 for ( jcolX = 0 ; jcolX < ncolX ; jcolX++ ) {
222 SubMtx_fillColumnDV(mtxX, jcolX, colDV) ;
223 for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
224 SubMtx_fillRowDV(mtxA, irowA, rowDV) ;
225 if ( ncolA == nrowX ) {
226 for ( ii = 0, sum = 0.0 ; ii < ncolA ; ii++ ) {
227 sum += rowA[ii] * colX[ii] ;
228 }
229 } else {
230 for ( ii = 0, sum = 0.0 ; ii < ncolA ; ii++ ) {
231 sum += rowA[ii] * colX[colindA[ii]] ;
232 }
233 }
234 if ( nrowA == nrowY ) {
235 SubMtx_locationOfRealEntry(mtxY, irowA, jcolX, &pYij) ;
236 } else {
237 SubMtx_locationOfRealEntry(mtxY, rowindA[irowA], jcolX,
238 &pYij) ;
239 }
240 *pYij = sum ;
241 }
242 }
243 DV_free(colDV) ;
244 DV_free(rowDV) ;
245 } else if ( type == SPOOLES_COMPLEX ) {
246 double *colX, *pYIij, *pYRij, *rowA ;
247 double idot, rdot ;
248 ZV *colZV, *rowZV ;
249
250 ops = 8*nrowA*ncolA*ncolX ;
251 colZV = ZV_new() ;
252 ZV_init(colZV, nrowX, NULL) ;
253 colX = ZV_entries(colZV) ;
254 rowZV = ZV_new() ;
255 ZV_init(rowZV, ncolA, NULL) ;
256 rowA = ZV_entries(rowZV) ;
257 for ( jcolX = 0 ; jcolX < ncolX ; jcolX++ ) {
258 SubMtx_fillColumnZV(mtxX, jcolX, colZV) ;
259 for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
260 SubMtx_fillRowZV(mtxA, irowA, rowZV) ;
261 if ( ncolA == nrowX ) {
262 ZVdotU(ncolA, colX, rowA, &rdot, &idot) ;
263 } else {
264 ZVdotiU(ncolA, colX, colindA, rowA, &rdot, &idot) ;
265 }
266 if ( nrowA == nrowY ) {
267 SubMtx_locationOfComplexEntry(mtxY,
268 irowA, jcolX, &pYRij, &pYIij) ;
269 } else {
270 SubMtx_locationOfComplexEntry(mtxY, rowindA[irowA], jcolX,
271 &pYRij, &pYIij) ;
272 }
273 *pYRij = rdot ;
274 *pYIij = idot ;
275 }
276 }
277 ZV_free(colZV) ;
278 ZV_free(rowZV) ;
279 }
280 MARKTIME(t2) ;
281 fprintf(msgFile, "\n %% CPU : %.3f to compute m-m, %.3f mflops",
282 t2 - t1, ops*1.e-6/(t2 - t1)) ;
283 if ( msglvl > 1 ) {
284 fprintf(msgFile, "\n\n %% Z SubMtx object") ;
285 fprintf(msgFile, "\n Z = zeros(%d,%d) ;", nrowY, ncolY) ;
286 SubMtx_writeForMatlab(mtxY, "Z", msgFile) ;
287 fflush(msgFile) ;
288 }
289 /*
290 ----------------------
291 print out the matrices
292 ----------------------
293 */
294 if ( msglvl > 1 ) {
295 fprintf(msgFile, "\n\n %% Y SubMtx object") ;
296 fprintf(msgFile, "\n Y = zeros(%d,%d) ;", nrowY, ncolY) ;
297 SubMtx_writeForMatlab(mtxY, "Y", msgFile) ;
298 fprintf(msgFile, "\n\n %% A SubMtx object") ;
299 fprintf(msgFile, "\n A = zeros(%d,%d) ;", nrowY, nrowX) ;
300 SubMtx_writeForMatlab(mtxA, "A", msgFile) ;
301 fprintf(msgFile, "\n\n %% X SubMtx object") ;
302 fprintf(msgFile, "\n X = zeros(%d,%d) ;", nrowX, ncolY) ;
303 SubMtx_writeForMatlab(mtxX, "X", msgFile) ;
304 fflush(msgFile) ;
305 }
306 /*
307 -----------------
308 check with matlab
309 -----------------
310 */
311 if ( msglvl > 1 ) {
312 fprintf(msgFile,
313 "\n\n emtx = abs(Y - A*X) ;"
314 "\n\n maxabs = max(max(emtx)) ") ;
315 fflush(msgFile) ;
316 }
317 /*
318 -------------------------------
319 compute the update Y := Y - A*X
320 (Y should now be zero)
321 -------------------------------
322 */
323 SubMtx_solveupd(mtxY, mtxA, mtxX) ;
324 /*
325 ----------------------
326 print out the Y matrix
327 ----------------------
328 */
329 if ( msglvl > 1 ) {
330 fprintf(msgFile, "\n\n %% Z SubMtx object") ;
331 fprintf(msgFile, "\n Z = zeros(%d,%d) ;", nrowY, ncolY) ;
332 SubMtx_writeForMatlab(mtxY, "Z", msgFile) ;
333 fflush(msgFile) ;
334 }
335 fprintf(msgFile, "\n RES %4d %4d %4d %4d %4d %4d %4d %4d %12.4e",
336 type, mode, nrowY, ncolY, nrowA, ncolA, nrowX, ncolX,
337 SubMtx_maxabs(mtxY)) ;
338 if ( msglvl > 1 ) {
339 fprintf(msgFile,
340 "\n\n emtx = abs(Z) ;"
341 "\n\n maxerr = max(max(emtx)) ") ;
342 fflush(msgFile) ;
343 }
344 /*
345 ------------------------
346 free the working storage
347 ------------------------
348 */
349 SubMtx_free(mtxA) ;
350 SubMtx_free(mtxX) ;
351 SubMtx_free(mtxY) ;
352 Drand_free(drand) ;
353
354 fprintf(msgFile, "\n") ;
355
356 return(1) ; }
357
358 /*--------------------------------------------------------------------*/
359