fei_EqnCommMgr.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_macros.hpp>
00010 
00011 #include <fei_defs.h>
00012 
00013 #include <fei_EqnCommMgr.hpp>
00014 
00015 #include <fei_CommUtils.hpp>
00016 
00017 #include <fei_ProcEqns.hpp>
00018 #include <fei_EqnBuffer.hpp>
00019 #include <fei_TemplateUtils.hpp>
00020 
00021 #include <algorithm>
00022 
00023 #undef fei_file
00024 #define fei_file "EqnCommMgr.cpp"
00025 #include <fei_ErrMacros.hpp>
00026 
00027 //-----Constructor--------------------------------------------------------------
00028 EqnCommMgr::EqnCommMgr(MPI_Comm comm, bool accumulate)
00029  : accumulate_(accumulate),
00030    localProc_(0),
00031    recvProcEqns_(NULL),
00032    exchangeIndicesCalled_(false),
00033    recvEqns_(NULL),
00034    solnValues_(),
00035    sendProcEqns_(NULL),
00036    sendEqns_(NULL),
00037    sendEqnSoln_(),
00038    essBCEqns_(NULL),
00039    comm_(comm)
00040 {
00041   localProc_ = fei::localProc(comm_);
00042   recvEqns_ = new EqnBuffer();
00043   sendEqns_ = new EqnBuffer();
00044   recvProcEqns_ = new ProcEqns();
00045   sendProcEqns_ = new ProcEqns();
00046 
00047   essBCEqns_ = new EqnBuffer();
00048 }
00049 
00050 //-----CopyConstructor--------------------------------------------------------------
00051 EqnCommMgr::EqnCommMgr(const EqnCommMgr& src)
00052  : accumulate_(src.accumulate_),
00053    localProc_(0),
00054    recvProcEqns_(NULL),
00055    exchangeIndicesCalled_(false),
00056    recvEqns_(NULL),
00057    solnValues_(),
00058    sendProcEqns_(NULL),
00059    sendEqns_(NULL),
00060    sendEqnSoln_(),
00061    essBCEqns_(NULL),
00062    comm_(src.comm_)
00063 {
00064   *this = src;
00065 }
00066 
00067 //-----Destructor---------------------------------------------------------------
00068 EqnCommMgr::~EqnCommMgr()
00069 {
00070    delete recvEqns_;
00071    delete sendEqns_;
00072    delete recvProcEqns_;
00073    delete sendProcEqns_;
00074 
00075    delete essBCEqns_;
00076 }
00077 
00078 //------------------------------------------------------------------------------
00079 EqnCommMgr& EqnCommMgr::operator=(const EqnCommMgr& src)
00080 {
00081    delete recvProcEqns_;
00082    recvProcEqns_ = src.recvProcEqns_->deepCopy();
00083 
00084    exchangeIndicesCalled_ = src.exchangeIndicesCalled_;
00085    accumulate_ = src.accumulate_;
00086 
00087    delete recvEqns_;
00088    recvEqns_ = src.recvEqns_->deepCopy();
00089 
00090    int len = src.solnValues_.size();
00091 
00092    solnValues_.resize(len);
00093 
00094    for(int i=0; i<len; i++) {
00095       solnValues_[i] = src.solnValues_[i];
00096    }
00097 
00098    delete sendProcEqns_;
00099    sendProcEqns_ = src.sendProcEqns_->deepCopy();
00100 
00101    delete sendEqns_;
00102    sendEqns_ = src.sendEqns_->deepCopy();
00103 
00104    len = src.sendEqnSoln_.size();
00105    sendEqnSoln_.resize(len);
00106 
00107    for(int i=0; i<len; i++) {
00108       sendEqnSoln_[i] = src.sendEqnSoln_[i];
00109    }
00110 
00111    delete essBCEqns_;
00112    essBCEqns_ = src.essBCEqns_->deepCopy();
00113 
00114    return(*this);
00115 }
00116 
00117 //------------------------------------------------------------------------------
00118 EqnCommMgr* EqnCommMgr::deepCopy()
00119 {
00120 //
00121 //This function is like a copy-constructor. This function produces not just a
00122 //structural copy of the current object, it copies EVERYTHING, including data
00123 //contents.
00124 //
00125    EqnCommMgr* dest = new EqnCommMgr(comm_);
00126    *dest = *this;
00127    return(dest);
00128 }
00129 
00130 //------------------------------------------------------------------------------
00131 void EqnCommMgr::addLocalEqn(int eqnNumber, int srcProc) {
00132 //
00133 //This function adds the eqnNumber, srcProc pair to the recvProcsEqns_
00134 //object, which does the right thing as far as only putting it in lists that
00135 //it isn't already in, preserving order, etc. 
00136 //
00137    if (srcProc == localProc_) {
00138       FEI_CERR << "EqnCommMgr::addRecvEqn: ERROR, srcProc == localProc_, "
00139           << "which is a recipe for a deadlock." << FEI_ENDL;
00140       std::abort();
00141    }
00142 
00143    recvProcEqns_->addEqn(eqnNumber, srcProc);
00144 }
00145 
00146 //------------------------------------------------------------------------------
00147 void EqnCommMgr::addSolnValues(int* eqnNumbers, double* values, int num)
00148 {
00149   if (!exchangeIndicesCalled_) {
00150     FEI_CERR << "EqnCommMgr::addSolnValues: ERROR, you may not call this until"
00151       " after exchangeIndices has been called." << FEI_ENDL;
00152     std::abort();
00153   }
00154 
00155   std::vector<int>& recvEqnNumbers = recvEqns_->eqnNumbers();
00156 
00157   double* solnValuesPtr = &solnValues_[0];
00158   for(int i=0; i<num; i++) {
00159     int index = fei::binarySearch(eqnNumbers[i], recvEqnNumbers);
00160 
00161     if (index < 0) continue;
00162 
00163     solnValuesPtr[index] = values[i];
00164   }
00165 }
00166 
00167 //------------------------------------------------------------------------------
00168 int EqnCommMgr::exchangeIndices(FEI_OSTREAM* dbgOut) {
00169 //
00170 //This function performs the communication necessary to exchange remote
00171 //contributions to the matrix structure (indices only, not coefficients)
00172 //among all participating processors.
00173 //
00174 //Most of this function is #ifdef'd according to whether FEI_SER is
00175 //defined...
00176 //
00177 #ifndef FEI_SER
00178 
00179   int numSendEqns = sendEqns_->getNumEqns();
00180   fei::CSVec** sendEqnsPtr = numSendEqns>0 ?&(sendEqns_->eqns()[0]) : NULL;
00181   std::vector<int>& sendEqnNumbers = sendEqns_->eqnNumbers();
00182   std::vector<int> sendEqnLengths(numSendEqns);
00183   for(int i=0; i<numSendEqns; ++i) {
00184     sendEqnLengths[i] = sendEqnsPtr[i]->size();
00185   }
00186 
00187   if (sendEqnNumbers.size() > 0) {
00188     sendProcEqns_->setProcEqnLengths(&sendEqnNumbers[0],
00189                                      &sendEqnLengths[0],
00190                                      numSendEqns);
00191   }
00192 
00193   CHK_ERR( mirrorProcEqns(*sendProcEqns_, *recvProcEqns_) );
00194   CHK_ERR( mirrorProcEqnLengths(*sendProcEqns_, *recvProcEqns_) );
00195 
00196   sendEqns_->setNumRHSs(sendEqns_->getNumRHSs());
00197   recvEqns_->setNumRHSs(sendEqns_->getNumRHSs());
00198 
00199   size_t numRecvProcs = recvProcEqns_->getNumProcs();
00200   size_t numSendProcs = sendProcEqns_->getNumProcs();
00201   std::vector<int>& sendProcs       = sendProcEqns_->procsPtr();
00202 
00203   //while we're here, let's allocate the array into which we will (later)
00204   //recv the soln values for the remote equations we've contributed to.
00205   sendEqnSoln_.resize(numSendEqns);
00206 
00207   //First we'll figure out the expected receive-lengths and the expected
00208   //send lengths.
00209 
00210   std::vector<int> recvProcTotalLengths(numRecvProcs);
00211   std::vector<int> sendProcTotalLengths(numSendProcs);
00212 
00213   std::vector<int>& recvProcs = recvProcEqns_->procsPtr();
00214   std::vector<int>& eqnsPerRecvProc = recvProcEqns_->eqnsPerProcPtr();
00215   std::vector<std::vector<int>*>& recvProcEqnLengths =
00216     recvProcEqns_->procEqnLengthsPtr();
00217   std::vector<std::vector<int>*>& recvProcEqnNumbers = 
00218     recvProcEqns_->procEqnNumbersPtr();
00219 
00220   for(unsigned i=0; i<numRecvProcs; i++) {
00221 
00222     //first we need the total of eqn-lengths for this recv-proc
00223     int totalLength = 0;
00224     for(int j=0; j<eqnsPerRecvProc[i]; j++) {
00225       totalLength += (*(recvProcEqnLengths[i]))[j];
00226     }
00227     recvProcTotalLengths[i] = totalLength;
00228   }
00229 
00230   std::vector<int>& eqnsPerSendProc = sendProcEqns_->eqnsPerProcPtr();
00231   std::vector<std::vector<int>*>& sendProcEqnNumbers = 
00232     sendProcEqns_->procEqnNumbersPtr();
00233   std::vector<std::vector<int>*>& sendProcLengths =
00234     sendProcEqns_->procEqnLengthsPtr();
00235 
00236   for(unsigned i=0; i<numSendProcs; i++) {
00237     int totalLength = 0;
00238     for(int j=0; j<eqnsPerSendProc[i]; j++) {
00239       totalLength += (*(sendProcLengths[i]))[j];
00240     }
00241     sendProcTotalLengths[i] = totalLength;
00242   }
00243 
00244   //Before we get too carried away here, lets make sure that the messages we
00245   //expect to receive from each other processor are the same length as the
00246   //other processor plans to send. (Many times I've had crashes in MPI_Wait
00247   //below, which are always mysterious to figure out. They usually result
00248   //from mis-matched send/recv lengths.)
00249 
00250   CHK_ERR( consistencyCheck("exchangeIndices", recvProcs, recvProcTotalLengths,
00251                             sendProcs, sendProcTotalLengths) );
00252 
00253   int** recvProcEqnIndices = NULL;
00254 
00255   if (numRecvProcs > 0) {
00256     recvProcEqnIndices = new int*[numRecvProcs];
00257   }
00258 
00259   MPI_Request* indRequests = NULL;
00260 
00261   if (numRecvProcs > 0) {
00262     indRequests = new MPI_Request[numRecvProcs];
00263   }
00264 
00265   int numRecvsStarted = 0;
00266   int indTag = 199904;
00267 
00268   //now, let's start the recvs for the incoming indices.
00269   for(unsigned i=0; i<numRecvProcs; i++) {
00270 
00271     int totalLength = recvProcTotalLengths[i];
00272 
00273     recvProcEqnIndices[i] = new int[totalLength];
00274 
00275     //launch the recvs for the indices now.
00276     if (MPI_Irecv(recvProcEqnIndices[i], totalLength, MPI_INT,
00277                   recvProcs[i], indTag, comm_, &indRequests[i]) != MPI_SUCCESS) {
00278        ERReturn(-1);
00279     }
00280 
00281     numRecvsStarted++;
00282   }
00283 
00284   //now we need to build the lists of outgoing indices and send those.
00285   std::vector<int> indices;
00286 
00287   for(unsigned i=0; i<numSendProcs; i++) {
00288     int totalLength = sendProcTotalLengths[i];
00289     int j;
00290 
00291     indices.resize(totalLength);
00292     int* indicesPtr = &indices[0];
00293 
00294     int offset = 0;
00295 
00296     for(j=0; j<eqnsPerSendProc[i]; j++) {
00297       int eqnLoc = sendEqns_->getEqnIndex((*(sendProcEqnNumbers[i]))[j]);
00298       std::vector<int>& sendIndices = sendEqns_->eqns()[eqnLoc]->indices();
00299       int* sendIndicesPtr = &sendIndices[0];
00300 
00301       for(int k=0; k<(*(sendProcLengths[i]))[j]; k++) {
00302         indicesPtr[offset++] = sendIndicesPtr[k];
00303       }
00304     }
00305 
00306     if (MPI_Send(indicesPtr, totalLength, MPI_INT, sendProcs[i],
00307                  indTag, comm_) != MPI_SUCCESS) ERReturn(-1)
00308   }
00309 
00310   //and finally, we're ready to complete the irecvs for the indices and put
00311   //them away.
00312   int numCompleted = 0;
00313   for(unsigned i=0; i<numRecvProcs; i++) {
00314     MPI_Status status;
00315     int index = i;
00316     MPI_Wait(&indRequests[i], &status);
00317     numCompleted++;
00318 
00319     int offset = 0;
00320     for(int j=0; j<eqnsPerRecvProc[index]; j++) {
00321       int eqn = (*(recvProcEqnNumbers[index]))[j];
00322       int* indxs = &(recvProcEqnIndices[index][offset]);
00323       int len = (*(recvProcEqnLengths[index]))[j];
00324 
00325       recvEqns_->addIndices(eqn, indxs, len);
00326 
00327       offset += len;
00328     }
00329 
00330     delete [] recvProcEqnIndices[index];
00331   }
00332 
00333   delete [] recvProcEqnIndices;
00334   delete [] indRequests;
00335 
00336   if (numRecvsStarted != numCompleted) {
00337     FEI_CERR << "EqnCommMgr::exchangeIndices: recv-send mismatch; "
00338           << "numRecvsStarted: " << numRecvsStarted << ", numCompleted: "
00339          << numCompleted << FEI_ENDL;
00340     std::abort();
00341   }
00342 
00343   //allocate the solnValue_ list, which is of size numRecvEqns.
00344   int numRecvEqns = recvEqns_->getNumEqns();
00345   solnValues_.resize(numRecvEqns);
00346 
00347   exchangeIndicesCalled_ = true;
00348 
00349   if (dbgOut != NULL) {
00350     FEI_OSTREAM& os = *dbgOut;
00351     os << "#ereb exchangeIndices, sendEqns_:"<<FEI_ENDL;
00352 //    os << *sendEqns_<<FEI_ENDL;
00353   }
00354 
00355 #else
00356   (void)dbgOut;
00357 #endif // #ifndef FEI_SER
00358 
00359    return(0);
00360 }
00361 
00362 //------------------------------------------------------------------------------
00363 int EqnCommMgr::consistencyCheck(const char* caller,
00364          std::vector<int>& recvProcs,
00365          std::vector<int>& recvProcTotalLengths,
00366          std::vector<int>& sendProcs,
00367          std::vector<int>& sendProcTotalLengths)
00368 {
00369   int err = 0;
00370   //First we'll gather each processors send-lengths onto all other processors.
00371   std::vector<int> globalProcSendLengths, globalSendProcs;
00372   std::vector<int> gatherSizes;
00373 
00374   CHK_ERR( fei::Allgatherv(comm_, sendProcTotalLengths,
00375          gatherSizes, globalProcSendLengths) );
00376 
00377   CHK_ERR( fei::Allgatherv(comm_, sendProcs,
00378          gatherSizes, globalSendProcs) );
00379 
00380   //Now check the consistency of the global send-lengths against local
00381   //receive-lengths.
00382   int offset = 0;
00383   for(unsigned i=0; i<gatherSizes.size(); i++) {
00384     int size = gatherSizes[i];
00385     if ((int)i==localProc_) {offset += size; continue; }
00386 
00387     //now we're ready to stride through processor i's sendProcs. Only do this
00388     //if processor i is one of our recvProcs.
00389     std::vector<int>::iterator rp_iter = 
00390       std::lower_bound(recvProcs.begin(), recvProcs.end(), (int)i);
00391     int rpIndex = -1;
00392     if (rp_iter != recvProcs.end() && (int)i == *rp_iter) {
00393       rpIndex = (int)(rp_iter - recvProcs.begin());
00394     }
00395 
00396     if (rpIndex < 0) {
00397       // proc i is not one of our recvProcs. Let's make sure
00398       //that we're not one of proc i's sendProcs.
00399       for(int j=0; j<size; j++) {
00400   if (globalSendProcs[offset+j] == localProc_) {
00401     FEI_CERR << "EqnCommMgr::"<<caller<<" ERROR: proc " << localProc_
00402          << " is not expecting to receive from proc " << i << " but proc "
00403          << i << " is expecting to send to proc " << localProc_ << FEI_ENDL;
00404     err = -1;
00405   }
00406       }
00407       if (err == -1) break;
00408 
00409       //if err != -1, simply jump to the next for(i... iteration.
00410       offset += size;
00411       continue;
00412     }
00413 
00414     for(int j=0; j<size; j++) {
00415       if (globalSendProcs[offset+j] == localProc_) {
00416   int sendLength = globalProcSendLengths[offset+j];
00417   int recvLength = recvProcTotalLengths[rpIndex];
00418   if (sendLength != recvLength) {
00419     FEI_CERR << "EqnCommMgr::"<<caller<<" ERROR: proc " << localProc_
00420          << " is expecting to receive " << recvLength << " indices from "
00421          << "proc " << i << " but proc " << i << " is expecting to send "
00422          << sendLength << " indices to proc " << localProc_ << FEI_ENDL;
00423     err = -1;
00424   }
00425       }
00426     }
00427 
00428     offset += size;
00429   }
00430 
00431   int globalErr = 0;
00432   CHK_ERR( fei::GlobalSum(comm_, err, globalErr) );
00433 
00434   return(globalErr);
00435 }
00436 
00437 //------------------------------------------------------------------------------
00438 int EqnCommMgr::exchangeEqns(FEI_OSTREAM* dbgOut)
00439 {
00440   //
00441   //This function performs the communication necessary to exchange remote
00442   //equations (both indices and coefficients) among all participating processors.
00443   //
00444 
00445   recvEqns_->resetCoefs();
00446 
00447   if (dbgOut != NULL) {
00448     FEI_OSTREAM& os = *dbgOut;
00449     os << "#ereb exchangeEqns begin, sendEqns_:"<<FEI_ENDL;
00450 //    os << *sendEqns_<<FEI_ENDL;
00451   }
00452 
00453   CHK_ERR( exchangeEqnBuffers(comm_, sendProcEqns_,
00454             sendEqns_, recvProcEqns_, recvEqns_, accumulate_));
00455 
00456   if (dbgOut != NULL) {
00457     FEI_OSTREAM& os = *dbgOut;
00458     os << "#ereb exchangeEqns end, sendEqns_:"<<FEI_ENDL;
00459 //    os << *sendEqns_<<FEI_ENDL;
00460   }
00461 
00462   return(0);
00463 }
00464 
00465 //------------------------------------------------------------------------------
00466 int EqnCommMgr::exchangeEqnBuffers(MPI_Comm comm, ProcEqns* sendProcEqns,
00467                               EqnBuffer* sendEqns, ProcEqns* recvProcEqns,
00468                               EqnBuffer* recvEqns, bool accumulate)
00469 {
00470 //
00471 //This function performs the communication necessary to exchange remote
00472 //equations (both indices and coefficients) among all participating processors.
00473 //
00474 //Most of this function is #ifdef'd according to whether FEI_SER is
00475 //defined.
00476 #ifdef FEI_SER
00477   (void)comm;
00478   (void)sendProcEqns;
00479   (void)sendEqns;
00480   (void)recvProcEqns;
00481   (void)recvEqns;
00482   (void)accumulate;
00483 #else
00484    int indTag = 19991130, coefTag = 19991131;
00485 
00486    size_t numRecvProcs = recvProcEqns->getNumProcs();
00487    size_t numSendProcs = sendProcEqns->getNumProcs();
00488    if ((numRecvProcs == 0) && (numSendProcs == 0)) return(0);
00489 
00490    //we assume that sendProcEqns and recvProcEqns are fully populated with
00491    //equation-numbers and their lengths, and that this data is consistent with
00492    //the data in sendEqns...
00493 
00494    MPI_Request* indRequests = NULL;
00495    MPI_Request* coefRequests = NULL;
00496    int** recvProcEqnIndices = NULL;
00497    double** recvProcEqnCoefs = NULL;
00498 
00499    if (numRecvProcs > 0) {
00500       indRequests = new MPI_Request[numRecvProcs];
00501       coefRequests = new MPI_Request[numRecvProcs];
00502       recvProcEqnIndices = new int*[numRecvProcs];
00503       recvProcEqnCoefs = new double*[numRecvProcs];
00504    }
00505 
00506    int numRHSs = sendEqns->getNumRHSs();
00507 
00508    //now, let's allocate the space for the incoming equations.
00509    //each row of recvProcEqnIndices will be of length
00510    //sum-of-recvProcEqnLengths[i], and each row of recvProcEqnCoefs will be
00511    //of length sum-of-recvProcEqnLengths[i] + numRHSs*eqnsPerRecvProc[i].
00512 
00513    std::vector<int>& recvProcs = recvProcEqns->procsPtr();
00514    std::vector<int>& eqnsPerRecvProc = recvProcEqns->eqnsPerProcPtr();
00515    std::vector<std::vector<int>*>& recvProcEqnNumbers =
00516      recvProcEqns->procEqnNumbersPtr();
00517    std::vector<std::vector<int>*>& recvProcEqnLengths =
00518      recvProcEqns->procEqnLengthsPtr();
00519 
00520    int padding = 2;
00521    for(unsigned i=0; i<numRecvProcs; i++) {
00522       int totalLength = 0;
00523 
00524       for(int j=0; j<eqnsPerRecvProc[i]; j++) {
00525          totalLength += (*(recvProcEqnLengths[i]))[j];
00526       }
00527 
00528       //in case we're only exchanging rhs coefs, (i.e., there are no
00529       //column-indices) let's make the length totalLength+padding so we can
00530       //send the new*Data_ indicators.
00531       recvProcEqnIndices[i] = new int[totalLength+padding];
00532 
00533       int coefLength = totalLength + numRHSs*eqnsPerRecvProc[i];
00534 
00535       recvProcEqnCoefs[i] = new double[coefLength];
00536       //let's go ahead and launch the recvs for the indices and coefs now.
00537       MPI_Irecv(recvProcEqnIndices[i], totalLength+padding, MPI_INT,
00538                 recvProcs[i], indTag, comm, &indRequests[i]);
00539 
00540       MPI_Irecv(recvProcEqnCoefs[i], coefLength, MPI_DOUBLE,
00541                 recvProcs[i], coefTag, comm, &coefRequests[i]);
00542    }
00543 
00544    //ok, now we need to build the lists of outgoing indices and coefs, and
00545    //send those.
00546    std::vector<int>& sendProcs = sendProcEqns->procsPtr();
00547    std::vector<int>& eqnsPerSendProc = sendProcEqns->eqnsPerProcPtr();
00548    std::vector<std::vector<int>*>& sendProcEqnNumbers =
00549      sendProcEqns->procEqnNumbersPtr();
00550    std::vector<std::vector<int>*>& sendProcEqnLengths =
00551      sendProcEqns->procEqnLengthsPtr();
00552 
00553    std::vector<std::vector<double>*>& sendRHS = *(sendEqns->rhsCoefsPtr());
00554 
00555    for(unsigned i=0; i<numSendProcs; i++) {
00556       int totalLength = 0;
00557       int* sendProcEqnNumbers_i = &(*(sendProcEqnNumbers[i]))[0];
00558       int* sendProcEqnLengths_i = &(*(sendProcEqnLengths[i]))[0];
00559       int j;
00560       for(j=0; j<eqnsPerSendProc[i]; j++)
00561          totalLength += sendProcEqnLengths_i[j];
00562 
00563       //As with the recv code above, let's make the indices length
00564       //be totalLength+padding...
00565       std::vector<int> indices(totalLength+padding);
00566       int* indicesPtr = &indices[0];
00567       int coefLength = totalLength + numRHSs*eqnsPerSendProc[i];
00568 
00569       std::vector<double> coefs(coefLength);
00570       double* coefsPtr = &coefs[0];
00571 
00572       int offset = 0;
00573 
00574       //first pack up the coefs and indices
00575       for(j=0; j<eqnsPerSendProc[i]; j++) {
00576          int eqnLoc = sendEqns->getEqnIndex(sendProcEqnNumbers_i[j]);
00577    int* sendIndices = &(sendEqns->eqns()[eqnLoc]->indices())[0];
00578    double* sendCoefs= &(sendEqns->eqns()[eqnLoc]->coefs())[0];
00579 
00580          for(int k=0; k<sendProcEqnLengths_i[j]; k++) {
00581             indicesPtr[offset] = sendIndices[k];
00582             coefsPtr[offset++] = sendCoefs[k];
00583          }
00584       }
00585 
00586       //now append the new*Data_ indicators to the end of the indices array
00587       indicesPtr[offset] = sendEqns->newCoefData_;
00588       indicesPtr[offset+1] = sendEqns->newRHSData_;
00589 
00590       //now append the RHS coefs to the end of the coefs array
00591       for(j=0; j<eqnsPerSendProc[i]; j++) {
00592          int eqnLoc = sendEqns->getEqnIndex(sendProcEqnNumbers_i[j]);
00593 
00594          for(int k=0; k<numRHSs; k++) {
00595             coefsPtr[offset++] = (*(sendRHS[eqnLoc]))[k];
00596          }
00597       }
00598 
00599       MPI_Send(&indices[0], (int)indices.size(), MPI_INT, sendProcs[i],
00600                indTag, comm);
00601       MPI_Send(&coefs[0], coefLength, MPI_DOUBLE, sendProcs[i],
00602          coefTag, comm);
00603    }
00604 
00605    //and finally, we're ready to complete the irecvs for the indices and coefs,
00606    //and put them away.
00607    for(unsigned i=0; i<numRecvProcs; i++) {
00608       MPI_Status status;
00609       int index = i;
00610       MPI_Wait(&indRequests[index], &status);
00611       MPI_Wait(&coefRequests[index], &status);
00612 
00613       int j, offset = 0;
00614       for(j=0; j<eqnsPerRecvProc[index]; j++) {
00615          int eqn = (*(recvProcEqnNumbers[index]))[j];
00616          int* indices = &(recvProcEqnIndices[index][offset]);
00617          double* coefs = &(recvProcEqnCoefs[index][offset]);
00618          int len = (*(recvProcEqnLengths[index]))[j];
00619 
00620          recvEqns->addEqn(eqn, coefs, indices, len, accumulate);
00621 
00622          offset += len;
00623       }
00624 
00625       recvEqns->newCoefData_ += recvProcEqnIndices[index][offset];
00626       recvEqns->newRHSData_  += recvProcEqnIndices[index][offset+1];
00627       delete [] recvProcEqnIndices[index];
00628    }
00629 
00630    //now unpack the RHS entries
00631 
00632    recvEqns->setNumRHSs(numRHSs);
00633 
00634    for(unsigned i=0; i<numRecvProcs; i++) {
00635       int j, offset = 0;
00636       for(j=0; j<eqnsPerRecvProc[i]; j++) {
00637          offset += (*(recvProcEqnLengths[i]))[j];
00638       }
00639 
00640       for(j=0; j<eqnsPerRecvProc[i]; j++) {
00641          int eqn = (*(recvProcEqnNumbers[i]))[j];
00642 
00643          for(int k=0; k<numRHSs; k++) {
00644             CHK_ERR( recvEqns->addRHS(eqn, k, recvProcEqnCoefs[i][offset++],
00645               accumulate));
00646          }
00647       }
00648 
00649       delete [] recvProcEqnCoefs[i];
00650    }
00651 
00652    delete [] recvProcEqnIndices;
00653    delete [] recvProcEqnCoefs;
00654    delete [] indRequests;
00655    delete [] coefRequests;
00656 
00657 #endif //#ifndef FEI_SER
00658 
00659    return(0);
00660 }
00661 
00662 //------------------------------------------------------------------------------
00663 void EqnCommMgr::exchangeSoln()
00664 {
00665   //Most of this function is #ifdef'd according to whether FEI_SER 
00666   //is defined...
00667 #ifndef FEI_SER
00668 
00669    int solnTag = 199906;
00670 
00671    MPI_Request* solnRequests = NULL;
00672    double** solnBuffer = NULL;
00673 
00674    size_t numSendProcs = sendProcEqns_->getNumProcs();
00675    std::vector<int>& sendProcs = sendProcEqns_->procsPtr();
00676    std::vector<int>& eqnsPerSendProc = sendProcEqns_->eqnsPerProcPtr();
00677 
00678    if (numSendProcs > 0) {
00679       solnRequests = new MPI_Request[numSendProcs];
00680       solnBuffer = new double*[numSendProcs];
00681    }
00682 
00683    MPI_Comm comm = comm_;
00684 
00685    //let's launch the recv's for the incoming solution values.
00686    for(unsigned i=0; i<numSendProcs; i++) {
00687       solnBuffer[i] = new double[eqnsPerSendProc[i]];
00688 
00689       MPI_Irecv(solnBuffer[i], eqnsPerSendProc[i], MPI_DOUBLE, sendProcs[i],
00690                 solnTag, comm, &solnRequests[i]);
00691    }
00692 
00693    size_t numRecvProcs = recvProcEqns_->getNumProcs();
00694    std::vector<int>& recvProcs = recvProcEqns_->procsPtr();
00695    std::vector<int>& eqnsPerRecvProc = recvProcEqns_->eqnsPerProcPtr();
00696    std::vector<std::vector<int>*>& recvProcEqnNumbers =
00697      recvProcEqns_->procEqnNumbersPtr();
00698 
00699    //now let's send the outgoing solutions.
00700    for(unsigned i=0; i<numRecvProcs; i++) {
00701       double* solnBuff = new double[eqnsPerRecvProc[i]];
00702 
00703       for(int j=0; j<eqnsPerRecvProc[i]; j++) {
00704   int eqnNumber = (*(recvProcEqnNumbers[i]))[j];
00705 
00706   int index = recvEqns_->getEqnIndex(eqnNumber);
00707   solnBuff[j] = solnValues_[index];
00708       }
00709 
00710       MPI_Send(solnBuff, eqnsPerRecvProc[i], MPI_DOUBLE, recvProcs[i],
00711                solnTag, comm);
00712 
00713       delete [] solnBuff;
00714    }
00715 
00716    std::vector<std::vector<int>*>& sendProcEqnNumbers =
00717      sendProcEqns_->procEqnNumbersPtr();
00718 
00719    //make sure the sendEqnSoln_ array is big enough...
00720    sendEqnSoln_.resize(sendEqns_->getNumEqns());
00721 
00722    //ok, complete the above recvs and store the soln values.
00723    for(unsigned i=0; i<numSendProcs; i++) {
00724       int index;
00725       MPI_Status status;
00726       MPI_Waitany((int)numSendProcs, solnRequests, &index, &status);
00727 
00728       for(int j=0; j<eqnsPerSendProc[index]; j++) {
00729   int eqnNumber = (*(sendProcEqnNumbers[index]))[j];
00730   int ind = sendEqns_->getEqnIndex(eqnNumber);
00731 
00732   sendEqnSoln_[ind] = solnBuffer[index][j];
00733       }
00734 
00735       delete [] solnBuffer[index];
00736    }
00737 
00738    delete [] solnRequests;
00739    delete [] solnBuffer;
00740 #endif //#ifndef FEI_SER
00741 }
00742 
00743 //------------------------------------------------------------------------------
00744 int EqnCommMgr::mirrorProcEqns(ProcEqns& inProcEqns, ProcEqns& outProcEqns)
00745 {
00746   //Beginning assumption: we (the local processor) have a populated ProcEqns
00747   //object (inProcEqns) which contains information pairing certain equations
00748   //with certain remote processors. These can be equations that we will be
00749   //receiving in an all-to-all exchange, or equations that we will be sending in
00750   //an all-to-all exchange. In either case, the "mirror" of that information is
00751   //needed before the exchange can be performed. i.e., if we know which eqns
00752   //we'll be sending, then the mirror info concerns the eqns we'll be recv'ing,
00753   //and vice-versa if we already know which eqns we'll be recv'ing.
00754   //
00755   //This function is to obtain that mirror info, and return it in outProcEqns.
00756   //
00757   //Given a populated ProcEqns object, we want to populate the mirror ProcEqns
00758   //object.
00759   //
00760   //First figure out how many procs belong in outProcEqns.
00761   //Then receive, from each of those procs, the list of associated equations.
00762   //Then we'll have the info necessary to populate the 'outProcEqns' object.
00763 
00764   //Most of this function is #ifdef'd according to whether FEI_SER 
00765   //is defined.
00766 #ifdef FEI_SER
00767   (void)inProcEqns;
00768   (void)outProcEqns;
00769 #else
00770   int numProcs = fei::numProcs(comm_);
00771 
00772   if (numProcs < 2) return(0);
00773 
00774   std::vector<int> buf(numProcs*2, 0);
00775 
00776   std::vector<int>& inProcs = inProcEqns.procsPtr();
00777   std::vector<int> outProcs;
00778   fei::mirrorProcs(comm_, inProcs, outProcs);
00779 
00780   std::vector<int>& eqnsPerInProc = inProcEqns.eqnsPerProcPtr();
00781 
00782   size_t numOutProcs = outProcs.size();
00783 
00784   std::vector<int> recvbuf(numOutProcs, 0);
00785 
00786   //now send a length (the contents of buf[i]) to each "in-proc"
00787   //(length of the list of equation data that is to follow).
00788   MPI_Request* requests = new MPI_Request[numOutProcs];
00789   MPI_Status* statuses = new MPI_Status[numOutProcs];
00790   int firsttag = 1014;
00791   int offset = 0;
00792   for(size_t i=0; i<numOutProcs; ++i) {
00793     if (MPI_Irecv(&(recvbuf[i]), 1, MPI_INT, outProcs[i], firsttag,
00794       comm_, &requests[offset++]) != MPI_SUCCESS) ERReturn(-1);
00795   }
00796 
00797   size_t numInProcs = inProcs.size();
00798   for(size_t i=0; i<numInProcs; ++i) {
00799     if (MPI_Send(&(eqnsPerInProc[i]), 1, MPI_INT, inProcs[i], firsttag,
00800      comm_) != MPI_SUCCESS) ERReturn(-1);
00801   }
00802 
00803   MPI_Waitall((int)numOutProcs, requests, statuses);
00804 
00805   delete [] requests;
00806   delete [] statuses;
00807 
00808   std::vector<int> lengths(numOutProcs);
00809 
00810   offset = 0;
00811   for(unsigned i=0; i<numOutProcs; ++i) {
00812     if (recvbuf[i] > 0) {
00813       lengths[offset++] = recvbuf[i];
00814     }
00815   }
00816 
00817   //now we need to create 'numOutProcs' lists, into which we'll receive the
00818   //equation-numbers that those procs send to us.
00819   std::vector<std::vector<int>*>* outEqns = NULL;
00820   if (numOutProcs > 0) {
00821     outEqns = new std::vector<std::vector<int>*>(numOutProcs);
00822   }
00823 
00824   for(unsigned i=0; i<numOutProcs; i++) {
00825     (*outEqns)[i] = new std::vector<int>(lengths[i]);
00826   }
00827 
00828   //finally we're ready to exchange lists of equations.
00829 
00830   CHK_ERR( fei::exchangeData(comm_, inProcs,
00831                              inProcEqns.procEqnNumbersPtr(),
00832                              outProcs, true, *outEqns) );
00833 
00834   //now we've completed all the communication, so we're ready to put the data
00835   //we received into the outProcEqns object.
00836   for(unsigned i=0; i<numOutProcs; i++) {
00837     std::vector<int>* eqnArray = (*outEqns)[i];
00838     int* eqns = &(*eqnArray)[0];
00839     size_t len = eqnArray->size();
00840     for(unsigned j=0; j<len; j++) {
00841       outProcEqns.addEqn(eqns[j], outProcs[i]);
00842     }
00843     delete eqnArray;
00844   }
00845 
00846   delete outEqns;
00847 
00848 #endif //#ifndef FEI_SER
00849 
00850   return(0);
00851 }
00852 
00853 //------------------------------------------------------------------------------
00854 int EqnCommMgr::mirrorProcEqnLengths(ProcEqns& inProcEqns,
00855              ProcEqns& outProcEqns)
00856 {
00857   //Support function to set up information needed for exchanging equation
00858   //info among processors.
00859   //This function plays a similar role to that of the above 'mirrorProcEqns'
00860   //function, but exchanges the length information rather than the
00861   //equation-number information. THIS FUNCTION ASSUMES that both inProcEqns and
00862   //outProcEqns are already populated with eqn-number information.
00863 
00864   //Most of this function is #ifdef'd according to whether FEI_SER 
00865   //is defined.
00866 #ifdef FEI_SER
00867   (void)inProcEqns;
00868   (void)outProcEqns;
00869 #else
00870   if (fei::numProcs(comm_) == 1) return(0);
00871 
00872   CHK_ERR( fei::exchangeData( comm_, inProcEqns.procsPtr(),
00873             inProcEqns.procEqnLengthsPtr(),
00874             outProcEqns.procsPtr(),
00875             true,
00876             outProcEqns.procEqnLengthsPtr() ) );
00877 #endif //#ifndef FEI_SER
00878 
00879   return(0);
00880 }
00881 
00882 //------------------------------------------------------------------------------
00883 int EqnCommMgr::addRemoteEqn(int eqnNumber, int destProc,
00884                             const double* coefs, const int* indices, int num) {
00885    (void)destProc;
00886    sendEqns_->newCoefData_ = 1;
00887 
00888    return(sendEqns_->addEqn(eqnNumber, coefs, indices, num, accumulate_));
00889 }
00890 
00891 //------------------------------------------------------------------------------
00892 int EqnCommMgr::addRemoteEqn(int eqnNumber, const double* coefs,
00893            const int* indices, int num)
00894 {
00895    sendEqns_->newCoefData_ = 1;
00896 
00897    return(sendEqns_->addEqn(eqnNumber, coefs, indices, num, accumulate_));
00898 }
00899 
00900 //------------------------------------------------------------------------------
00901 void EqnCommMgr::setNumRHSs(int numRHSs) {
00902    sendEqns_->setNumRHSs(numRHSs);
00903    recvEqns_->setNumRHSs(numRHSs);
00904 }
00905 
00906 //------------------------------------------------------------------------------
00907 int EqnCommMgr::addRemoteRHS(int eqnNumber, int destProc, int rhsIndex,
00908                             double value)
00909 {
00910    (void)destProc;
00911    sendEqns_->newRHSData_ = 1;
00912    return(sendEqns_->addRHS(eqnNumber, rhsIndex, value));
00913 }
00914 
00915 //------------------------------------------------------------------------------
00916 int EqnCommMgr::addRemoteRHS(int eqnNumber, int rhsIndex, double value)
00917 {
00918    sendEqns_->newRHSData_ = 1;
00919    return(sendEqns_->addRHS(eqnNumber, rhsIndex, value));
00920 }
00921 
00922 //------------------------------------------------------------------------------
00923 void EqnCommMgr::addRemoteIndices(int eqnNumber, int destProc,
00924                                 int* indices, int num)
00925 {
00926   if (destProc < 0) {
00927     FEI_CERR << "fei: EqnCommMgr::addRemoteIndices ERROR, destProc < 0" << FEI_ENDL;
00928     std::abort();
00929   }
00930 
00931   sendEqns_->addIndices(eqnNumber, indices, num);
00932 
00933   sendProcEqns_->addEqn(eqnNumber, destProc);
00934 }
00935 
00936 //------------------------------------------------------------------------------
00937 void EqnCommMgr::resetCoefs() {
00938    recvEqns_->resetCoefs();
00939    sendEqns_->resetCoefs();
00940    essBCEqns_->resetCoefs();
00941 
00942    sendEqns_->newCoefData_ = 0;
00943    sendEqns_->newRHSData_ = 0;
00944    recvEqns_->newCoefData_ = 0;
00945    recvEqns_->newRHSData_ = 0;
00946 
00947    int numRecvEqns = recvEqns_->getNumEqns();
00948    for(int i=0; i<numRecvEqns; i++) {
00949       solnValues_[i] = 0.0;
00950    }
00951 }
00952 
00953 //----------------------------------------------------------------------------
00954 int EqnCommMgr::gatherSharedBCs(EqnBuffer& bcEqns)
00955 {
00956   //Gather boundary-condition equations from all sharing processors to the
00957   //owning processors.
00958   //
00959   //On entry to this function, bcEqns contains all boundary-condition equations
00960   //that were specified by the finite-element application on this processor.
00961   //On exit, bcEqns will also include any boundary-condition equations that were
00962   //specified by the finite-element application on processors that share nodes
00963   //that are owned by this processor.
00964   //
00965   //We'll set up two new equation buffers: sendBCs and recvBCs. From the input
00966   //bcEqns, we'll put each eqn that is in sendEqns_ into sendBCs because these
00967   //are the equations that are shared and remotely owned. We'll then do a
00968   //gather where all processors send their shared-bc eqns to the owner of those
00969   //equations. The result of this gather will be in recvBCs, which we will then
00970   //merge into bcEqns.
00971 
00972   if (fei::numProcs(comm_) == 1) return(0);
00973 
00974   int i;
00975   EqnBuffer sendBCs, recvBCs;
00976   ProcEqns sendBCProcEqns, recvBCProcEqns;
00977 
00978   //now loop through the equations in bcEqns, checking whether they're in
00979   //our sendEqns_ object.
00980   std::vector<int>& bcEqnNumbers = bcEqns.eqnNumbers();
00981   int numBCeqns = bcEqnNumbers.size();
00982 
00983   std::vector<int>& sendEqnNumbers = sendEqns_->eqnNumbers();
00984   std::vector<int>& sendProcs = sendProcEqns_->procsPtr();
00985   std::vector<std::vector<int>*>& sendProcEqnNumbers = 
00986     sendProcEqns_->procEqnNumbersPtr();
00987   size_t numSendProcs = sendProcs.size();
00988 
00989   for(i=0; i<numBCeqns; i++) {
00990     int eqn = bcEqnNumbers[i];
00991 
00992     int index = fei::binarySearch(eqn, sendEqnNumbers);
00993     if (index<0) continue;
00994 
00995     std::vector<double>& coefs = bcEqns.eqns()[i]->coefs();
00996     std::vector<int>& indices = bcEqns.eqns()[i]->indices();
00997     CHK_ERR( sendBCs.addEqn(eqn, &coefs[0], &indices[0],
00998           indices.size(), false) );
00999 
01000     for(unsigned p=0; p<numSendProcs; p++) {
01001       if (std::binary_search(sendProcEqnNumbers[p]->begin(),
01002                              sendProcEqnNumbers[p]->end(), eqn)) {
01003 
01004         sendBCProcEqns.addEqn(eqn, indices.size(), sendProcs[p]);
01005       }
01006     }
01007   }
01008 
01009   //now set up the required mirror info, then perform the exchange among procs.
01010   CHK_ERR( mirrorProcEqns(sendBCProcEqns, recvBCProcEqns) );
01011   CHK_ERR( mirrorProcEqnLengths(sendBCProcEqns, recvBCProcEqns) );
01012   CHK_ERR( exchangeEqnBuffers(comm_, &sendBCProcEqns,
01013             &sendBCs, &recvBCProcEqns, &recvBCs, false) );
01014 
01015   //finally, merge the recvBCs into the bcEqns buffer.
01016   CHK_ERR( bcEqns.addEqns(recvBCs, false) );
01017 
01018   return(0);
01019 }
01020 
01021 //------------------------------------------------------------------------------
01022 int EqnCommMgr::exchangeRemEssBCs(int* essEqns, int numEssEqns,double* essAlpha,
01023           double* essGamma, MPI_Comm comm,
01024           FEI_OSTREAM* dbgOut)
01025 {
01026   delete essBCEqns_;
01027   essBCEqns_ = new EqnBuffer();
01028 
01029   EqnBuffer* sendEssEqns = new EqnBuffer();
01030   ProcEqns* essSendProcEqns = new ProcEqns();
01031 
01032   std::vector<fei::CSVec*>& _sendEqns = sendEqns_->eqns();
01033   std::vector<int>& _sendEqnNumbers = sendEqns_->eqnNumbers();
01034   int _numSendEqns = sendEqns_->getNumEqns();
01035 
01036   //check to see if any of the essEqns are in the _sendEqns indices.
01037   //the ones that are, will need to be sent to other processors.
01038 
01039   if (dbgOut != NULL) {
01040     FEI_OSTREAM& os = *dbgOut;
01041     os << "#ereb: num-remote-rows: " << _numSendEqns
01042        << ", numEssEqns: " << numEssEqns << FEI_ENDL;
01043 //    for(int ii=0; ii<numEssEqns; ++ii) {
01044 //      os << "#ereb, essEqns["<<ii<<"]: "<<essEqns[ii]<<FEI_ENDL;
01045 //    }
01046     os << "#ereb sendEqns_:"<<FEI_ENDL;
01047     os << *sendEqns_<<FEI_ENDL;
01048   }
01049 
01050   bool accumulate = false;
01051 
01052   std::vector<int> offsets(numEssEqns);
01053   int* offsetsPtr = numEssEqns>0 ? &offsets[0] : NULL;
01054 
01055   for(int j=0; j<_numSendEqns; j++) {
01056 
01057     std::vector<int>& indices = _sendEqns[j]->indices();
01058 
01059     fei::binarySearch(numEssEqns, essEqns, offsetsPtr,
01060         &indices[0], indices.size());
01061 
01062     int sendEqn_j = _sendEqnNumbers[j];
01063 
01064     int proc = getSendProcNumber(sendEqn_j);
01065 
01066     const int* sendEqnsPtr_j = &indices[0];
01067 
01068     if (dbgOut != NULL) {
01069       FEI_OSTREAM& os = *dbgOut;
01070       os << "#ereb sendeqns["<<j<<"].length: "
01071          <<_sendEqns[j]->size()<<", numEssEqns: " << numEssEqns<<FEI_ENDL;
01072     }
01073 
01074     for(int i=0; i<numEssEqns; i++) {
01075 
01076       int index = offsetsPtr[i];
01077 
01078       if (index < 0) continue;
01079 
01080       essSendProcEqns->addEqn(sendEqn_j, proc);
01081 
01082       double coef = essGamma[i]/essAlpha[i];
01083       int essEqns_i = essEqns[i];
01084 
01085       int* essEqns_i_ptr = &essEqns_i;
01086 
01087       sendEssEqns->addEqn(sendEqn_j, &coef,
01088         essEqns_i_ptr, 1, accumulate);
01089 
01090       for(size_t k=0; k<_sendEqns[j]->size(); k++) {
01091 
01092   int row = sendEqnsPtr_j[k];
01093 
01094   essBCEqns_->addEqn(row, &coef, essEqns_i_ptr, 1, accumulate);
01095       }
01096     }
01097   }
01098 
01099   if (dbgOut != NULL) {
01100     FEI_OSTREAM& os = *dbgOut;
01101     os << "#ereb sendEssEqns:"<<FEI_ENDL;
01102     os << *sendEssEqns<<FEI_ENDL;
01103   }
01104 
01105   ProcEqns* essRecvProcEqns = new ProcEqns();
01106 
01107   CHK_ERR( mirrorProcEqns(*essSendProcEqns, *essRecvProcEqns) );
01108 
01109   std::vector<int>& eqnNumbers = sendEssEqns->eqnNumbers();
01110   std::vector<fei::CSVec*>& sendEssEqnsVec = sendEssEqns->eqns();
01111   std::vector<int> eqnLengths(eqnNumbers.size());
01112   for(size_t i=0; i<eqnNumbers.size(); ++i) {
01113     eqnLengths[i] = sendEssEqnsVec[i]->size();
01114   }
01115 
01116   if (eqnNumbers.size() > 0) {
01117     essSendProcEqns->setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
01118                                        eqnNumbers.size());
01119   }
01120 
01121   CHK_ERR( mirrorProcEqnLengths(*essSendProcEqns, *essRecvProcEqns) );
01122 
01123 
01124   CHK_ERR( exchangeEqnBuffers(comm, essSendProcEqns, sendEssEqns,
01125             essRecvProcEqns, essBCEqns_, accumulate) );
01126 
01127   if (dbgOut != NULL) {
01128     FEI_OSTREAM& os = *dbgOut;
01129     os << "essBCEqns:"<<FEI_ENDL;
01130     os << *essBCEqns_ << FEI_ENDL;
01131   }
01132 
01133   delete sendEssEqns;
01134   delete essSendProcEqns;
01135   delete essRecvProcEqns;
01136 
01137   return(0);
01138 }
01139 
01140 //------------------------------------------------------------------------------
01141 int EqnCommMgr::exchangePtToBlkInfo(snl_fei::PointBlockMap& blkEqnMapper)
01142 {
01143   std::set<int> sendIndices;
01144   std::vector<fei::CSVec*>& sendeqns = sendEqns_->eqns();
01145   for(size_t i=0; i<sendeqns.size(); ++i) {
01146     std::vector<int>& indices = sendeqns[i]->indices();
01147     int len = indices.size();
01148     if (len < 1) continue;
01149     int* indicesPtr = &indices[0];
01150     for(int j=0; j<len; ++j) {
01151       sendIndices.insert(indicesPtr[j]);
01152     }
01153   }
01154 
01155   std::set<int> recvIndices;
01156   std::vector<fei::CSVec*>& recveqns = recvEqns_->eqns();
01157   for(size_t i=0; i<recveqns.size(); ++i) {
01158     std::vector<int>& indices = recveqns[i]->indices();
01159     int len = indices.size();
01160     if (len < 1) continue;
01161     int* indicesPtr = &indices[0];
01162     for(int j=0; j<len; ++j) {
01163       recvIndices.insert(indicesPtr[j]);
01164     }
01165   }
01166 
01167   std::map<int,int>* ptEqns =  blkEqnMapper.getPtEqns();
01168   size_t numPtEqns = ptEqns->size();
01169 
01170   std::map<int,int>::const_iterator
01171     pteq = ptEqns->begin(),
01172     pteq_end = ptEqns->end();
01173 
01174   std::vector<int> ptBlkInfo(numPtEqns*3);
01175   int* ptBlkInfoPtr = &ptBlkInfo[0];
01176 
01177   int offset = 0;
01178   for(; pteq!=pteq_end; ++pteq) {
01179     int ptEqn = (*pteq).first;
01180     if (sendIndices.find(ptEqn) == sendIndices.end()) continue;
01181 
01182     int blkEqn = blkEqnMapper.eqnToBlkEqn(ptEqn);
01183     int blkSize = blkEqnMapper.getBlkEqnSize(blkEqn);
01184 
01185     ptBlkInfoPtr[offset++] = ptEqn;
01186     ptBlkInfoPtr[offset++] = blkEqn;
01187     ptBlkInfoPtr[offset++] = blkSize;
01188   }
01189 
01190   ptBlkInfo.resize(offset);
01191 
01192   size_t numRecvProcs = recvProcEqns_->getNumProcs();
01193   std::vector<int>& recvProcs = recvProcEqns_->procsPtr();
01194   size_t numSendProcs = sendProcEqns_->getNumProcs();
01195   std::vector<int>& sendProcs = sendProcEqns_->procsPtr();
01196 
01197   std::vector<std::vector<int>* > recvData(numRecvProcs);
01198   for(unsigned i=0; i<numRecvProcs; ++i) {
01199     recvData[i] = new std::vector<int>;
01200   }
01201 
01202   std::vector<std::vector<int>* > sendData(numSendProcs);
01203   for(unsigned i=0; i<numSendProcs; ++i) {
01204     sendData[i] = &ptBlkInfo;
01205   }
01206 
01207   CHK_ERR( fei::exchangeData(comm_, sendProcs, sendData,
01208             recvProcs, false, recvData) );
01209 
01210   for(unsigned i=0; i<numRecvProcs; ++i) {
01211     size_t len = recvData[i]->size()/3;
01212     int* dataPtr = &(*(recvData[i]))[0];
01213 
01214     offset = 0;
01215     for(unsigned eq=0; eq<len; ++eq) {
01216       int ptEqn = dataPtr[offset++];
01217 
01218       if (recvIndices.find(ptEqn) == recvIndices.end()) {
01219         offset += 2; continue;
01220       }
01221 
01222       int blkEqn = dataPtr[offset++];
01223       int blkSize = dataPtr[offset++];
01224 
01225       blkEqnMapper.setEqn(ptEqn, blkEqn, blkSize);
01226     }
01227     delete recvData[i];
01228   }
01229 
01230   return(0);
01231 }
01232 
01233 //------------------------------------------------------------------------------
01234 int EqnCommMgr::addRemoteEqns(fei::CSRMat& mat, bool onlyIndices)
01235 {
01236   std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
01237   std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
01238   std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
01239   std::vector<double>& pckCoefs = mat.getPackedCoefs();
01240 
01241   for(size_t i=0; i<rowNumbers.size(); ++i) {
01242     int row = rowNumbers[i];
01243     int offset = rowOffsets[i];
01244     int rowlen = rowOffsets[i+1]-offset;
01245     int* indices = &pckColInds[offset];
01246     double* coefs = &pckCoefs[offset];
01247 
01248     int proc = getSendProcNumber(row);
01249     if (proc == localProc_ || proc < 0) continue;
01250 
01251     if (onlyIndices) {
01252       addRemoteIndices(row, proc, indices, rowlen);
01253     }
01254     else {
01255       CHK_ERR( addRemoteEqn(row, proc, coefs, indices, rowlen) );
01256     }
01257   }
01258 
01259   return(0);
01260 }
01261 
01262 //------------------------------------------------------------------------------
01263 int EqnCommMgr::addRemoteRHS(fei::CSVec& vec, int rhsIndex)
01264 {
01265   std::vector<int>& indices = vec.indices();
01266   std::vector<double>& coefs = vec.coefs();
01267 
01268   for(size_t i=0; i<indices.size(); i++) {
01269     int proc = getSendProcNumber(indices[i]);
01270 
01271     if (proc == localProc_ || proc < 0) continue;
01272 
01273     CHK_ERR( addRemoteRHS(indices[i], proc, rhsIndex, coefs[i]) );
01274   }
01275 
01276   return(0);
01277 }
01278 
01279 //------------------------------------------------------------------------------
01280 int EqnCommMgr::getSendProcNumber(int eqn)
01281 {
01282   size_t numSendProcs = sendProcEqns_->getNumProcs();
01283   std::vector<int>& sendProcs = sendProcEqns_->procsPtr();
01284   std::vector<std::vector<int>*>& sendProcEqnNumbers = 
01285     sendProcEqns_->procEqnNumbersPtr();
01286 
01287   for(unsigned i=0; i<numSendProcs; i++) {
01288     if (std::binary_search(sendProcEqnNumbers[i]->begin(),
01289                            sendProcEqnNumbers[i]->end(), eqn)) {
01290 
01291       return(sendProcs[i]);
01292     }
01293   }
01294 
01295   return(-1);
01296 }
01297 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends
Generated on Wed Apr 13 10:08:23 2011 for FEI by  doxygen 1.6.3