1 /* test_solveupdT.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_solveupdT() 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, mode, msglvl, ncolA, ncolAT,
24 nentA, nrowA, nrowAT, 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 ( type < 0 || type > 6
91 || nrowA <= 0 || nrowA > nrowX
92 || ncolA <= 0 || ncolA > nrowY
93 || nentA > nrowA*ncolA
94 || nrowX <= 0 ) {
95 fprintf(stderr, "\n invalid input\n") ;
96 exit(-1) ;
97 }
98 switch ( type ) {
99 case SPOOLES_REAL :
100 case SPOOLES_COMPLEX :
101 break ;
102 default :
103 fprintf(stderr, "\n invalid type %d\n", type) ;
104 exit(-1) ;
105 }
106 switch ( mode ) {
107 case SUBMTX_DENSE_ROWS :
108 case SUBMTX_DENSE_COLUMNS :
109 case SUBMTX_SPARSE_ROWS :
110 case SUBMTX_SPARSE_COLUMNS :
111 break ;
112 default :
113 fprintf(stderr, "\n invalid mode %d\n", mode) ;
114 exit(-1) ;
115 }
116 /*
117 --------------------------------------
118 initialize the random number generator
119 --------------------------------------
120 */
121 drand = Drand_new() ;
122 Drand_init(drand) ;
123 Drand_setSeed(drand, seed) ;
124 Drand_setNormal(drand, 0.0, 1.0) ;
125 /*
126 ----------------------------
127 initialize the X SubMtx object
128 ----------------------------
129 */
130 MARKTIME(t1) ;
131 mtxX = SubMtx_new() ;
132 SubMtx_init(mtxX, type, SUBMTX_DENSE_COLUMNS, 0, 0,
133 nrowX, ncolX, nrowX*ncolX) ;
134 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
135 if ( SUBMTX_IS_REAL(mtxX) ) {
136 Drand_fillDvector(drand, nrowX*ncolX, entX) ;
137 } else if ( SUBMTX_IS_COMPLEX(mtxX) ) {
138 Drand_fillDvector(drand, 2*nrowX*ncolX, entX) ;
139 }
140 MARKTIME(t2) ;
141 fprintf(msgFile, "\n %% CPU : %.3f to initialize X SubMtx object",
142 t2 - t1) ;
143 /*
144 ----------------------------
145 initialize the Y SubMtx object
146 ----------------------------
147 */
148 MARKTIME(t1) ;
149 mtxY = SubMtx_new() ;
150 SubMtx_init(mtxY, type, SUBMTX_DENSE_COLUMNS, 0, 0,
151 nrowY, ncolY, nrowY*ncolY) ;
152 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
153 if ( SUBMTX_IS_REAL(mtxX) ) {
154 DVzero(nrowY*ncolY, entY) ;
155 } else if ( SUBMTX_IS_COMPLEX(mtxX) ) {
156 DVzero(2*nrowY*ncolY, entY) ;
157 }
158 MARKTIME(t2) ;
159 fprintf(msgFile, "\n %% CPU : %.3f to initialize Y SubMtx object",
160 t2 - t1) ;
161 /*
162 -----------------------------------
163 initialize the A matrix SubMtx object
164 -----------------------------------
165 */
166 mtxA = SubMtx_new() ;
167 SubMtx_initRandom(mtxA, type, mode, 0, 0, nrowA, ncolA, nentA, seed) ;
168 /*
169 -------------------------
170 load the row indices of A
171 -------------------------
172 */
173 SubMtx_rowIndices(mtxA, &nrowA, &rowindA) ;
174 ivec = IVinit(nrowX, -1) ;
175 IVramp(nrowX, ivec, 0, 1) ;
176 IVshuffle(nrowX, ivec, seed+1) ;
177 IVcopy(nrowA, rowindA, ivec) ;
178 IVqsortUp(nrowA, rowindA) ;
179 IVfree(ivec) ;
180 if ( msglvl > 1 ) {
181 fprintf(msgFile, "\n %% row indices of A") ;
182 IVfprintf(msgFile, nrowA, rowindA) ;
183 fflush(msgFile) ;
184 }
185 /*
186 ----------------------------
187 load the column indices of A
188 ----------------------------
189 */
190 SubMtx_columnIndices(mtxA, &ncolA, &colindA) ;
191 ivec = IVinit(nrowY, -1) ;
192 IVramp(nrowY, ivec, 0, 1) ;
193 IVshuffle(nrowY, ivec, seed+2) ;
194 IVcopy(ncolA, colindA, ivec) ;
195 IVqsortUp(ncolA, colindA) ;
196 IVfree(ivec) ;
197 if ( msglvl > 1 ) {
198 fprintf(msgFile, "\n %% column indices of A") ;
199 IVfprintf(msgFile, ncolA, colindA) ;
200 fflush(msgFile) ;
201 }
202 /*
203 ----------------------------------
204 compute the matrix-matrix multiply
205 ----------------------------------
206 */
207 nrowAT = ncolA ;
208 ncolAT = nrowA ;
209 if ( type == SPOOLES_REAL ) {
210 double *colX, *pYij, *rowAT ;
211 double sum ;
212 DV *colDV, *rowDV ;
213 int ii, irowAT, jcolX ;
214
215 ops = 2*nrowA*ncolA*ncolX ;
216 colDV = DV_new() ;
217 DV_init(colDV, nrowX, NULL) ;
218 colX = DV_entries(colDV) ;
219 rowDV = DV_new() ;
220 DV_init(rowDV, ncolAT, NULL) ;
221 rowAT = DV_entries(rowDV) ;
222 for ( jcolX = 0 ; jcolX < ncolX ; jcolX++ ) {
223 SubMtx_fillColumnDV(mtxX, jcolX, colDV) ;
224 for ( irowAT = 0 ; irowAT < nrowAT ; irowAT++ ) {
225 SubMtx_fillColumnDV(mtxA, irowAT, rowDV) ;
226 if ( ncolAT == nrowX ) {
227 for ( ii = 0, sum = 0.0 ; ii < ncolAT ; ii++ ) {
228 sum += rowAT[ii] * colX[ii] ;
229 }
230 } else {
231 for ( ii = 0, sum = 0.0 ; ii < ncolAT ; ii++ ) {
232 sum += rowAT[ii] * colX[rowindA[ii]] ;
233 }
234 }
235 if ( nrowAT == nrowY ) {
236 SubMtx_locationOfRealEntry(mtxY, irowAT, jcolX, &pYij) ;
237 } else {
238 SubMtx_locationOfRealEntry(mtxY, colindA[irowAT], jcolX,
239 &pYij) ;
240 }
241 *pYij = sum ;
242 }
243 }
244 DV_free(colDV) ;
245 DV_free(rowDV) ;
246 } else if ( type == SPOOLES_COMPLEX ) {
247 double *colX, *pYIij, *pYRij, *rowAT ;
248 double idot, rdot ;
249 ZV *colZV, *rowZV ;
250 int irowAT, jcolX ;
251
252 ops = 8*nrowA*ncolA*ncolX ;
253 colZV = ZV_new() ;
254 ZV_init(colZV, nrowX, NULL) ;
255 colX = ZV_entries(colZV) ;
256 rowZV = ZV_new() ;
257 ZV_init(rowZV, ncolAT, NULL) ;
258 rowAT = ZV_entries(rowZV) ;
259 for ( jcolX = 0 ; jcolX < ncolX ; jcolX++ ) {
260 SubMtx_fillColumnZV(mtxX, jcolX, colZV) ;
261 for ( irowAT = 0 ; irowAT < nrowAT ; irowAT++ ) {
262 SubMtx_fillColumnZV(mtxA, irowAT, rowZV) ;
263 if ( ncolAT == nrowX ) {
264 ZVdotU(nrowA, colX, rowAT, &rdot, &idot) ;
265 } else {
266 ZVdotiU(nrowA, colX, rowindA, rowAT, &rdot, &idot) ;
267 }
268 if ( nrowAT == nrowY ) {
269 SubMtx_locationOfComplexEntry(mtxY,
270 irowAT, jcolX, &pYRij, &pYIij) ;
271 } else {
272 SubMtx_locationOfComplexEntry(mtxY, colindA[irowAT], jcolX,
273 &pYRij, &pYIij) ;
274 }
275 *pYRij = rdot ;
276 *pYIij = idot ;
277 }
278 }
279 ZV_free(colZV) ;
280 ZV_free(rowZV) ;
281 }
282 MARKTIME(t2) ;
283 fprintf(msgFile, "\n %% CPU : %.3f to compute m-m, %.3f mflops",
284 t2 - t1, ops*1.e-6/(t2 - t1)) ;
285 if ( msglvl > 1 ) {
286 fprintf(msgFile, "\n\n %% Z SubMtx object") ;
287 fprintf(msgFile, "\n Z = zeros(%d,%d) ;", nrowY, ncolY) ;
288 SubMtx_writeForMatlab(mtxY, "Z", msgFile) ;
289 fflush(msgFile) ;
290 }
291 /*
292 ----------------------
293 print out the matrices
294 ----------------------
295 */
296 if ( msglvl > 1 ) {
297 fprintf(msgFile, "\n\n %% Y SubMtx object") ;
298 fprintf(msgFile, "\n Y = zeros(%d,%d) ;", nrowY, ncolY) ;
299 SubMtx_writeForMatlab(mtxY, "Y", msgFile) ;
300 fprintf(msgFile, "\n\n %% A SubMtx object") ;
301 fprintf(msgFile, "\n A = zeros(%d,%d) ;", nrowX, nrowY) ;
302 SubMtx_writeForMatlab(mtxA, "A", msgFile) ;
303 fprintf(msgFile, "\n\n %% X SubMtx object") ;
304 fprintf(msgFile, "\n X = zeros(%d,%d) ;", nrowX, ncolY) ;
305 SubMtx_writeForMatlab(mtxX, "X", msgFile) ;
306 fflush(msgFile) ;
307 }
308 /*
309 -----------------
310 check with matlab
311 -----------------
312 */
313 if ( msglvl > 1 ) {
314 fprintf(msgFile,
315 "\n\n emtx = abs(Y - transpose(A)*X) ;"
316 "\n\n maxabs = max(max(emtx)) ") ;
317 fflush(msgFile) ;
318 }
319 /*
320 -------------------------------
321 compute the update Y := Y - A*X
322 (Y should now be zero)
323 -------------------------------
324 */
325 SubMtx_solveupdT(mtxY, mtxA, mtxX) ;
326 /*
327 ----------------------
328 print out the Y matrix
329 ----------------------
330 */
331 if ( msglvl > 1 ) {
332 fprintf(msgFile, "\n\n %% Y SubMtx object") ;
333 fprintf(msgFile, "\n Z = zeros(%d,%d) ;", nrowY, ncolY) ;
334 SubMtx_writeForMatlab(mtxY, "Z", msgFile) ;
335 fflush(msgFile) ;
336 }
337 fprintf(msgFile, "\n %% RES %4d %4d %4d %4d %4d %4d %4d %4d %12.4e",
338 type, mode, nrowY, ncolY, nrowA, ncolA, nrowX, ncolX,
339 SubMtx_maxabs(mtxY)) ;
340 if ( msglvl > 1 ) {
341 fprintf(msgFile,
342 "\n\n emtx = abs(Z) ;"
343 "\n\n maxerr = max(max(emtx)) ") ;
344 fflush(msgFile) ;
345 }
346 /*
347 ------------------------
348 free the working storage
349 ------------------------
350 */
351 SubMtx_free(mtxA) ;
352 SubMtx_free(mtxX) ;
353 SubMtx_free(mtxY) ;
354 Drand_free(drand) ;
355
356 fprintf(msgFile, "\n") ;
357
358 return(1) ; }
359
360 /*--------------------------------------------------------------------*/
361