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