/* split.c */
#include "../spoolesMPI.h"
#define MYDEBUG 0
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------------------
purpose -- to split a DenseMtx object by rows
mtx -- DenseMtx object
rowmapIV -- map from rows to owning processes
firsttag -- first tag to be used in these messages
stats[4] -- statistics vector
stats[0] -- # of messages sent
stats[1] -- # of messages received
stats[2] -- # of bytes sent
stats[3] -- # of bytes received
msglvl -- message level
msgFile -- message file
comm -- MPI communicator
return value -- a new DenseMtx object filled with the owned rows
created -- 98may16, cca
modified -- 98sep26, cca
mtx is not modified
-----------------------------------------------------------------
*/
DenseMtx *
DenseMtx_MPI_splitByRows (
DenseMtx *mtx,
IV *rowmapIV,
int stats[],
int msglvl,
FILE *msgFile,
int firsttag,
MPI_Comm comm
) {
DenseMtx *inmtx, *keepmtx, *outmtx ;
double *inbuffer, *outbuffer ;
int destination, ii, inbuffersize, incount, iproc, irow,
lasttag, left, myid, ncol, ndouble, neqns, nkeep,
nmoved, nowned, nproc, nrecv, nrow, nsend, tagbound,
offset, outbuffersize, outcount, right, source, tag, type ;
int *head, *link, *rowind, *rowmap, *rowsToRecv, *rowsToSend ;
MPI_Status status ;
/*
-------------------------------------------------
get id of self, # of processes and # of equations
-------------------------------------------------
*/
MPI_Comm_rank(comm, &myid) ;
MPI_Comm_size(comm, &nproc) ;
/*--------------------------------------------------------------------*/
{
int rc = 1 ;
int *rcs = IVinit(nproc, -1) ;
/*
--------------
check the data
--------------
*/
if ( mtx == NULL ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n mtx is NULL\n") ;
rc = -1 ;
}
if ( rowmapIV == NULL ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n rowmapIV is NULL\n") ;
rc = -2 ;
}
if ( msglvl > 0 && msgFile == NULL ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n msglvl > 0 and msgFile is NULL\n") ;
rc = -3 ;
}
if ( firsttag < 0 ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n firsttag = %d\n", firsttag) ;
rc = -4 ;
}
lasttag = firsttag + nproc ;
if ( lasttag > (tagbound = maxTagMPI(comm)) ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n lasttag = %d, tag_bound = %d", lasttag, tagbound) ;
rc = -5 ;
}
MPI_Allgather((void *) &rc, 1, MPI_INT,
(void *) rcs, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( rcs[iproc] != 1 ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_splitByRows()"
"\n trouble with return code") ;
IVfprintf(msgFile, nproc, rcs) ;
MPI_Finalize() ;
exit(rc) ;
}
}
}
IVfree(rcs) ;
}
/*--------------------------------------------------------------------*/
/*
-----------------------
get type and dimensions
-----------------------
*/
type = mtx->type ;
IV_sizeAndEntries(rowmapIV, &neqns, &rowmap) ;
DenseMtx_dimensions(mtx, &nrow, &ncol) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitByRows"
"\n nproc = %d, myid = %d, neqns = %d, nrow = %d, ncol = %d",
nproc, myid, neqns, nrow, ncol) ;
fflush(msgFile) ;
}
/*
-------------------------------------------------------
communicate the type's and ncol's to all the processors
-------------------------------------------------------
*/
{
int *ivec = IVinit(nproc, -1) ;
MPI_Allgather((void *) &type, 1, MPI_INT,
(void *) ivec, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( ivec[iproc] != type ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_splitByRows()"
"\n trouble with types") ;
IVfprintf(msgFile, nproc, ivec) ;
MPI_Finalize() ;
exit(-1) ;
}
}
}
MPI_Allgather((void *) &ncol, 1, MPI_INT,
(void *) ivec, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( ivec[iproc] != ncol ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_splitByRows()"
"\n trouble with ncols") ;
IVfprintf(msgFile, nproc, ivec) ;
MPI_Finalize() ;
exit(-1) ;
}
}
}
IVfree(ivec) ;
}
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------------
get the counts of the entries to send to the other processors
-------------------------------------------------------------
*/
DenseMtx_rowIndices(mtx, &nrow, &rowind) ;
rowsToSend = IVinit(2*nproc, 0) ;
rowsToRecv = rowsToSend + nproc ;
head = IVinit(nproc, -1) ;
link = IVinit(nrow, -1) ;
for ( ii = 0, nkeep = 0 ; ii < nrow ; ii++ ) {
irow = rowind[ii] ;
iproc = rowmap[irow] ;
link[ii] = head[iproc] ;
head[iproc] = ii ;
if ( iproc != myid ) {
rowsToSend[iproc]++ ;
} else {
nkeep++ ;
}
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n nkeep = %d, row send counts ", nkeep) ;
IVfprintf(msgFile, nproc, rowsToSend) ;
fflush(msgFile) ;
}
/*
-------------------------------
do an all-to-all gather/scatter
-------------------------------
*/
MPI_Alltoall((void *) rowsToSend, 1, MPI_INT,
(void *) rowsToRecv, 1, MPI_INT, comm) ;
nowned = nkeep + IVsum(nproc, rowsToRecv) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n nkeep = %d, row receive counts ", nkeep) ;
IVfprintf(msgFile, nproc, rowsToRecv) ;
fflush(msgFile) ;
}
/*
-------------------------
determine the buffer size
-------------------------
*/
nsend = IVmax(nproc, rowsToSend, &iproc) ;
nrecv = IVmax(nproc, rowsToRecv, &iproc) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n nsend = %d, nrecv = %d", nsend, nrecv) ;
fflush(msgFile) ;
}
/*
-------------------------------------------
allocate the send/receive DenseMtx objects
-------------------------------------------
*/
if ( nsend > 0 ) {
outmtx = DenseMtx_new() ;
if ( mtx->inc1 == 1 ) {
DenseMtx_init(outmtx, type, myid, -1, nsend, ncol, 1, nsend) ;
} else {
DenseMtx_init(outmtx, type, myid, -1, nsend, ncol, ncol, 1) ;
}
} else {
outmtx = NULL ;
}
if ( nrecv > 0 ) {
inmtx = DenseMtx_new() ;
if ( mtx->inc1 == 1 ) {
DenseMtx_init(inmtx, type, myid, -1, nrecv, ncol, 1, nrecv) ;
} else {
DenseMtx_init(inmtx, type, myid, -1, nrecv, ncol, ncol, 1) ;
}
} else {
inmtx = NULL ;
}
/*
-------------------------------------
allocate the DenseMtx object to keep
-------------------------------------
*/
keepmtx = DenseMtx_new() ;
if ( mtx->inc1 == 1 ) {
DenseMtx_init(keepmtx, type, myid, -1, nowned, ncol, 1, nowned) ;
} else {
DenseMtx_init(keepmtx, type, myid, -1, nowned, ncol, ncol, 1) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n keepmtx object allocated") ;
fflush(msgFile) ;
}
/*
----------------------------------------------------------------
copy all rows to keep from the input matrix into the keep matrix
----------------------------------------------------------------
*/
for ( ii = head[myid], nmoved = 0 ; ii != -1 ; ii = link[ii] ) {
DenseMtx_copyRowAndIndex(keepmtx, nmoved, mtx, ii) ;
nmoved++ ;
}
if ( nmoved > 0 ) {
/*
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n keepmtx") ;
DenseMtx_writeForHumanEye(keepmtx, msgFile) ;
fflush(msgFile) ;
}
*/
}
/*
--------------------------------------------------------------
loop over the processes, gather their values and send them off
--------------------------------------------------------------
*/
for ( offset = 1, tag = firsttag ; offset < nproc ; offset++, tag++ ) {
right = (myid + offset) % nproc ;
left = (nproc + myid - offset) % nproc ;
outcount = rowsToSend[right] ;
incount = rowsToRecv[left] ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n\n ### process %d, send %d to right %d, recv %d from left %d",
myid, outcount, right, incount, left) ;
fflush(msgFile) ;
}
if ( outcount > 0 ) {
/*
--------------------------
load the out matrix object
--------------------------
*/
if ( mtx->inc1 == 1 ) {
DenseMtx_init(outmtx, type,
myid, -1, outcount, ncol, 1, outcount) ;
} else {
DenseMtx_init(outmtx, type,
myid, -1, outcount, ncol, ncol, 1) ;
}
for ( ii = head[right], nmoved = 0 ; ii != -1 ; ii = link[ii] ) {
DenseMtx_copyRowAndIndex(outmtx, nmoved, mtx, ii) ;
nmoved++ ;
}
destination = right ;
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n outmtx for process %d", destination) ;
DenseMtx_writeForHumanEye(outmtx, msgFile) ;
fflush(msgFile) ;
}
stats[0]++ ;
stats[2] += sizeof(double)*DV_size(&outmtx->wrkDV) ;
} else {
/*
------------------------------------------
set the destination to be the NULL process
------------------------------------------
*/
destination = MPI_PROC_NULL ;
}
if ( incount > 0 ) {
/*
----------------------------------
initialize the input matrix object
----------------------------------
*/
if ( mtx->inc1 == 1 ) {
DenseMtx_init(inmtx, type,
myid, -1, incount, ncol, 1, incount) ;
} else {
DenseMtx_init(inmtx, type,
myid, -1, incount, ncol, ncol, 1) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n inmtx initialized to have %d rows",
incount) ;
fflush(msgFile) ;
}
source = left ;
stats[1]++ ;
stats[3] += sizeof(double)*DV_size(&inmtx->wrkDV) ;
} else {
source = MPI_PROC_NULL ;
}
/*
-----------------
do a send/receive
-----------------
*/
inbuffersize = outbuffersize = 0 ;
inbuffer = outbuffer = NULL ;
if ( outmtx != NULL ) {
outbuffersize = DV_size(&outmtx->wrkDV) ;
outbuffer = DV_entries(&outmtx->wrkDV) ;
}
if ( inmtx != NULL ) {
inbuffersize = DV_size(&inmtx->wrkDV) ;
inbuffer = DV_entries(&inmtx->wrkDV) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n inbuffersize = %d, inbuffer = %p"
"\n outbuffersize = %d, outbuffer = %p",
inbuffersize, inbuffer, outbuffersize, outbuffer) ;
fflush(msgFile) ;
}
MPI_Sendrecv((void*) outbuffer, outbuffersize, MPI_DOUBLE,
destination, tag, (void*) inbuffer, inbuffersize,
MPI_DOUBLE, source, tag, comm, &status) ;
if ( msglvl > 2 ) {
MPI_Get_count(&status, MPI_DOUBLE, &ndouble) ;
fprintf(msgFile,
"\n\n message received, source %d, tag %d, double count %d",
status.MPI_SOURCE, status.MPI_TAG, ndouble) ;
fprintf(msgFile, "\n MPI_ERROR = %d", status.MPI_ERROR) ;
fflush(msgFile) ;
}
if ( source != MPI_PROC_NULL ) {
/*
-------------------------------------
initialize the object from its buffer
-------------------------------------
*/
DenseMtx_initFromBuffer(inmtx) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n DenseMtx object intialized from its buffer") ;
fflush(msgFile) ;
}
if ( msglvl > 4 ) {
DenseMtx_writeForHumanEye(inmtx, msgFile) ;
fflush(msgFile) ;
}
}
if ( incount > 0 ) {
if ( nkeep + incount > nowned ) {
fprintf(msgFile, "\n fatal error in DenseMtx_splitByRows()"
"\n nkeep = %d, nrecv = %d, nowned = %d",
nkeep, nrecv, nowned) ;
exit(-1) ;
}
for ( irow = 0 ; irow < incount ; irow++, nkeep++ ) {
DenseMtx_copyRowAndIndex(keepmtx, nkeep, inmtx, irow) ;
}
}
}
/*
-------------------------
sort the rows and columns
-------------------------
*/
DenseMtx_sort(keepmtx) ;
/*
------------------------------------------------------
check that the matrix contains only the rows it should
------------------------------------------------------
*/
nrow = keepmtx->nrow ;
rowind = keepmtx->rowind ;
for ( ii = 0 ; ii < nrow ; ii++ ) {
irow = rowind[ii] ;
if ( irow < 0 || irow >= neqns ) {
fprintf(stderr,
"\n process %d : local row %d, global row %d, neqns = %d\n",
myid, ii, irow, neqns) ;
exit(-1) ;
}
if ( rowmap[irow] != myid ) {
fprintf(stderr,
"\n process %d : local row %d, global row %d, map = %d\n",
myid, ii, irow, rowmap[irow]) ;
exit(-1) ;
}
}
/*
------------------------
free the working storage
------------------------
*/
if ( outmtx != NULL ) {
DenseMtx_free(outmtx) ;
}
if ( inmtx != NULL ) {
DenseMtx_free(inmtx) ;
}
IVfree(rowsToSend) ;
IVfree(head) ;
IVfree(link) ;
return(keepmtx) ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------------------
purpose -- to scatter a DenseMtx object by rows
Xglobal -- global DenseMtx object, significant only for root
Xlocal -- local DenseMtx object, if NULL on input, will
be created if necessary
rowmapIV -- map from rows to owning processes
firsttag -- first tag to be used in these messages
stats[4] -- statistics vector
stats[0] -- # of messages sent
stats[1] -- # of messages received
stats[2] -- # of bytes sent
stats[3] -- # of bytes received
msglvl -- message level
msgFile -- message file
comm -- MPI communicator
return value -- Xlocal, a local DenseMtx object
created -- 98may16, cca
modified -- 98sep26, cca
mtx is not modified
-----------------------------------------------------------------
*/
DenseMtx *
DenseMtx_MPI_splitFromGlobalByRows (
DenseMtx *Xglobal,
DenseMtx *Xlocal,
IV *rowmapIV,
int root,
int stats[],
int msglvl,
FILE *msgFile,
int firsttag,
MPI_Comm comm
) {
DenseMtx *tempmtx ;
double *buffer ;
int ii, iproc, irow, maxnrow, myid, ncolX, nkeep,
nproc, nrowloc, nrowmap, nrowX, nsend, size, type ;
int *counts, *head, *link, *rowind, *rowmap ;
MPI_Status status ;
/*
-------------------------------------------------
get id of self, # of processes and # of equations
-------------------------------------------------
*/
MPI_Comm_rank(comm, &myid) ;
MPI_Comm_size(comm, &nproc) ;
if ( root < 0 || root >= nproc ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n root = %d, nproc = %d\n", root, nproc) ;
MPI_Finalize() ;
exit(-1) ;
}
/*--------------------------------------------------------------------*/
/*
--------------
check the data
--------------
*/
{
int rc = 1 ;
int *rcs = IVinit(nproc, -1) ;
if ( myid == root ) {
if ( Xglobal == NULL ) {
fprintf(stderr,
"\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
"\n Xglobal is NULL\n") ;
rc = -1 ;
}
if ( rowmapIV == NULL ) {
fprintf(stderr,
"\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
"\n rowmapIV is NULL\n") ;
rc = -2 ;
}
}
if ( msglvl > 0 && msgFile == NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
"\n msglvl > 0 and msgFile = NULL\n") ;
rc = -3 ;
}
if ( firsttag < 0 ) {
fprintf(stderr,
"\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
"\n firsttag = %d\n", firsttag) ;
rc = -4 ;
}
MPI_Allgather((void *) &rc, 1, MPI_INT,
(void *) rcs, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( rcs[iproc] != 1 ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_splitFromGlobalByRows()"
"\n trouble with return code") ;
IVfprintf(msgFile, nproc, rcs) ;
MPI_Finalize() ;
exit(rc) ;
}
}
}
IVfree(rcs) ;
}
/*--------------------------------------------------------------------*/
if ( myid == root ) {
type = Xglobal->type ;
IV_sizeAndEntries(rowmapIV, &nrowmap, &rowmap) ;
DenseMtx_dimensions(Xglobal, &nrowX, &ncolX) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitFromGlobalByRows"
"\n nproc = %d, myid = %d, nrowmap = %d, nrowX = %d, ncolX = %d",
nproc, myid, nrowmap, nrowX, ncolX) ;
fflush(msgFile) ;
}
}
/*
----------------------------------------
broadcast the type of entries and number
of right hand sides to all processors
----------------------------------------
*/
MPI_Bcast((void *) &type, 1, MPI_INT, root, comm) ;
MPI_Bcast((void *) &ncolX, 1, MPI_INT, root, comm) ;
stats[0] += 2 ;
stats[1] += 2 ;
stats[2] += 2*sizeof(int) ;
stats[3] += 2*sizeof(int) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n inside DenseMtx_MPI_splitFromGlobalByRows"
"\n type %d, ncolX %d", type, ncolX) ;
fflush(msgFile) ;
}
if ( myid == root ) {
/*
------------------------------------------------
create a head/link structure for the matrix rows
------------------------------------------------
*/
DenseMtx_rowIndices(Xglobal, &nrowX, &rowind) ;
counts = IVinit(nproc, 0) ;
head = IVinit(nproc, -1) ;
link = IVinit(nrowX, -1) ;
for ( ii = nrowX - 1 ; ii >= 0 ; ii-- ) {
irow = rowind[ii] ;
iproc = rowmap[irow] ;
link[ii] = head[iproc] ;
head[iproc] = ii ;
counts[iproc]++ ;
}
} else {
counts = NULL ;
}
/*
-------------------------------------------------
communicate the number of rows for each processor
-------------------------------------------------
*/
MPI_Scatter((void *) counts, 1, MPI_INT,
(void *) &nrowloc, 1, MPI_INT, root, comm) ;
stats[0] += 1 ;
stats[1] += 1 ;
stats[2] += (nproc-1)*sizeof(int) ;
stats[3] += (nproc-1)*sizeof(int) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n nrowloc = %d", nrowloc) ;
fflush(msgFile) ;
}
/*
-----------------------------------------------------------
if necessary, create the local Xloc matrix, then initialize
-----------------------------------------------------------
*/
if ( Xlocal == NULL ) {
Xlocal = DenseMtx_new() ;
}
DenseMtx_init(Xlocal, type, -1, -1, nrowloc, ncolX, 1, nrowloc) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n Xlocal after initialization") ;
DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
fflush(msgFile) ;
}
if ( myid == root ) {
/*
---------------------------------
load local matrix with owned rows
---------------------------------
*/
if ( nrowloc > 0 ) {
int iglob, iloc = 0 ;
for ( iglob = head[root] ; iglob != -1 ; iglob = link[iglob] ) {
DenseMtx_copyRowAndIndex(Xlocal, iloc, Xglobal, iglob) ;
iloc++ ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n Xlocal on root") ;
DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
fflush(msgFile) ;
}
}
/*
-----------------------------------------------------
create a temporary matrix to send to other processors
-----------------------------------------------------
*/
counts[myid] = 0 ;
maxnrow = IVmax(nproc, counts, &iproc) ;
if ( maxnrow > 0 ) {
DenseMtx *tempmtx = DenseMtx_new() ;
DenseMtx_init(tempmtx, type, -1, -1, maxnrow, ncolX, 1, maxnrow) ;
/*
----------------------------------
loop over the processors
collect owned rows into tempmtx
send tempmtx to processor
----------------------------------
*/
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( iproc != root && (nrowloc = counts[iproc]) > 0 ) {
DenseMtx_init(tempmtx, type, -1, -1,
nrowloc, ncolX, 1, nrowloc) ;
nsend = 0 ;
for ( ii = head[iproc] ; ii != -1 ; ii = link[ii] ) {
DenseMtx_copyRowAndIndex(tempmtx, nsend, Xglobal, ii) ;
nsend++ ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n tempmtx for proc %d", iproc) ;
DenseMtx_writeForHumanEye(tempmtx, msgFile) ;
fflush(msgFile) ;
}
size = DV_size(&tempmtx->wrkDV) ;
buffer = DV_entries(&tempmtx->wrkDV) ;
MPI_Send((void *) buffer, size, MPI_DOUBLE,
iproc, firsttag, comm) ;
stats[0] += 1 ;
stats[2] += size*sizeof(double) ;
}
}
DenseMtx_free(tempmtx) ;
}
/*
------------------------
free the working storage
------------------------
*/
IVfree(head) ;
IVfree(link) ;
IVfree(counts) ;
} else {
/*
------------------
non-root processor
------------------
*/
if ( nrowloc > 0 ) {
size = DV_size(&Xlocal->wrkDV) ;
buffer = DV_entries(&Xlocal->wrkDV) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n size = %d, buffer = %p", size, buffer) ;
fflush(msgFile) ;
}
MPI_Recv((void *) buffer, size, MPI_DOUBLE,
root, firsttag, comm, &status);
stats[1] += 1 ;
stats[3] += size*sizeof(double) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n Xlocal rec'd from root %d", root) ;
fflush(msgFile) ;
}
DenseMtx_initFromBuffer(Xlocal) ;
if ( msglvl > 2 ) {
DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
fflush(msgFile) ;
}
} else {
Xlocal = NULL ;
}
}
if ( msglvl > 3 ) {
if ( Xlocal != NULL ) {
fprintf(msgFile, "\n\n Xlocal") ;
DenseMtx_writeForHumanEye(Xlocal, msgFile) ;
} else {
fprintf(msgFile, "\n\n Xlocal is NULL") ;
}
fflush(msgFile) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n leaving DenseMtx_splitFromGlobalByRows()") ;
fflush(msgFile) ;
}
return(Xlocal) ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------------------
purpose -- to merge a DenseMtx object by rows
Xlocal -- DenseMtx object, can be NULL
Xglobal -- DenseMtx object, can be NULL
significant only for root
firsttag -- first tag to be used in these messages
stats[4] -- statistics vector
stats[0] -- # of messages sent
stats[1] -- # of messages received
stats[2] -- # of bytes sent
stats[3] -- # of bytes received
msglvl -- message level
msgFile -- message file
comm -- MPI communicator
return value --
if processor is root
Xglobal is returned, if was NULL on input, it is created
else
NULL
endif
Xlocal is not modified
created -- 98sep27, cca
-----------------------------------------------------------------
*/
DenseMtx *
DenseMtx_MPI_mergeToGlobalByRows (
DenseMtx *Xglobal,
DenseMtx *Xlocal,
int root,
int stats[],
int msglvl,
FILE *msgFile,
int firsttag,
MPI_Comm comm
) {
double *buffer ;
int iproc, maxnrow, myid, ncolX, nproc, nrowXloc, size, type ;
int *incounts ;
/*
-------------------------------------------------
get id of self, # of processes and # of equations
-------------------------------------------------
*/
MPI_Comm_rank(comm, &myid) ;
MPI_Comm_size(comm, &nproc) ;
if ( root < 0 || root >= nproc ) {
fprintf(stderr, "\n fatal error in DenseMtx_MPI_splitByRows()"
"\n root = %d, nproc = %d\n", root, nproc) ;
MPI_Finalize() ;
exit(-1) ;
}
/*--------------------------------------------------------------------*/
/*
--------------
check the data
--------------
*/
{
int rc = 1 ;
int *ivec = IVinit(nproc, -1) ;
if ( msglvl > 0 && msgFile == NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
"\n msglvl > 0 and msgFile = NULL\n") ;
rc = -2 ;
}
if ( firsttag < 0 ) {
fprintf(stderr,
"\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
"\n firsttag = %d\n", firsttag) ;
rc = -3 ;
}
MPI_Allgather((void *) &rc, 1, MPI_INT,
(void *) ivec, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( ivec[iproc] != 1 ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
"\n trouble with return code") ;
IVfprintf(msgFile, nproc, ivec) ;
MPI_Finalize() ;
exit(rc) ;
}
}
}
/*
------------------------------
check the type of the matrices
------------------------------
*/
type = (Xlocal != NULL) ? Xlocal->type : -1 ;
MPI_Allgather((void *) &type, 1, MPI_INT,
(void *) ivec, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( ivec[iproc] != -1 ) {
if ( type == -1 ) {
type = ivec[iproc] ;
} else if ( type != ivec[iproc] ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
"\n trouble with types\n") ;
IVfprintf(msgFile, nproc, ivec) ;
MPI_Finalize() ;
exit(-1) ;
}
}
}
}
/*
--------------------------------------
check the # of columns of the matrices
--------------------------------------
*/
ncolX = (Xlocal != NULL) ? Xlocal->ncol : 0 ;
MPI_Allgather((void *) &ncolX, 1, MPI_INT,
(void *) ivec, 1, MPI_INT, comm) ;
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( ivec[iproc] != 0 ) {
if ( ncolX == 0 ) {
ncolX = ivec[iproc] ;
} else if ( ncolX != ivec[iproc] ) {
if ( msgFile != NULL ) {
fprintf(msgFile,
"\n fatal error in DenseMtx_MPI_mergeToGlobalByRows()"
"\n trouble with ncolX\n") ;
IVfprintf(msgFile, nproc, ivec) ;
MPI_Finalize() ;
exit(-1) ;
}
}
}
}
IVfree(ivec) ;
}
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------
gather the number of incoming rows onto processor root
------------------------------------------------------
*/
nrowXloc = (Xlocal != NULL) ? Xlocal->nrow : 0 ;
if ( myid == root ) {
incounts = IVinit(nproc, 0) ;
} else {
incounts = NULL ;
}
MPI_Gather(&nrowXloc, 1, MPI_INT,
(void *) incounts, 1, MPI_INT, root, comm) ;
if ( myid == root ) {
DenseMtx *tempmtx ;
int count, iglob, iloc, incount, nrowXglobal ;
MPI_Status status ;
/*
----------------------------------------------------------
if necessary, create the global matrix, then initialize it
----------------------------------------------------------
*/
nrowXglobal = IVsum(nproc, incounts) ;
if ( Xglobal == NULL ) {
Xglobal = DenseMtx_new() ;
}
DenseMtx_init(Xglobal, type, -1, -1,
nrowXglobal, ncolX, 1, nrowXglobal) ;
/*
---------------------------------
load local matrix with owned rows
---------------------------------
*/
for ( iloc = iglob = 0 ; iloc < nrowXloc ; iloc++, iglob++ ) {
DenseMtx_copyRowAndIndex(Xglobal, iglob, Xlocal, iloc) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n after loading Xlocal on root") ;
DenseMtx_writeForHumanEye(Xglobal, msgFile) ;
fflush(msgFile) ;
}
/*
----------------------------------------------------------
create a temporary matrix to receive from other processors
----------------------------------------------------------
*/
incounts[root] = 0 ;
maxnrow = IVmax(nproc, incounts, &iproc) ;
tempmtx = DenseMtx_new() ;
DenseMtx_init(tempmtx, type, -1, -1, maxnrow, ncolX, 1, maxnrow) ;
size = DV_size(&tempmtx->wrkDV) ;
buffer = DV_entries(&tempmtx->wrkDV) ;
/*
----------------------------
loop over the processors
receive rows into tempmtx
----------------------------
*/
for ( iproc = 0 ; iproc < nproc ; iproc++ ) {
if ( iproc != root && (incount = incounts[iproc]) > 0 ) {
MPI_Recv((void *) buffer, size, MPI_DOUBLE,
iproc, firsttag, comm, &status) ;
MPI_Get_count(&status, MPI_DOUBLE, &count) ;
stats[1] += 1 ;
stats[3] += count*sizeof(double) ;
DenseMtx_initFromBuffer(tempmtx) ;
for ( iloc = 0 ; iloc < incount ; iloc++, iglob++ ) {
DenseMtx_copyRowAndIndex(Xglobal, iglob, tempmtx, iloc) ;
}
}
}
/*
------------------------
free the working storage
------------------------
*/
IVfree(incounts) ;
DenseMtx_free(tempmtx) ;
} else {
/*
------------------
non-root processor
------------------
*/
if ( nrowXloc > 0 ) {
size = DV_size(&Xlocal->wrkDV) ;
buffer = DV_entries(&Xlocal->wrkDV) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n size = %d, buffer = %p", size, buffer) ;
fflush(msgFile) ;
}
MPI_Send((void *) buffer, size, MPI_DOUBLE,
root, firsttag, comm) ;
stats[0] += 1 ;
stats[2] += size*sizeof(double) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n Xlocal sent to root %d", root) ;
fflush(msgFile) ;
}
}
Xglobal = NULL ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n leaving DenseMtx_mergeToGlobalByRows()") ;
fflush(msgFile) ;
}
return(Xglobal) ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1