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