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