1 /* allInOneMPI.c */
2
3 #include "../spoolesMPI.h"
4 #include "../../timings.h"
5
6 /*--------------------------------------------------------------------*/
7 int
main(int argc,char * argv[])8 main ( int argc, char *argv[] ) {
9 /*
10 ------------------------------------------------------------
11 all-in-one MPI program for each process
12
13 order, factor and solve A X = Y
14
15 ( 1) read in matrix entries and form InpMtx object for A
16 ( 2) order the system using minimum degree
17 ( 3) permute the front tree
18 ( 4) create the owners map IV object
19 ( 5) permute the matrix A and redistribute
20 ( 6) compute the symbolic factorization
21 ( 7) compute the numeric factorization
22 ( 8) split the factors into submatrices
23 ( 9) create the submatrix map and redistribute
24 (10) read in right hand side entries
25 and form dense matrix DenseMtx object for Y
26 (11) permute and redistribute Y
27 (12) solve the linear system
28 (13) gather X on processor 0
29
30 created -- 98jun13, cca
31 ------------------------------------------------------------
32 */
33 /*--------------------------------------------------------------------*/
34 char buffer[128] ;
35 Chv *rootchv ;
36 ChvManager *chvmanager ;
37 DenseMtx *mtxX, *mtxY, *newY ;
38 SubMtxManager *mtxmanager, *solvemanager ;
39 FrontMtx *frontmtx ;
40 InpMtx *mtxA, *newA ;
41 double cutoff, droptol = 0.0, minops, tau = 100. ;
42 double cpus[20] ;
43 double *opcounts ;
44 DV *cumopsDV ;
45 ETree *frontETree ;
46 FILE *inputFile, *msgFile ;
47 Graph *graph ;
48 int error, firsttag, ient, irow, jcol, lookahead = 0,
49 msglvl, myid, nedges, nent, neqns, nmycol, nproc, nrhs,
50 nrow, pivotingflag, root, seed, symmetryflag, type ;
51 int stats[20] ;
52 int *rowind ;
53 IV *oldToNewIV, *ownedColumnsIV, *ownersIV,
54 *newToOldIV, *vtxmapIV ;
55 IVL *adjIVL, *symbfacIVL ;
56 SolveMap *solvemap ;
57 /*--------------------------------------------------------------------*/
58 /*
59 ---------------------------------------------------------------
60 find out the identity of this process and the number of process
61 ---------------------------------------------------------------
62 */
63 MPI_Init(&argc, &argv) ;
64 MPI_Comm_rank(MPI_COMM_WORLD, &myid) ;
65 MPI_Comm_size(MPI_COMM_WORLD, &nproc) ;
66 /*--------------------------------------------------------------------*/
67 /*
68 --------------------
69 get input parameters
70 --------------------
71 */
72 if ( argc != 7 ) {
73 fprintf(stdout,
74 "\n usage: %s msglvl msgFile type symmetryflag pivotingflag seed"
75 "\n msglvl -- message level"
76 "\n msgFile -- message file"
77 "\n type -- type of entries"
78 "\n 1 (SPOOLES_REAL) -- real entries"
79 "\n 2 (SPOOLES_COMPLEX) -- complex entries"
80 "\n symmetryflag -- type of matrix"
81 "\n 0 (SPOOLES_SYMMETRIC) -- symmetric entries"
82 "\n 1 (SPOOLES_HERMITIAN) -- Hermitian entries"
83 "\n 2 (SPOOLES_NONSYMMETRIC) -- nonsymmetric entries"
84 "\n pivotingflag -- type of pivoting"
85 "\n 0 (SPOOLES_NO_PIVOTING) -- no pivoting used"
86 "\n 1 (SPOOLES_PIVOTING) -- pivoting used"
87 "\n seed -- random number seed"
88 "\n "
89 "\n note: matrix entries are read in from matrix.k.input"
90 "\n where k is the process number"
91 "\n note: rhs entries are read in from rhs.k.input"
92 "\n where k is the process number"
93 "\n", argv[0]) ;
94 return(0) ;
95 }
96 msglvl = atoi(argv[1]) ;
97 if ( strcmp(argv[2], "stdout") == 0 ) {
98 msgFile = stdout ;
99 } else {
100 sprintf(buffer, "res.%d", myid) ;
101 if ( (msgFile = fopen(buffer, "w")) == NULL ) {
102 fprintf(stderr, "\n fatal error in %s"
103 "\n unable to open file %s\n",
104 argv[0], buffer) ;
105 return(-1) ;
106 }
107 }
108 type = atoi(argv[3]) ;
109 symmetryflag = atoi(argv[4]) ;
110 pivotingflag = atoi(argv[5]) ;
111 seed = atoi(argv[6]) ;
112 IVzero(20, stats) ;
113 DVzero(20, cpus) ;
114 fprintf(msgFile,
115 "\n\n input data"
116 "\n msglvl = %d"
117 "\n msgFile = %s"
118 "\n type = %d"
119 "\n symmetryflag = %d"
120 "\n pivotingflag = %d"
121 "\n seed = %d",
122 msglvl, argv[2], type, symmetryflag, pivotingflag, seed) ;
123 fflush(msgFile) ;
124 MPI_Barrier(MPI_COMM_WORLD) ;
125 /*--------------------------------------------------------------------*/
126 /*
127 --------------------------------------------
128 STEP 1: read the entries from the input file
129 and create the InpMtx object
130 --------------------------------------------
131 */
132 sprintf(buffer, "haggar.mtx.%d.input", myid) ;
133 inputFile = fopen(buffer, "r") ;
134 fscanf(inputFile, "%d %d %d", &neqns, &neqns, &nent) ;
135 MPI_Barrier(MPI_COMM_WORLD) ;
136 mtxA = InpMtx_new() ;
137 InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nent, 0) ;
138 if ( type == SPOOLES_REAL ) {
139 double value ;
140 for ( ient = 0 ; ient < nent ; ient++ ) {
141 fscanf(inputFile, "%d %d %le", &irow, &jcol, &value) ;
142 InpMtx_inputRealEntry(mtxA, irow, jcol, value) ;
143 }
144 } else if ( type == SPOOLES_COMPLEX ) {
145 double imag, real ;
146 for ( ient = 0 ; ient < nent ; ient++ ) {
147 fscanf(inputFile, "%d %d %le %le", &irow, &jcol, &real, &imag) ;
148 InpMtx_inputComplexEntry(mtxA, irow, jcol, real, imag) ;
149 }
150 }
151 fclose(inputFile) ;
152 InpMtx_sortAndCompress(mtxA) ;
153 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
154 if ( msglvl > 2 ) {
155 fprintf(msgFile, "\n\n input matrix") ;
156 InpMtx_writeForHumanEye(mtxA, msgFile) ;
157 fflush(msgFile) ;
158 }
159 /*--------------------------------------------------------------------*/
160 /*
161 ----------------------------------------------------
162 STEP 2: read the rhs entries from the rhs input file
163 and create the DenseMtx object for Y
164 ----------------------------------------------------
165 */
166 sprintf(buffer, "haggar.rhs.%d.input", myid) ;
167 inputFile = fopen(buffer, "r") ;
168 fscanf(inputFile, "%d %d", &nrow, &nrhs) ;
169 mtxY = DenseMtx_new() ;
170 DenseMtx_init(mtxY, type, 0, 0, nrow, nrhs, 1, nrow) ;
171 DenseMtx_rowIndices(mtxY, &nrow, &rowind) ;
172 if ( type == SPOOLES_REAL ) {
173 double value ;
174 for ( irow = 0 ; irow < nrow ; irow++ ) {
175 fscanf(inputFile, "%d", rowind + irow) ;
176 for ( jcol = 0 ; jcol < nrhs ; jcol++ ) {
177 fscanf(inputFile, "%le", &value) ;
178 DenseMtx_setRealEntry(mtxY, irow, jcol, value) ;
179 }
180 }
181 } if ( type == SPOOLES_COMPLEX ) {
182 double imag, real ;
183 for ( irow = 0 ; irow < nrow ; irow++ ) {
184 fscanf(inputFile, "%d", rowind + irow) ;
185 for ( jcol = 0 ; jcol < nrhs ; jcol++ ) {
186 fscanf(inputFile, "%le %le", &real, &imag) ;
187 DenseMtx_setComplexEntry(mtxY, irow, jcol, real, imag) ;
188 }
189 }
190 }
191 fclose(inputFile) ;
192 if ( msglvl > 2 ) {
193 fprintf(msgFile, "\n\n rhs matrix in original ordering") ;
194 DenseMtx_writeForHumanEye(mtxY, msgFile) ;
195 fflush(msgFile) ;
196 }
197 /*--------------------------------------------------------------------*/
198 /*
199 -------------------------------------------------------
200 STEP 2 : find a low-fill ordering
201 (1) create the Graph object
202 (2) order the graph using multiple minimum degree
203 (3) find out who has the best ordering w.r.t. op count,
204 and broadcast that front tree object
205 -------------------------------------------------------
206 */
207 graph = Graph_new() ;
208 adjIVL = InpMtx_MPI_fullAdjacency(mtxA, stats,
209 msglvl, msgFile, MPI_COMM_WORLD) ;
210 nedges = IVL_tsize(adjIVL) ;
211 Graph_init2(graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL,
212 NULL, NULL) ;
213 if ( msglvl > 2 ) {
214 fprintf(msgFile, "\n\n graph of the input matrix") ;
215 Graph_writeForHumanEye(graph, msgFile) ;
216 fflush(msgFile) ;
217 }
218 frontETree = orderViaMMD(graph, seed + myid, msglvl, msgFile) ;
219 Graph_free(graph) ;
220 if ( msglvl > 2 ) {
221 fprintf(msgFile, "\n\n front tree from ordering") ;
222 ETree_writeForHumanEye(frontETree, msgFile) ;
223 fflush(msgFile) ;
224 }
225 opcounts = DVinit(nproc, 0.0) ;
226 opcounts[myid] = ETree_nFactorOps(frontETree, type, symmetryflag) ;
227 MPI_Allgather((void *) &opcounts[myid], 1, MPI_DOUBLE,
228 (void *) opcounts, 1, MPI_DOUBLE, MPI_COMM_WORLD) ;
229 minops = DVmin(nproc, opcounts, &root) ;
230 DVfree(opcounts) ;
231 frontETree = ETree_MPI_Bcast(frontETree, root,
232 msglvl, msgFile, MPI_COMM_WORLD) ;
233 if ( msglvl > 2 ) {
234 fprintf(msgFile, "\n\n best front tree") ;
235 ETree_writeForHumanEye(frontETree, msgFile) ;
236 fflush(msgFile) ;
237 }
238 /*--------------------------------------------------------------------*/
239 /*
240 -------------------------------------------------------
241 STEP 3: get the permutations, permute the front tree,
242 permute the matrix and right hand side.
243 -------------------------------------------------------
244 */
245 oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
246 newToOldIV = ETree_newToOldVtxPerm(frontETree) ;
247 ETree_permuteVertices(frontETree, oldToNewIV) ;
248 InpMtx_permute(mtxA, IV_entries(oldToNewIV), IV_entries(oldToNewIV)) ;
249 if ( symmetryflag == SPOOLES_SYMMETRIC
250 || symmetryflag == SPOOLES_HERMITIAN ) {
251 InpMtx_mapToUpperTriangle(mtxA) ;
252 }
253 InpMtx_changeCoordType(mtxA, INPMTX_BY_CHEVRONS) ;
254 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
255 DenseMtx_permuteRows(mtxY, oldToNewIV) ;
256 if ( msglvl > 2 ) {
257 fprintf(msgFile, "\n\n rhs matrix in new ordering") ;
258 DenseMtx_writeForHumanEye(mtxY, msgFile) ;
259 fflush(msgFile) ;
260 }
261 /*--------------------------------------------------------------------*/
262 /*
263 -------------------------------------------
264 STEP 4: generate the owners map IV object
265 and the map from vertices to owners
266 -------------------------------------------
267 */
268 cutoff = 1./(2*nproc) ;
269 cumopsDV = DV_new() ;
270 DV_init(cumopsDV, nproc, NULL) ;
271 ownersIV = ETree_ddMap(frontETree,
272 type, symmetryflag, cumopsDV, cutoff) ;
273 DV_free(cumopsDV) ;
274 vtxmapIV = IV_new() ;
275 IV_init(vtxmapIV, neqns, NULL) ;
276 IVgather(neqns, IV_entries(vtxmapIV),
277 IV_entries(ownersIV), ETree_vtxToFront(frontETree)) ;
278 if ( msglvl > 2 ) {
279 fprintf(msgFile, "\n\n map from fronts to owning processes") ;
280 IV_writeForHumanEye(ownersIV, msgFile) ;
281 fprintf(msgFile, "\n\n map from vertices to owning processes") ;
282 IV_writeForHumanEye(vtxmapIV, msgFile) ;
283 fflush(msgFile) ;
284 }
285 if ( myid == 0 ) {
286 IV_writeToFile(ownersIV, "../../Tree/drivers/haggar.ivf") ;
287 }
288 /*--------------------------------------------------------------------*/
289 /*
290 ---------------------------------------------------
291 STEP 5: redistribute the matrix and right hand side
292 ---------------------------------------------------
293 */
294 firsttag = 0 ;
295 newA = InpMtx_MPI_split(mtxA, vtxmapIV, stats,
296 msglvl, msgFile, firsttag, MPI_COMM_WORLD) ;
297 firsttag++ ;
298 InpMtx_free(mtxA) ;
299 mtxA = newA ;
300 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
301 if ( msglvl > 2 ) {
302 fprintf(msgFile, "\n\n split InpMtx") ;
303 InpMtx_writeForHumanEye(mtxA, msgFile) ;
304 fflush(msgFile) ;
305 }
306 newY = DenseMtx_MPI_splitByRows(mtxY, vtxmapIV, stats, msglvl,
307 msgFile, firsttag, MPI_COMM_WORLD) ;
308 DenseMtx_free(mtxY) ;
309 mtxY = newY ;
310 firsttag += nproc ;
311 if ( msglvl > 2 ) {
312 fprintf(msgFile, "\n\n split DenseMtx Y") ;
313 DenseMtx_writeForHumanEye(mtxY, msgFile) ;
314 fflush(msgFile) ;
315 }
316 /*--------------------------------------------------------------------*/
317 /*
318 ------------------------------------------
319 STEP 6: compute the symbolic factorization
320 ------------------------------------------
321 */
322 symbfacIVL = SymbFac_MPI_initFromInpMtx(frontETree, ownersIV, mtxA,
323 stats, msglvl, msgFile, firsttag, MPI_COMM_WORLD) ;
324 firsttag += frontETree->nfront ;
325 if ( msglvl > 2 ) {
326 fprintf(msgFile, "\n\n local symbolic factorization") ;
327 IVL_writeForHumanEye(symbfacIVL, msgFile) ;
328 fflush(msgFile) ;
329 }
330 /*--------------------------------------------------------------------*/
331 /*
332 -----------------------------------
333 STEP 7: initialize the front matrix
334 -----------------------------------
335 */
336 mtxmanager = SubMtxManager_new() ;
337 SubMtxManager_init(mtxmanager, NO_LOCK, 0) ;
338 frontmtx = FrontMtx_new() ;
339 FrontMtx_init(frontmtx, frontETree, symbfacIVL, type, symmetryflag,
340 FRONTMTX_DENSE_FRONTS, pivotingflag, NO_LOCK, myid,
341 ownersIV, mtxmanager, msglvl, msgFile) ;
342 /*--------------------------------------------------------------------*/
343 /*
344 ---------------------------------
345 STEP 8: compute the factorization
346 ---------------------------------
347 */
348 chvmanager = ChvManager_new() ;
349 ChvManager_init(chvmanager, NO_LOCK, 0) ;
350 rootchv = FrontMtx_MPI_factorInpMtx(frontmtx, mtxA, tau, droptol,
351 chvmanager, ownersIV, lookahead, &error, cpus,
352 stats, msglvl, msgFile, firsttag, MPI_COMM_WORLD) ;
353 ChvManager_free(chvmanager) ;
354 firsttag += 3*frontETree->nfront + 2 ;
355 if ( msglvl > 2 ) {
356 fprintf(msgFile, "\n\n numeric factorization") ;
357 FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
358 fflush(msgFile) ;
359 }
360 if ( error >= 0 ) {
361 fprintf(stderr,
362 "\n proc %d : factorization error at front %d", myid, error) ;
363 MPI_Finalize() ;
364 exit(-1) ;
365 }
366 /*--------------------------------------------------------------------*/
367 /*
368 ------------------------------------------------
369 STEP 9: post-process the factorization and split
370 the factor matrices into submatrices
371 ------------------------------------------------
372 */
373 FrontMtx_MPI_postProcess(frontmtx, ownersIV, stats, msglvl,
374 msgFile, firsttag, MPI_COMM_WORLD) ;
375 firsttag += 5*nproc ;
376 if ( msglvl > 2 ) {
377 fprintf(msgFile, "\n\n numeric factorization after post-processing");
378 FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
379 fflush(msgFile) ;
380 }
381 /*--------------------------------------------------------------------*/
382 /*
383 -----------------------------------
384 STEP 10: create the solve map object
385 -----------------------------------
386 */
387 solvemap = SolveMap_new() ;
388 SolveMap_ddMap(solvemap, frontmtx->symmetryflag,
389 FrontMtx_upperBlockIVL(frontmtx),
390 FrontMtx_lowerBlockIVL(frontmtx),
391 nproc, ownersIV, FrontMtx_frontTree(frontmtx),
392 seed, msglvl, msgFile);
393 if ( msglvl > 2 ) {
394 SolveMap_writeForHumanEye(solvemap, msgFile) ;
395 fflush(msgFile) ;
396 }
397 /*--------------------------------------------------------------------*/
398 /*
399 ----------------------------------------------------
400 STEP 11: redistribute the submatrices of the factors
401 ----------------------------------------------------
402 */
403 FrontMtx_MPI_split(frontmtx, solvemap,
404 stats, msglvl, msgFile, firsttag, MPI_COMM_WORLD) ;
405 if ( msglvl > 2 ) {
406 fprintf(msgFile, "\n\n numeric factorization after split") ;
407 FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
408 fflush(msgFile) ;
409 }
410 /*--------------------------------------------------------------------*/
411 /*
412 ------------------------------------------------
413 STEP 13: permute and redistribute Y if necessary
414 ------------------------------------------------
415 */
416 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
417 IV *rowmapIV ;
418 /*
419 ----------------------------------------------------------
420 pivoting has taken place, redistribute the right hand side
421 to match the final rows and columns in the fronts
422 ----------------------------------------------------------
423 */
424 rowmapIV = FrontMtx_MPI_rowmapIV(frontmtx, ownersIV, msglvl,
425 msgFile, MPI_COMM_WORLD) ;
426 newY = DenseMtx_MPI_splitByRows(mtxY, rowmapIV, stats, msglvl,
427 msgFile, firsttag, MPI_COMM_WORLD) ;
428 DenseMtx_free(mtxY) ;
429 mtxY = newY ;
430 IV_free(rowmapIV) ;
431 }
432 if ( msglvl > 2 ) {
433 fprintf(msgFile, "\n\n rhs matrix after split") ;
434 DenseMtx_writeForHumanEye(mtxY, msgFile) ;
435 fflush(msgFile) ;
436 }
437 /*--------------------------------------------------------------------*/
438 /*
439 ------------------------------------------
440 STEP 14: create a solution DenseMtx object
441 ------------------------------------------
442 */
443 ownedColumnsIV = FrontMtx_ownedColumnsIV(frontmtx, myid, ownersIV,
444 msglvl, msgFile) ;
445 nmycol = IV_size(ownedColumnsIV) ;
446 mtxX = DenseMtx_new() ;
447 if ( nmycol > 0 ) {
448 DenseMtx_init(mtxX, type, 0, 0, nmycol, nrhs, 1, nmycol) ;
449 DenseMtx_rowIndices(mtxX, &nrow, &rowind) ;
450 IVcopy(nmycol, rowind, IV_entries(ownedColumnsIV)) ;
451 }
452 /*--------------------------------------------------------------------*/
453 /*
454 --------------------------------
455 STEP 15: solve the linear system
456 --------------------------------
457 */
458 solvemanager = SubMtxManager_new() ;
459 SubMtxManager_init(solvemanager, NO_LOCK, 0) ;
460 FrontMtx_MPI_solve(frontmtx, mtxX, mtxY, solvemanager, solvemap, cpus,
461 stats, msglvl, msgFile, firsttag, MPI_COMM_WORLD) ;
462 SubMtxManager_free(solvemanager) ;
463 if ( msglvl > 2 ) {
464 fprintf(msgFile, "\n solution in new ordering") ;
465 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
466 }
467 /*--------------------------------------------------------------------*/
468 /*
469 --------------------------------------------------------
470 STEP 15: permute the solution into the original ordering
471 and assemble the solution onto processor zero
472 --------------------------------------------------------
473 */
474 DenseMtx_permuteRows(mtxX, newToOldIV) ;
475 if ( msglvl > 2 ) {
476 fprintf(msgFile, "\n\n solution in old ordering") ;
477 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
478 fflush(msgFile) ;
479 }
480 IV_fill(vtxmapIV, 0) ;
481 firsttag++ ;
482 mtxX = DenseMtx_MPI_splitByRows(mtxX, vtxmapIV, stats, msglvl, msgFile,
483 firsttag, MPI_COMM_WORLD) ;
484 if ( myid == 0 && msglvl > 0 ) {
485 fprintf(msgFile, "\n\n complete solution in old ordering") ;
486 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
487 fflush(msgFile) ;
488 }
489 /*--------------------------------------------------------------------*/
490 MPI_Finalize() ;
491
492 return(1) ; }
493 /*--------------------------------------------------------------------*/
494