1 /* test_factorupd.c */
2
3 #include "../Chv.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 Chv_update{H,S,N}() methods.
13 T := T - U^T * D * U
14 T := T - U^H * D * U
15 T := T - L * D * U
16
17 created -- 98apr23, cca
18 -------------------------------------
19 */
20 {
21 Chv *chvT ;
22 SubMtx *mtxD, *mtxL, *mtxU ;
23 double imag, ops, real, t1, t2 ;
24 Drand *drand ;
25 DV *tempDV ;
26 FILE *msgFile ;
27 int irow, msglvl, ncolT, nDT, ncolU, nentT, nentU, nrowD,
28 nrowL, nrowT, offset, seed, size, sparsityflag, symflag, type ;
29 int *colindT, *colindU, *ivec, *rowindL, *rowindT ;
30
31 if ( argc != 13 ) {
32 fprintf(stdout,
33 "\n\n usage : %s msglvl msgFile type symflag sparsityflag"
34 "\n ncolT ncolU nrowD nentU offset seed"
35 "\n msglvl -- message level"
36 "\n msgFile -- message file"
37 "\n type -- entries type"
38 "\n 1 -- real"
39 "\n 2 -- complex"
40 "\n symflag -- type of matrix U"
41 "\n 0 -- symmetric"
42 "\n 1 -- hermitian"
43 "\n 2 -- nonsymmetric"
44 "\n sparsityflag -- dense or sparse"
45 "\n 0 -- dense"
46 "\n 1 -- sparse"
47 "\n ncolT -- # of rows and columns in matrix T"
48 "\n nDT -- # of internal rows and columns in matrix T"
49 "\n ncolU -- # of rows and columns in matrix U"
50 "\n nrowD -- # of rows and columns in matrix D"
51 "\n nentU -- # of entries in matrix U"
52 "\n offset -- distance between D_I and T"
53 "\n seed -- random number seed"
54 "\n", argv[0]) ;
55 return(0) ;
56 }
57 if ( (msglvl = atoi(argv[1])) < 0 ) {
58 fprintf(stderr, "\n message level must be positive\n") ;
59 exit(-1) ;
60 }
61 if ( strcmp(argv[2], "stdout") == 0 ) {
62 msgFile = stdout ;
63 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
64 fprintf(stderr, "\n unable to open file %s\n", argv[2]) ;
65 return(-1) ;
66 }
67 type = atoi(argv[3]) ;
68 symflag = atoi(argv[4]) ;
69 sparsityflag = atoi(argv[5]) ;
70 ncolT = atoi(argv[6]) ;
71 nDT = atoi(argv[7]) ;
72 ncolU = atoi(argv[8]) ;
73 nrowD = atoi(argv[9]) ;
74 nentU = atoi(argv[10]) ;
75 offset = atoi(argv[11]) ;
76 seed = atoi(argv[12]) ;
77 fprintf(msgFile, "\n %% %s:"
78 "\n %% msglvl = %d"
79 "\n %% msgFile = %s"
80 "\n %% type = %d"
81 "\n %% symflag = %d"
82 "\n %% sparsityflag = %d"
83 "\n %% ncolT = %d"
84 "\n %% nDT = %d"
85 "\n %% ncolU = %d"
86 "\n %% nrowD = %d"
87 "\n %% nentU = %d"
88 "\n %% offset = %d"
89 "\n %% seed = %d",
90 argv[0], msglvl, argv[2], type, symflag, sparsityflag,
91 ncolT, nDT, ncolU, nrowD, nentU, offset, seed) ;
92 /*
93 -----------------------------
94 check for errors in the input
95 -----------------------------
96 */
97 if ( (type != SPOOLES_REAL
98 && type != SPOOLES_COMPLEX)
99 || (symflag != SPOOLES_SYMMETRIC
100 && symflag != SPOOLES_HERMITIAN
101 && symflag != SPOOLES_NONSYMMETRIC)
102 || (sparsityflag < 0 || sparsityflag > 1)
103 || ncolT <= 0 || ncolU > (ncolT + offset) || nrowD <= 0 ) {
104 fprintf(stderr, "\n invalid input\n") ;
105 exit(-1) ;
106 }
107 /*
108 --------------------------------------
109 initialize the random number generator
110 --------------------------------------
111 */
112 drand = Drand_new() ;
113 Drand_init(drand) ;
114 Drand_setSeed(drand, ++seed) ;
115 Drand_setNormal(drand, 0.0, 1.0) ;
116 /*
117 -----------------------
118 get a vector of indices
119 -----------------------
120 */
121 size = nrowD + offset + ncolT ;
122 ivec = IVinit(size, -1) ;
123 IVramp(size, ivec, 0, 1) ;
124 /*
125 ----------------------------
126 initialize the T Chv object
127 ----------------------------
128 */
129 fprintf(msgFile, "\n\n %% symflag = %d", symflag) ;
130 MARKTIME(t1) ;
131 chvT = Chv_new() ;
132 Chv_init(chvT, 0, nDT, ncolT - nDT, ncolT - nDT, type, symflag) ;
133 nentT = Chv_nent(chvT) ;
134 if ( CHV_IS_REAL(chvT) ) {
135 Drand_fillDvector(drand, nentT, Chv_entries(chvT)) ;
136 } else if ( CHV_IS_COMPLEX(chvT) ) {
137 Drand_fillDvector(drand, 2*nentT, Chv_entries(chvT)) ;
138 }
139 Chv_columnIndices(chvT, &ncolT, &colindT) ;
140 IVcopy(ncolT, colindT, ivec + nrowD + offset) ;
141 if ( CHV_IS_NONSYMMETRIC(chvT) ) {
142 Chv_rowIndices(chvT, &nrowT, &rowindT) ;
143 IVcopy(nrowT, rowindT, colindT) ;
144 }
145 IVfree(ivec) ;
146 if ( CHV_IS_HERMITIAN(chvT) ) {
147 fprintf(msgFile, "\n\n %% hermitian\n") ;
148 /*
149 ---------------------------------------------------------
150 hermitian example, set imaginary part of diagonal to zero
151 ---------------------------------------------------------
152 */
153 for ( irow = 0 ; irow < nDT ; irow++ ) {
154 Chv_complexEntry(chvT, irow, irow, &real, &imag) ;
155 Chv_setComplexEntry(chvT, irow, irow, real, 0.0) ;
156 }
157 }
158 MARKTIME(t2) ;
159 fprintf(msgFile, "\n %% CPU : %.3f to initialize chvT Chv object",
160 t2 - t1) ;
161 fprintf(msgFile, "\n T = zeros(%d,%d); ", size, size) ;
162 Chv_writeForMatlab(chvT, "T", msgFile) ;
163 /*
164 ---------------------------
165 initialize the D Mtx object
166 ---------------------------
167 */
168 MARKTIME(t1) ;
169 mtxD = SubMtx_new() ;
170 if ( CHV_IS_REAL(chvT) ) {
171 if ( CHV_IS_SYMMETRIC(chvT) ) {
172 SubMtx_initRandom(mtxD, SPOOLES_REAL, SUBMTX_BLOCK_DIAGONAL_SYM,
173 0, 0, nrowD, nrowD, nrowD*nrowD, ++seed) ;
174 } else {
175 SubMtx_initRandom(mtxD, SPOOLES_REAL, SUBMTX_DIAGONAL,
176 0, 0, nrowD, nrowD, nrowD*nrowD, ++seed) ;
177 }
178 } else if ( CHV_IS_COMPLEX(chvT) ) {
179 if ( CHV_IS_HERMITIAN(chvT) ) {
180 SubMtx_initRandom(mtxD,SPOOLES_COMPLEX,SUBMTX_BLOCK_DIAGONAL_HERM,
181 0, 0, nrowD, nrowD, nrowD*nrowD, ++seed) ;
182 } else if ( CHV_IS_SYMMETRIC(chvT) ) {
183 SubMtx_initRandom(mtxD,SPOOLES_COMPLEX, SUBMTX_BLOCK_DIAGONAL_SYM,
184 0, 0, nrowD, nrowD, nrowD*nrowD, ++seed) ;
185 } else {
186 SubMtx_initRandom(mtxD, SPOOLES_COMPLEX, SUBMTX_DIAGONAL,
187 0, 0, nrowD, nrowD, nrowD*nrowD, ++seed) ;
188 }
189 }
190 MARKTIME(t2) ;
191 fprintf(msgFile, "\n %% CPU : %.3f to initialize D SubMtx object",
192 t2 - t1) ;
193 fprintf(msgFile, "\n D = zeros(%d,%d) ;", nrowD, nrowD) ;
194 SubMtx_writeForMatlab(mtxD, "D", msgFile) ;
195 /*
196 ----------------------------
197 initialize the U SubMtx object
198 ----------------------------
199 */
200 MARKTIME(t1) ;
201 mtxU = SubMtx_new() ;
202 if ( CHV_IS_REAL(chvT) ) {
203 if ( sparsityflag == 0 ) {
204 SubMtx_initRandom(mtxU, SPOOLES_REAL, SUBMTX_DENSE_COLUMNS,
205 0, 0, nrowD, ncolU, nentU, ++seed) ;
206 } else {
207 SubMtx_initRandom(mtxU, SPOOLES_REAL, SUBMTX_SPARSE_COLUMNS,
208 0, 0, nrowD, ncolU, nentU, ++seed) ;
209 }
210 } else if ( CHV_IS_COMPLEX(chvT) ) {
211 if ( sparsityflag == 0 ) {
212 SubMtx_initRandom(mtxU, SPOOLES_COMPLEX, SUBMTX_DENSE_COLUMNS,
213 0, 0, nrowD, ncolU, nentU, ++seed) ;
214 } else {
215 SubMtx_initRandom(mtxU, SPOOLES_COMPLEX, SUBMTX_SPARSE_COLUMNS,
216 0, 0, nrowD, ncolU, nentU, ++seed) ;
217 }
218 }
219 ivec = IVinit(offset + ncolT, -1) ;
220 IVramp(offset + ncolT, ivec, nrowD, 1) ;
221 IVshuffle(offset + ncolT, ivec, ++seed) ;
222 SubMtx_columnIndices(mtxU, &ncolU, &colindU) ;
223 IVcopy(ncolU, colindU, ivec) ;
224 IVqsortUp(ncolU, colindU) ;
225 IVfree(ivec) ;
226 MARKTIME(t2) ;
227 fprintf(msgFile, "\n %% CPU : %.3f to initialize U SubMtx object",
228 t2 - t1) ;
229 fprintf(msgFile, "\n U = zeros(%d,%d) ;", nrowD, size) ;
230 SubMtx_writeForMatlab(mtxU, "U", msgFile) ;
231 if ( CHV_IS_NONSYMMETRIC(chvT) ) {
232 /*
233 ----------------------------
234 initialize the L SubMtx object
235 ----------------------------
236 */
237 MARKTIME(t1) ;
238 mtxL = SubMtx_new() ;
239 if ( CHV_IS_REAL(chvT) ) {
240 if ( sparsityflag == 0 ) {
241 SubMtx_initRandom(mtxL, SPOOLES_REAL, SUBMTX_DENSE_ROWS,
242 0, 0, ncolU, nrowD, nentU, ++seed) ;
243 } else {
244 SubMtx_initRandom(mtxL, SPOOLES_REAL, SUBMTX_SPARSE_ROWS,
245 0, 0, ncolU, nrowD, nentU, ++seed) ;
246 }
247 } else if ( CHV_IS_COMPLEX(chvT) ) {
248 if ( sparsityflag == 0 ) {
249 SubMtx_initRandom(mtxL, SPOOLES_COMPLEX, SUBMTX_DENSE_ROWS,
250 0, 0, ncolU, nrowD, nentU, ++seed) ;
251 } else {
252 SubMtx_initRandom(mtxL, SPOOLES_COMPLEX, SUBMTX_SPARSE_ROWS,
253 0, 0, ncolU, nrowD, nentU, ++seed) ;
254 }
255 }
256 SubMtx_rowIndices(mtxL, &nrowL, &rowindL) ;
257 IVcopy(nrowL, rowindL, colindU) ;
258 MARKTIME(t2) ;
259 fprintf(msgFile, "\n %% CPU : %.3f to initialize L SubMtx object",
260 t2 - t1) ;
261 fprintf(msgFile, "\n L = zeros(%d,%d) ;", size, nrowD) ;
262 SubMtx_writeForMatlab(mtxL, "L", msgFile) ;
263 } else {
264 mtxL = NULL ;
265 }
266 /*
267 --------------------------------
268 compute the matrix-matrix update
269 --------------------------------
270 */
271 tempDV = DV_new() ;
272 ops = 8*nrowD*nrowD*ncolU ;
273 if ( CHV_IS_SYMMETRIC(chvT) ) {
274 Chv_updateS(chvT, mtxD, mtxU, tempDV) ;
275 } else if ( CHV_IS_HERMITIAN(chvT) ) {
276 Chv_updateH(chvT, mtxD, mtxU, tempDV) ;
277 } else if ( CHV_IS_NONSYMMETRIC(chvT) ) {
278 Chv_updateN(chvT, mtxL, mtxD, mtxU, tempDV) ;
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 Chv object") ;
285 fprintf(msgFile, "\n Z = zeros(%d,%d); ", size, size) ;
286 Chv_writeForMatlab(chvT, "Z", msgFile) ;
287 fflush(msgFile) ;
288 }
289 /*
290 -----------------
291 check with matlab
292 -----------------
293 */
294 if ( msglvl > 1 ) {
295 if ( CHV_IS_HERMITIAN(chvT) ) {
296 fprintf(msgFile, "\n\n B = ctranspose(U) * D * U ;") ;
297 } else if ( CHV_IS_SYMMETRIC(chvT) ) {
298 fprintf(msgFile, "\n\n B = transpose(U) * D * U ;") ;
299 } else {
300 fprintf(msgFile, "\n\n B = L * D * U ;") ;
301 }
302 fprintf(msgFile,
303 "\n\n for irow = 1:%d"
304 "\n for jcol = 1:%d"
305 "\n if T(irow,jcol) ~= 0.0"
306 "\n T(irow,jcol) = T(irow,jcol) - B(irow,jcol) ;"
307 "\n end"
308 "\n end"
309 "\n end"
310 "\n emtx = abs(Z - T) ;",
311 size, size) ;
312 fprintf(msgFile, "\n\n maxabs = max(max(emtx)) ") ;
313 fflush(msgFile) ;
314 }
315 /*
316 ------------------------
317 free the working storage
318 ------------------------
319 */
320 if ( mtxL != NULL ) {
321 SubMtx_free(mtxL) ;
322 }
323 Chv_free(chvT) ;
324 SubMtx_free(mtxD) ;
325 SubMtx_free(mtxU) ;
326 DV_free(tempDV) ;
327 Drand_free(drand) ;
328
329 fprintf(msgFile, "\n") ;
330
331 return(1) ; }
332
333 /*--------------------------------------------------------------------*/
334