1 /* split.c */
2
3 #include "../spoolesMPI.h"
4
5 #define MYDEBUG 0
6
7 /*--------------------------------------------------------------------*/
8 /*
9 -----------------------------------------------------------------
10 purpose -- to split a DenseMtx object by rows
11
12 mtx -- DenseMtx object
13 rowmapIV -- map from rows to owning processes
14 firsttag -- first tag to be used in these messages
15 stats[4] -- statistics vector
16 stats[0] -- # of messages sent
17 stats[1] -- # of messages received
18 stats[2] -- # of bytes sent
19 stats[3] -- # of bytes received
20 msglvl -- message level
21 msgFile -- message file
22 comm -- MPI communicator
23
24 return value -- a new DenseMtx object filled with the owned rows
25
26 created -- 98may16, cca
27 modified -- 98sep26, cca
28 mtx is not modified
29 -----------------------------------------------------------------
30 */
31 DenseMtx *
DenseMtx_MPI_splitByRows(DenseMtx * mtx,IV * rowmapIV,int stats[],int msglvl,FILE * msgFile,int firsttag,MPI_Comm comm)32 DenseMtx_MPI_splitByRows (
33 DenseMtx *mtx,
34 IV *rowmapIV,
35 int stats[],
36 int msglvl,
37 FILE *msgFile,
38 int firsttag,
39 MPI_Comm comm
40 ) {
41 DenseMtx *inmtx, *keepmtx, *outmtx ;
42 double *inbuffer, *outbuffer ;
43 int destination, ii, inbuffersize, incount, iproc, irow,
44 lasttag, left, myid, ncol, ndouble, neqns, nkeep,
45 nmoved, nowned, nproc, nrecv, nrow, nsend, tagbound,
46 offset, outbuffersize, outcount, right, source, tag, type ;
47 int *head, *link, *rowind, *rowmap, *rowsToRecv, *rowsToSend ;
48 MPI_Status status ;
49 /*
50 -------------------------------------------------
51 get id of self, # of processes and # of equations
52 -------------------------------------------------
53 */
54 MPI_Comm_rank(comm, &myid) ;
55 MPI_Comm_size(comm, &nproc) ;
56 /*--------------------------------------------------------------------*/
57 {
58 int rc = 1 ;
59 int *rcs = IVinit(nproc, -1) ;
60 /*
61 --------------
62 check the data
63 --------------
64 */
65 if ( mtx == NULL ) {
66 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
67 "\n mtx is NULL\n") ;
68 rc = -1 ;
69 }
70 if ( rowmapIV == NULL ) {
71 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
72 "\n rowmapIV is NULL\n") ;
73 rc = -2 ;
74 }
75 if ( msglvl > 0 && msgFile == NULL ) {
76 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
77 "\n msglvl > 0 and msgFile is NULL\n") ;
78 rc = -3 ;
79 }
80 if ( firsttag < 0 ) {
81 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
82 "\n firsttag = %d\n", firsttag) ;
83 rc = -4 ;
84 }
85 lasttag = firsttag + nproc ;
86 if ( lasttag > (tagbound = maxTagMPI(comm)) ) {
87 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
88 "\n lasttag = %d, tag_bound = %d", lasttag, tagbound) ;
89 rc = -5 ;
90 }
91 MPI_Allgather((void *) &rc, 1, MPI_INT,
92 (void *) rcs, 1, MPI_INT, comm) ;
93 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
94 if ( rcs[iproc] != 1 ) {
95 if ( msgFile != NULL ) {
96 fprintf(msgFile,
97 "\n fatal error in DenseMtx_MPI_splitByRows()"
98 "\n trouble with return code") ;
99 IVfprintf(msgFile, nproc, rcs) ;
100 MPI_Finalize() ;
101 exit(rc) ;
102 }
103 }
104 }
105 IVfree(rcs) ;
106 }
107 /*--------------------------------------------------------------------*/
108 /*
109 -----------------------
110 get type and dimensions
111 -----------------------
112 */
113 type = mtx->type ;
114 IV_sizeAndEntries(rowmapIV, &neqns, &rowmap) ;
115 DenseMtx_dimensions(mtx, &nrow, &ncol) ;
116 if ( msglvl > 2 ) {
117 fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitByRows"
118 "\n nproc = %d, myid = %d, neqns = %d, nrow = %d, ncol = %d",
119 nproc, myid, neqns, nrow, ncol) ;
120 fflush(msgFile) ;
121 }
122 /*
123 -------------------------------------------------------
124 communicate the type's and ncol's to all the processors
125 -------------------------------------------------------
126 */
127 {
128 int *ivec = IVinit(nproc, -1) ;
129 MPI_Allgather((void *) &type, 1, MPI_INT,
130 (void *) ivec, 1, MPI_INT, comm) ;
131 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
132 if ( ivec[iproc] != type ) {
133 if ( msgFile != NULL ) {
134 fprintf(msgFile,
135 "\n fatal error in DenseMtx_MPI_splitByRows()"
136 "\n trouble with types") ;
137 IVfprintf(msgFile, nproc, ivec) ;
138 MPI_Finalize() ;
139 exit(-1) ;
140 }
141 }
142 }
143 MPI_Allgather((void *) &ncol, 1, MPI_INT,
144 (void *) ivec, 1, MPI_INT, comm) ;
145 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
146 if ( ivec[iproc] != ncol ) {
147 if ( msgFile != NULL ) {
148 fprintf(msgFile,
149 "\n fatal error in DenseMtx_MPI_splitByRows()"
150 "\n trouble with ncols") ;
151 IVfprintf(msgFile, nproc, ivec) ;
152 MPI_Finalize() ;
153 exit(-1) ;
154 }
155 }
156 }
157 IVfree(ivec) ;
158 }
159 /*--------------------------------------------------------------------*/
160 /*
161 -------------------------------------------------------------
162 get the counts of the entries to send to the other processors
163 -------------------------------------------------------------
164 */
165 DenseMtx_rowIndices(mtx, &nrow, &rowind) ;
166 rowsToSend = IVinit(2*nproc, 0) ;
167 rowsToRecv = rowsToSend + nproc ;
168 head = IVinit(nproc, -1) ;
169 link = IVinit(nrow, -1) ;
170 for ( ii = 0, nkeep = 0 ; ii < nrow ; ii++ ) {
171 irow = rowind[ii] ;
172 iproc = rowmap[irow] ;
173 link[ii] = head[iproc] ;
174 head[iproc] = ii ;
175 if ( iproc != myid ) {
176 rowsToSend[iproc]++ ;
177 } else {
178 nkeep++ ;
179 }
180 }
181 if ( msglvl > 2 ) {
182 fprintf(msgFile, "\n nkeep = %d, row send counts ", nkeep) ;
183 IVfprintf(msgFile, nproc, rowsToSend) ;
184 fflush(msgFile) ;
185 }
186 /*
187 -------------------------------
188 do an all-to-all gather/scatter
189 -------------------------------
190 */
191 MPI_Alltoall((void *) rowsToSend, 1, MPI_INT,
192 (void *) rowsToRecv, 1, MPI_INT, comm) ;
193 nowned = nkeep + IVsum(nproc, rowsToRecv) ;
194 if ( msglvl > 2 ) {
195 fprintf(msgFile, "\n nkeep = %d, row receive counts ", nkeep) ;
196 IVfprintf(msgFile, nproc, rowsToRecv) ;
197 fflush(msgFile) ;
198 }
199 /*
200 -------------------------
201 determine the buffer size
202 -------------------------
203 */
204 nsend = IVmax(nproc, rowsToSend, &iproc) ;
205 nrecv = IVmax(nproc, rowsToRecv, &iproc) ;
206 if ( msglvl > 2 ) {
207 fprintf(msgFile, "\n nsend = %d, nrecv = %d", nsend, nrecv) ;
208 fflush(msgFile) ;
209 }
210 /*
211 -------------------------------------------
212 allocate the send/receive DenseMtx objects
213 -------------------------------------------
214 */
215 if ( nsend > 0 ) {
216 outmtx = DenseMtx_new() ;
217 if ( mtx->inc1 == 1 ) {
218 DenseMtx_init(outmtx, type, myid, -1, nsend, ncol, 1, nsend) ;
219 } else {
220 DenseMtx_init(outmtx, type, myid, -1, nsend, ncol, ncol, 1) ;
221 }
222 } else {
223 outmtx = NULL ;
224 }
225 if ( nrecv > 0 ) {
226 inmtx = DenseMtx_new() ;
227 if ( mtx->inc1 == 1 ) {
228 DenseMtx_init(inmtx, type, myid, -1, nrecv, ncol, 1, nrecv) ;
229 } else {
230 DenseMtx_init(inmtx, type, myid, -1, nrecv, ncol, ncol, 1) ;
231 }
232 } else {
233 inmtx = NULL ;
234 }
235 /*
236 -------------------------------------
237 allocate the DenseMtx object to keep
238 -------------------------------------
239 */
240 keepmtx = DenseMtx_new() ;
241 if ( mtx->inc1 == 1 ) {
242 DenseMtx_init(keepmtx, type, myid, -1, nowned, ncol, 1, nowned) ;
243 } else {
244 DenseMtx_init(keepmtx, type, myid, -1, nowned, ncol, ncol, 1) ;
245 }
246 if ( msglvl > 2 ) {
247 fprintf(msgFile, "\n keepmtx object allocated") ;
248 fflush(msgFile) ;
249 }
250 /*
251 ----------------------------------------------------------------
252 copy all rows to keep from the input matrix into the keep matrix
253 ----------------------------------------------------------------
254 */
255 for ( ii = head[myid], nmoved = 0 ; ii != -1 ; ii = link[ii] ) {
256 DenseMtx_copyRowAndIndex(keepmtx, nmoved, mtx, ii) ;
257 nmoved++ ;
258 }
259 if ( nmoved > 0 ) {
260 /*
261 if ( msglvl > 3 ) {
262 fprintf(msgFile, "\n\n keepmtx") ;
263 DenseMtx_writeForHumanEye(keepmtx, msgFile) ;
264 fflush(msgFile) ;
265 }
266 */
267 }
268 /*
269 --------------------------------------------------------------
270 loop over the processes, gather their values and send them off
271 --------------------------------------------------------------
272 */
273 for ( offset = 1, tag = firsttag ; offset < nproc ; offset++, tag++ ) {
274 right = (myid + offset) % nproc ;
275 left = (nproc + myid - offset) % nproc ;
276 outcount = rowsToSend[right] ;
277 incount = rowsToRecv[left] ;
278 if ( msglvl > 2 ) {
279 fprintf(msgFile,
280 "\n\n ### process %d, send %d to right %d, recv %d from left %d",
281 myid, outcount, right, incount, left) ;
282 fflush(msgFile) ;
283 }
284 if ( outcount > 0 ) {
285 /*
286 --------------------------
287 load the out matrix object
288 --------------------------
289 */
290 if ( mtx->inc1 == 1 ) {
291 DenseMtx_init(outmtx, type,
292 myid, -1, outcount, ncol, 1, outcount) ;
293 } else {
294 DenseMtx_init(outmtx, type,
295 myid, -1, outcount, ncol, ncol, 1) ;
296 }
297 for ( ii = head[right], nmoved = 0 ; ii != -1 ; ii = link[ii] ) {
298 DenseMtx_copyRowAndIndex(outmtx, nmoved, mtx, ii) ;
299 nmoved++ ;
300 }
301 destination = right ;
302 if ( msglvl > 3 ) {
303 fprintf(msgFile, "\n\n outmtx for process %d", destination) ;
304 DenseMtx_writeForHumanEye(outmtx, msgFile) ;
305 fflush(msgFile) ;
306 }
307 stats[0]++ ;
308 stats[2] += sizeof(double)*DV_size(&outmtx->wrkDV) ;
309 } else {
310 /*
311 ------------------------------------------
312 set the destination to be the NULL process
313 ------------------------------------------
314 */
315 destination = MPI_PROC_NULL ;
316 }
317 if ( incount > 0 ) {
318 /*
319 ----------------------------------
320 initialize the input matrix object
321 ----------------------------------
322 */
323 if ( mtx->inc1 == 1 ) {
324 DenseMtx_init(inmtx, type,
325 myid, -1, incount, ncol, 1, incount) ;
326 } else {
327 DenseMtx_init(inmtx, type,
328 myid, -1, incount, ncol, ncol, 1) ;
329 }
330 if ( msglvl > 2 ) {
331 fprintf(msgFile, "\n\n inmtx initialized to have %d rows",
332 incount) ;
333 fflush(msgFile) ;
334 }
335 source = left ;
336 stats[1]++ ;
337 stats[3] += sizeof(double)*DV_size(&inmtx->wrkDV) ;
338 } else {
339 source = MPI_PROC_NULL ;
340 }
341 /*
342 -----------------
343 do a send/receive
344 -----------------
345 */
346 inbuffersize = outbuffersize = 0 ;
347 inbuffer = outbuffer = NULL ;
348 if ( outmtx != NULL ) {
349 outbuffersize = DV_size(&outmtx->wrkDV) ;
350 outbuffer = DV_entries(&outmtx->wrkDV) ;
351 }
352 if ( inmtx != NULL ) {
353 inbuffersize = DV_size(&inmtx->wrkDV) ;
354 inbuffer = DV_entries(&inmtx->wrkDV) ;
355 }
356 if ( msglvl > 2 ) {
357 fprintf(msgFile,
358 "\n inbuffersize = %d, inbuffer = %p"
359 "\n outbuffersize = %d, outbuffer = %p",
360 inbuffersize, inbuffer, outbuffersize, outbuffer) ;
361 fflush(msgFile) ;
362 }
363 MPI_Sendrecv((void*) outbuffer, outbuffersize, MPI_DOUBLE,
364 destination, tag, (void*) inbuffer, inbuffersize,
365 MPI_DOUBLE, source, tag, comm, &status) ;
366 if ( msglvl > 2 ) {
367 MPI_Get_count(&status, MPI_DOUBLE, &ndouble) ;
368 fprintf(msgFile,
369 "\n\n message received, source %d, tag %d, double count %d",
370 status.MPI_SOURCE, status.MPI_TAG, ndouble) ;
371 fprintf(msgFile, "\n MPI_ERROR = %d", status.MPI_ERROR) ;
372 fflush(msgFile) ;
373 }
374 if ( source != MPI_PROC_NULL ) {
375 /*
376 -------------------------------------
377 initialize the object from its buffer
378 -------------------------------------
379 */
380 DenseMtx_initFromBuffer(inmtx) ;
381 if ( msglvl > 2 ) {
382 fprintf(msgFile,
383 "\n DenseMtx object intialized from its buffer") ;
384 fflush(msgFile) ;
385 }
386 if ( msglvl > 4 ) {
387 DenseMtx_writeForHumanEye(inmtx, msgFile) ;
388 fflush(msgFile) ;
389 }
390 }
391 if ( incount > 0 ) {
392 if ( nkeep + incount > nowned ) {
393 fprintf(msgFile, "\n fatal error in DenseMtx_splitByRows()"
394 "\n nkeep = %d, nrecv = %d, nowned = %d",
395 nkeep, nrecv, nowned) ;
396 exit(-1) ;
397 }
398 for ( irow = 0 ; irow < incount ; irow++, nkeep++ ) {
399 DenseMtx_copyRowAndIndex(keepmtx, nkeep, inmtx, irow) ;
400 }
401 }
402 }
403 /*
404 -------------------------
405 sort the rows and columns
406 -------------------------
407 */
408 DenseMtx_sort(keepmtx) ;
409 /*
410 ------------------------------------------------------
411 check that the matrix contains only the rows it should
412 ------------------------------------------------------
413 */
414 nrow = keepmtx->nrow ;
415 rowind = keepmtx->rowind ;
416 for ( ii = 0 ; ii < nrow ; ii++ ) {
417 irow = rowind[ii] ;
418 if ( irow < 0 || irow >= neqns ) {
419 fprintf(stderr,
420 "\n process %d : local row %d, global row %d, neqns = %d\n",
421 myid, ii, irow, neqns) ;
422 exit(-1) ;
423 }
424 if ( rowmap[irow] != myid ) {
425 fprintf(stderr,
426 "\n process %d : local row %d, global row %d, map = %d\n",
427 myid, ii, irow, rowmap[irow]) ;
428 exit(-1) ;
429 }
430 }
431 /*
432 ------------------------
433 free the working storage
434 ------------------------
435 */
436 if ( outmtx != NULL ) {
437 DenseMtx_free(outmtx) ;
438 }
439 if ( inmtx != NULL ) {
440 DenseMtx_free(inmtx) ;
441 }
442 IVfree(rowsToSend) ;
443 IVfree(head) ;
444 IVfree(link) ;
445
446 return(keepmtx) ; }
447
448 /*--------------------------------------------------------------------*/
449 /*
450 -----------------------------------------------------------------
451 purpose -- to scatter a DenseMtx object by rows
452
453 Xglobal -- global DenseMtx object, significant only for root
454 Xlocal -- local DenseMtx object, if NULL on input, will
455 be created if necessary
456 rowmapIV -- map from rows to owning processes
457 firsttag -- first tag to be used in these messages
458 stats[4] -- statistics vector
459 stats[0] -- # of messages sent
460 stats[1] -- # of messages received
461 stats[2] -- # of bytes sent
462 stats[3] -- # of bytes received
463 msglvl -- message level
464 msgFile -- message file
465 comm -- MPI communicator
466
467 return value -- Xlocal, a local DenseMtx object
468
469 created -- 98may16, cca
470 modified -- 98sep26, cca
471 mtx is not modified
472 -----------------------------------------------------------------
473 */
474 DenseMtx *
DenseMtx_MPI_splitFromGlobalByRows(DenseMtx * Xglobal,DenseMtx * Xlocal,IV * rowmapIV,int root,int stats[],int msglvl,FILE * msgFile,int firsttag,MPI_Comm comm)475 DenseMtx_MPI_splitFromGlobalByRows (
476 DenseMtx *Xglobal,
477 DenseMtx *Xlocal,
478 IV *rowmapIV,
479 int root,
480 int stats[],
481 int msglvl,
482 FILE *msgFile,
483 int firsttag,
484 MPI_Comm comm
485 ) {
486 DenseMtx *tempmtx ;
487 double *buffer ;
488 int ii, iproc, irow, maxnrow, myid, ncolX, nkeep,
489 nproc, nrowloc, nrowmap, nrowX, nsend, size, type ;
490 int *counts, *head, *link, *rowind, *rowmap ;
491 MPI_Status status ;
492 /*
493 -------------------------------------------------
494 get id of self, # of processes and # of equations
495 -------------------------------------------------
496 */
497 MPI_Comm_rank(comm, &myid) ;
498 MPI_Comm_size(comm, &nproc) ;
499 if ( root < 0 || root >= nproc ) {
500 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
501 "\n root = %d, nproc = %d\n", root, nproc) ;
502 MPI_Finalize() ;
503 exit(-1) ;
504 }
505 /*--------------------------------------------------------------------*/
506 /*
507 --------------
508 check the data
509 --------------
510 */
511 {
512 int rc = 1 ;
513 int *rcs = IVinit(nproc, -1) ;
514
515 if ( myid == root ) {
516 if ( Xglobal == NULL ) {
517 fprintf(stderr,
518 "\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
519 "\n Xglobal is NULL\n") ;
520 rc = -1 ;
521 }
522 if ( rowmapIV == NULL ) {
523 fprintf(stderr,
524 "\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
525 "\n rowmapIV is NULL\n") ;
526 rc = -2 ;
527 }
528 }
529 if ( msglvl > 0 && msgFile == NULL ) {
530 fprintf(msgFile,
531 "\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
532 "\n msglvl > 0 and msgFile = NULL\n") ;
533 rc = -3 ;
534 }
535 if ( firsttag < 0 ) {
536 fprintf(stderr,
537 "\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
538 "\n firsttag = %d\n", firsttag) ;
539 rc = -4 ;
540 }
541 MPI_Allgather((void *) &rc, 1, MPI_INT,
542 (void *) rcs, 1, MPI_INT, comm) ;
543 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
544 if ( rcs[iproc] != 1 ) {
545 if ( msgFile != NULL ) {
546 fprintf(msgFile,
547 "\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
548 "\n trouble with return code") ;
549 IVfprintf(msgFile, nproc, rcs) ;
550 MPI_Finalize() ;
551 exit(rc) ;
552 }
553 }
554 }
555 IVfree(rcs) ;
556 }
557 /*--------------------------------------------------------------------*/
558
559 if ( myid == root ) {
560 type = Xglobal->type ;
561 IV_sizeAndEntries(rowmapIV, &nrowmap, &rowmap) ;
562 DenseMtx_dimensions(Xglobal, &nrowX, &ncolX) ;
563 if ( msglvl > 2 ) {
564 fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitFromGlobalByRows"
565 "\n nproc = %d, myid = %d, nrowmap = %d, nrowX = %d, ncolX = %d",
566 nproc, myid, nrowmap, nrowX, ncolX) ;
567 fflush(msgFile) ;
568 }
569 }
570 /*
571 ----------------------------------------
572 broadcast the type of entries and number
573 of right hand sides to all processors
574 ----------------------------------------
575 */
576 MPI_Bcast((void *) &type, 1, MPI_INT, root, comm) ;
577 MPI_Bcast((void *) &ncolX, 1, MPI_INT, root, comm) ;
578 stats[0] += 2 ;
579 stats[1] += 2 ;
580 stats[2] += 2*sizeof(int) ;
581 stats[3] += 2*sizeof(int) ;
582 if ( msglvl > 2 ) {
583 fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitFromGlobalByRows"
584 "\n type %d, ncolX %d", type, ncolX) ;
585 fflush(msgFile) ;
586 }
587 if ( myid == root ) {
588 /*
589 ------------------------------------------------
590 create a head/link structure for the matrix rows
591 ------------------------------------------------
592 */
593 DenseMtx_rowIndices(Xglobal, &nrowX, &rowind) ;
594 counts = IVinit(nproc, 0) ;
595 head = IVinit(nproc, -1) ;
596 link = IVinit(nrowX, -1) ;
597 for ( ii = nrowX - 1 ; ii >= 0 ; ii-- ) {
598 irow = rowind[ii] ;
599 iproc = rowmap[irow] ;
600 link[ii] = head[iproc] ;
601 head[iproc] = ii ;
602 counts[iproc]++ ;
603 }
604 } else {
605 counts = NULL ;
606 }
607 /*
608 -------------------------------------------------
609 communicate the number of rows for each processor
610 -------------------------------------------------
611 */
612 MPI_Scatter((void *) counts, 1, MPI_INT,
613 (void *) &nrowloc, 1, MPI_INT, root, comm) ;
614 stats[0] += 1 ;
615 stats[1] += 1 ;
616 stats[2] += (nproc-1)*sizeof(int) ;
617 stats[3] += (nproc-1)*sizeof(int) ;
618 if ( msglvl > 2 ) {
619 fprintf(msgFile, "\n\n nrowloc = %d", nrowloc) ;
620 fflush(msgFile) ;
621 }
622 /*
623 -----------------------------------------------------------
624 if necessary, create the local Xloc matrix, then initialize
625 -----------------------------------------------------------
626 */
627 if ( Xlocal == NULL ) {
628 Xlocal = DenseMtx_new() ;
629 }
630 DenseMtx_init(Xlocal, type, -1, -1, nrowloc, ncolX, 1, nrowloc) ;
631 if ( msglvl > 2 ) {
632 fprintf(msgFile, "\n\n Xlocal after initialization") ;
633 DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
634 fflush(msgFile) ;
635 }
636 if ( myid == root ) {
637 /*
638 ---------------------------------
639 load local matrix with owned rows
640 ---------------------------------
641 */
642 if ( nrowloc > 0 ) {
643 int iglob, iloc = 0 ;
644 for ( iglob = head[root] ; iglob != -1 ; iglob = link[iglob] ) {
645 DenseMtx_copyRowAndIndex(Xlocal, iloc, Xglobal, iglob) ;
646 iloc++ ;
647 }
648 if ( msglvl > 2 ) {
649 fprintf(msgFile, "\n\n Xlocal on root") ;
650 DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
651 fflush(msgFile) ;
652 }
653 }
654 /*
655 -----------------------------------------------------
656 create a temporary matrix to send to other processors
657 -----------------------------------------------------
658 */
659 counts[myid] = 0 ;
660 maxnrow = IVmax(nproc, counts, &iproc) ;
661 if ( maxnrow > 0 ) {
662 DenseMtx *tempmtx = DenseMtx_new() ;
663 DenseMtx_init(tempmtx, type, -1, -1, maxnrow, ncolX, 1, maxnrow) ;
664 /*
665 ----------------------------------
666 loop over the processors
667 collect owned rows into tempmtx
668 send tempmtx to processor
669 ----------------------------------
670 */
671 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
672 if ( iproc != root && (nrowloc = counts[iproc]) > 0 ) {
673 DenseMtx_init(tempmtx, type, -1, -1,
674 nrowloc, ncolX, 1, nrowloc) ;
675 nsend = 0 ;
676 for ( ii = head[iproc] ; ii != -1 ; ii = link[ii] ) {
677 DenseMtx_copyRowAndIndex(tempmtx, nsend, Xglobal, ii) ;
678 nsend++ ;
679 }
680 if ( msglvl > 2 ) {
681 fprintf(msgFile, "\n\n tempmtx for proc %d", iproc) ;
682 DenseMtx_writeForHumanEye(tempmtx, msgFile) ;
683 fflush(msgFile) ;
684 }
685 size = DV_size(&tempmtx->wrkDV) ;
686 buffer = DV_entries(&tempmtx->wrkDV) ;
687 MPI_Send((void *) buffer, size, MPI_DOUBLE,
688 iproc, firsttag, comm) ;
689 stats[0] += 1 ;
690 stats[2] += size*sizeof(double) ;
691 }
692 }
693 DenseMtx_free(tempmtx) ;
694 }
695 /*
696 ------------------------
697 free the working storage
698 ------------------------
699 */
700 IVfree(head) ;
701 IVfree(link) ;
702 IVfree(counts) ;
703 } else {
704 /*
705 ------------------
706 non-root processor
707 ------------------
708 */
709 if ( nrowloc > 0 ) {
710 size = DV_size(&Xlocal->wrkDV) ;
711 buffer = DV_entries(&Xlocal->wrkDV) ;
712 if ( msglvl > 2 ) {
713 fprintf(msgFile, "\n\n size = %d, buffer = %p", size, buffer) ;
714 fflush(msgFile) ;
715 }
716 MPI_Recv((void *) buffer, size, MPI_DOUBLE,
717 root, firsttag, comm, &status);
718 stats[1] += 1 ;
719 stats[3] += size*sizeof(double) ;
720 if ( msglvl > 2 ) {
721 fprintf(msgFile, "\n\n Xlocal rec'd from root %d", root) ;
722 fflush(msgFile) ;
723 }
724 DenseMtx_initFromBuffer(Xlocal) ;
725 if ( msglvl > 2 ) {
726 DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
727 fflush(msgFile) ;
728 }
729 } else {
730 Xlocal = NULL ;
731 }
732 }
733 if ( msglvl > 3 ) {
734 if ( Xlocal != NULL ) {
735 fprintf(msgFile, "\n\n Xlocal") ;
736 DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
737 } else {
738 fprintf(msgFile, "\n\n Xlocal is NULL") ;
739 }
740 fflush(msgFile) ;
741 }
742 if ( msglvl > 2 ) {
743 fprintf(msgFile, "\n\n leaving DenseMtx_splitFromGlobalByRows()") ;
744 fflush(msgFile) ;
745 }
746 return(Xlocal) ; }
747
748 /*--------------------------------------------------------------------*/
749 /*
750 -----------------------------------------------------------------
751 purpose -- to merge a DenseMtx object by rows
752
753 Xlocal -- DenseMtx object, can be NULL
754 Xglobal -- DenseMtx object, can be NULL
755 significant only for root
756 firsttag -- first tag to be used in these messages
757 stats[4] -- statistics vector
758 stats[0] -- # of messages sent
759 stats[1] -- # of messages received
760 stats[2] -- # of bytes sent
761 stats[3] -- # of bytes received
762 msglvl -- message level
763 msgFile -- message file
764 comm -- MPI communicator
765
766 return value --
767 if processor is root
768 Xglobal is returned, if was NULL on input, it is created
769 else
770 NULL
771 endif
772
773 Xlocal is not modified
774
775 created -- 98sep27, cca
776 -----------------------------------------------------------------
777 */
778 DenseMtx *
DenseMtx_MPI_mergeToGlobalByRows(DenseMtx * Xglobal,DenseMtx * Xlocal,int root,int stats[],int msglvl,FILE * msgFile,int firsttag,MPI_Comm comm)779 DenseMtx_MPI_mergeToGlobalByRows (
780 DenseMtx *Xglobal,
781 DenseMtx *Xlocal,
782 int root,
783 int stats[],
784 int msglvl,
785 FILE *msgFile,
786 int firsttag,
787 MPI_Comm comm
788 ) {
789 double *buffer ;
790 int iproc, maxnrow, myid, ncolX, nproc, nrowXloc, size, type ;
791 int *incounts ;
792 /*
793 -------------------------------------------------
794 get id of self, # of processes and # of equations
795 -------------------------------------------------
796 */
797 MPI_Comm_rank(comm, &myid) ;
798 MPI_Comm_size(comm, &nproc) ;
799 if ( root < 0 || root >= nproc ) {
800 fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
801 "\n root = %d, nproc = %d\n", root, nproc) ;
802 MPI_Finalize() ;
803 exit(-1) ;
804 }
805 /*--------------------------------------------------------------------*/
806 /*
807 --------------
808 check the data
809 --------------
810 */
811 {
812 int rc = 1 ;
813 int *ivec = IVinit(nproc, -1) ;
814
815 if ( msglvl > 0 && msgFile == NULL ) {
816 fprintf(msgFile,
817 "\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
818 "\n msglvl > 0 and msgFile = NULL\n") ;
819 rc = -2 ;
820 }
821 if ( firsttag < 0 ) {
822 fprintf(stderr,
823 "\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
824 "\n firsttag = %d\n", firsttag) ;
825 rc = -3 ;
826 }
827 MPI_Allgather((void *) &rc, 1, MPI_INT,
828 (void *) ivec, 1, MPI_INT, comm) ;
829 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
830 if ( ivec[iproc] != 1 ) {
831 if ( msgFile != NULL ) {
832 fprintf(msgFile,
833 "\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
834 "\n trouble with return code") ;
835 IVfprintf(msgFile, nproc, ivec) ;
836 MPI_Finalize() ;
837 exit(rc) ;
838 }
839 }
840 }
841 /*
842 ------------------------------
843 check the type of the matrices
844 ------------------------------
845 */
846 type = (Xlocal != NULL) ? Xlocal->type : -1 ;
847 MPI_Allgather((void *) &type, 1, MPI_INT,
848 (void *) ivec, 1, MPI_INT, comm) ;
849 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
850 if ( ivec[iproc] != -1 ) {
851 if ( type == -1 ) {
852 type = ivec[iproc] ;
853 } else if ( type != ivec[iproc] ) {
854 if ( msgFile != NULL ) {
855 fprintf(msgFile,
856 "\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
857 "\n trouble with types\n") ;
858 IVfprintf(msgFile, nproc, ivec) ;
859 MPI_Finalize() ;
860 exit(-1) ;
861 }
862 }
863 }
864 }
865 /*
866 --------------------------------------
867 check the # of columns of the matrices
868 --------------------------------------
869 */
870 ncolX = (Xlocal != NULL) ? Xlocal->ncol : 0 ;
871 MPI_Allgather((void *) &ncolX, 1, MPI_INT,
872 (void *) ivec, 1, MPI_INT, comm) ;
873 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
874 if ( ivec[iproc] != 0 ) {
875 if ( ncolX == 0 ) {
876 ncolX = ivec[iproc] ;
877 } else if ( ncolX != ivec[iproc] ) {
878 if ( msgFile != NULL ) {
879 fprintf(msgFile,
880 "\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
881 "\n trouble with ncolX\n") ;
882 IVfprintf(msgFile, nproc, ivec) ;
883 MPI_Finalize() ;
884 exit(-1) ;
885 }
886 }
887 }
888 }
889 IVfree(ivec) ;
890 }
891 /*--------------------------------------------------------------------*/
892 /*
893 ------------------------------------------------------
894 gather the number of incoming rows onto processor root
895 ------------------------------------------------------
896 */
897 nrowXloc = (Xlocal != NULL) ? Xlocal->nrow : 0 ;
898 if ( myid == root ) {
899 incounts = IVinit(nproc, 0) ;
900 } else {
901 incounts = NULL ;
902 }
903 MPI_Gather(&nrowXloc, 1, MPI_INT,
904 (void *) incounts, 1, MPI_INT, root, comm) ;
905 if ( myid == root ) {
906 DenseMtx *tempmtx ;
907 int count, iglob, iloc, incount, nrowXglobal ;
908 MPI_Status status ;
909 /*
910 ----------------------------------------------------------
911 if necessary, create the global matrix, then initialize it
912 ----------------------------------------------------------
913 */
914 nrowXglobal = IVsum(nproc, incounts) ;
915 if ( Xglobal == NULL ) {
916 Xglobal = DenseMtx_new() ;
917 }
918 DenseMtx_init(Xglobal, type, -1, -1,
919 nrowXglobal, ncolX, 1, nrowXglobal) ;
920 /*
921 ---------------------------------
922 load local matrix with owned rows
923 ---------------------------------
924 */
925 for ( iloc = iglob = 0 ; iloc < nrowXloc ; iloc++, iglob++ ) {
926 DenseMtx_copyRowAndIndex(Xglobal, iglob, Xlocal, iloc) ;
927 }
928 if ( msglvl > 2 ) {
929 fprintf(msgFile, "\n\n after loading Xlocal on root") ;
930 DenseMtx_writeForHumanEye(Xglobal, msgFile) ;
931 fflush(msgFile) ;
932 }
933 /*
934 ----------------------------------------------------------
935 create a temporary matrix to receive from other processors
936 ----------------------------------------------------------
937 */
938 incounts[root] = 0 ;
939 maxnrow = IVmax(nproc, incounts, &iproc) ;
940 tempmtx = DenseMtx_new() ;
941 DenseMtx_init(tempmtx, type, -1, -1, maxnrow, ncolX, 1, maxnrow) ;
942 size = DV_size(&tempmtx->wrkDV) ;
943 buffer = DV_entries(&tempmtx->wrkDV) ;
944 /*
945 ----------------------------
946 loop over the processors
947 receive rows into tempmtx
948 ----------------------------
949 */
950 for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
951 if ( iproc != root && (incount = incounts[iproc]) > 0 ) {
952 MPI_Recv((void *) buffer, size, MPI_DOUBLE,
953 iproc, firsttag, comm, &status) ;
954 MPI_Get_count(&status, MPI_DOUBLE, &count) ;
955 stats[1] += 1 ;
956 stats[3] += count*sizeof(double) ;
957 DenseMtx_initFromBuffer(tempmtx) ;
958 for ( iloc = 0 ; iloc < incount ; iloc++, iglob++ ) {
959 DenseMtx_copyRowAndIndex(Xglobal, iglob, tempmtx, iloc) ;
960 }
961 }
962 }
963 /*
964 ------------------------
965 free the working storage
966 ------------------------
967 */
968 IVfree(incounts) ;
969 DenseMtx_free(tempmtx) ;
970 } else {
971 /*
972 ------------------
973 non-root processor
974 ------------------
975 */
976 if ( nrowXloc > 0 ) {
977 size = DV_size(&Xlocal->wrkDV) ;
978 buffer = DV_entries(&Xlocal->wrkDV) ;
979 if ( msglvl > 2 ) {
980 fprintf(msgFile, "\n\n size = %d, buffer = %p", size, buffer) ;
981 fflush(msgFile) ;
982 }
983 MPI_Send((void *) buffer, size, MPI_DOUBLE,
984 root, firsttag, comm) ;
985 stats[0] += 1 ;
986 stats[2] += size*sizeof(double) ;
987 if ( msglvl > 2 ) {
988 fprintf(msgFile, "\n\n Xlocal sent to root %d", root) ;
989 fflush(msgFile) ;
990 }
991 }
992 Xglobal = NULL ;
993 }
994 if ( msglvl > 2 ) {
995 fprintf(msgFile, "\n\n leaving DenseMtx_mergeToGlobalByRows()") ;
996 fflush(msgFile) ;
997 }
998 return(Xglobal) ; }
999
1000 /*--------------------------------------------------------------------*/
1001