FEI Version of the Day
fei_Vector_core.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_fstream.hpp"
00010 
00011 #include "fei_Vector_core.hpp"
00012 #include "fei_VectorSpace.hpp"
00013 #include "fei_CSVec.hpp"
00014 #include "snl_fei_RecordCollection.hpp"
00015 #include "fei_TemplateUtils.hpp"
00016 
00017 #undef fei_file
00018 #define fei_file "fei_Vector_core.cpp"
00019 
00020 #include "fei_ErrMacros.hpp"
00021 
00022 fei::Vector_core::Vector_core(fei::SharedPtr<fei::VectorSpace> vecSpace,
00023                               int numLocalEqns)
00024   : eqnComm_(),
00025     vecSpace_(vecSpace),
00026     comm_(vecSpace->getCommunicator()),
00027     firstLocalOffset_(0),
00028     lastLocalOffset_(0),
00029     numLocal_(0),
00030     work_indices_(),
00031     work_indices2_(),
00032     haveFEVector_(false),
00033     remotelyOwned_(),
00034     overlapAlreadySet_(false),
00035     dbgprefix_("Vcore: ")
00036 {
00037   eqnComm_.reset(new fei::EqnComm(comm_,numLocalEqns));
00038   remotelyOwned_.resize(fei::numProcs(comm_));
00039   for(unsigned i=0; i<remotelyOwned_.size(); ++i) {
00040     remotelyOwned_[i] = new CSVec;
00041   }
00042 
00043   const std::vector<int>& offsets = eqnComm_->getGlobalOffsets();
00044   firstLocalOffset_ = offsets[fei::localProc(comm_)];
00045   lastLocalOffset_ = offsets[fei::localProc(comm_)+1] - 1;
00046   numLocal_ = lastLocalOffset_ - firstLocalOffset_ + 1;
00047 
00048   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00049     FEI_OSTREAM& os = *output_stream_;
00050     os<<dbgprefix_<<" ctor firstLocal="<<firstLocalOffset_<<", lastLocal="
00051      <<lastLocalOffset_<<FEI_ENDL;
00052   }
00053 }
00054 
00055 fei::Vector_core::~Vector_core()
00056 {
00057   for(unsigned i=0; i<remotelyOwned_.size(); ++i) {
00058     delete remotelyOwned_[i];
00059   }
00060 }
00061 
00062 void fei::Vector_core::setOverlap(int numRemoteEqns,
00063                                   const int* remoteEqns)
00064 {
00065   if (numRemoteEqns == 0 && remoteEqns == NULL) {
00066     if (overlapAlreadySet_) return;
00067   }
00068 
00069   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00070     FEI_OSTREAM& os = *output_stream_;
00071     os << dbgprefix_<<"setOverlap"<<FEI_ENDL;
00072   }
00073 
00074   int local_proc = fei::localProc(comm_);
00075 
00076   if (numRemoteEqns != 0 && remoteEqns != NULL) {
00077     for(int i=0; i<numRemoteEqns; ++i) {
00078       int proc = eqnComm_->getOwnerProc(remoteEqns[i]);
00079       if (proc == local_proc) continue;
00080 
00081       fei::add_entry(*remotelyOwned_[proc], remoteEqns[i], 0.0);
00082     }
00083   }
00084   else {
00085     std::vector<int> eqns;
00086     vecSpace_->getIndices_SharedAndOwned(eqns);
00087 
00088     for(size_t i=0; i<eqns.size(); ++i) {
00089       int proc = eqnComm_->getOwnerProc(eqns[i]);
00090       if (proc == local_proc) continue;
00091 
00092       fei::add_entry(*remotelyOwned_[proc], eqns[i], 0.0);
00093     }
00094   }
00095 
00096   overlapAlreadySet_ = true;
00097 }
00098 
00099 int fei::Vector_core::scatterToOverlap()
00100 {
00101   if (fei::numProcs(comm_) == 1 || haveFEVector()) {
00102     return(0);
00103   }
00104 
00105 #ifndef FEI_SER
00106   if (!overlapAlreadySet_) {
00107     setOverlap();
00108   }
00109 
00110   //...and now the overlap is whatever is in our remotelyOwned_ vectors.
00111 
00112   //first find out which procs we'll be receiving from.
00113   std::vector<int> recvProcs;
00114   for(unsigned i=0; i<remotelyOwned_.size(); ++i) {
00115     if ((int)i == fei::localProc(comm_)) continue;
00116     if (remotelyOwned_[i]->size() == 0) continue;
00117 
00118     recvProcs.push_back((int)i);
00119   }
00120 
00121   //find out the send-procs.
00122   std::vector<int> sendProcs;
00123   fei::mirrorProcs(comm_, recvProcs, sendProcs);
00124 
00125   //declare arrays to send from, and corresponding sizes
00126   std::vector<std::vector<int> > send_ints(sendProcs.size());
00127   std::vector<std::vector<double> > send_doubles(sendProcs.size());
00128   std::vector<int> send_sizes(sendProcs.size());
00129 
00130   std::vector<MPI_Request> mpiReqs(sendProcs.size()+recvProcs.size());
00131   std::vector<MPI_Status> mpiStatuses(sendProcs.size()+recvProcs.size());
00132   int tag1 = 11111;
00133   int tag2 = 11112;
00134 
00135   //first, the procs we're going to send to, have to let us know
00136   //how much data we're supposed to send. So we have to receive
00137   //sizes and then indices from the "send"-procs.
00138   for(unsigned i=0; i<sendProcs.size(); ++i) {
00139     MPI_Irecv(&send_sizes[i], 1, MPI_INT, sendProcs[i],
00140               tag1, comm_, &mpiReqs[i]);
00141   }
00142 
00143   //now we'll send the sizes of our remotely-owned data to the
00144   //procs that we will be receiving the data from, and also the
00145   //indices that we want to receive.
00146   for(unsigned i=0; i<recvProcs.size(); ++i) {
00147     int proc = recvProcs[i];
00148 
00149     int size = remotelyOwned_[proc]->size();
00150     MPI_Send(&size, 1, MPI_INT, proc, tag1, comm_);
00151   }
00152  
00153   MPI_Waitall(sendProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
00154 
00155   //now resize our send_ints and send_doubles arrays, and post the recvs
00156   //for indices that we're supposed to pack.
00157   for(unsigned i=0; i<sendProcs.size(); ++i) {
00158     int proc = sendProcs[i];
00159     int size = send_sizes[i];
00160     send_ints[i].resize(size);
00161     MPI_Irecv(&(send_ints[i][0]), size, MPI_INT, proc, tag1,
00162               comm_, &mpiReqs[i]);
00163     send_doubles[i].resize(size);
00164   }
00165 
00166   //now send the indices that we want to receive data for.
00167   for(unsigned i=0; i<recvProcs.size(); ++i) {
00168     int proc = recvProcs[i];
00169     int size = remotelyOwned_[proc]->size();
00170     int* indices = &(remotelyOwned_[proc]->indices())[0];
00171     MPI_Send(indices, size, MPI_INT, proc, tag1, comm_);
00172   }
00173 
00174   MPI_Waitall(sendProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
00175 
00176   //now post our recvs.
00177   for(unsigned i=0; i<recvProcs.size(); ++i) {
00178     int proc = recvProcs[i];
00179     int size = remotelyOwned_[proc]->size();
00180     double* coefs = &(remotelyOwned_[proc]->coefs())[0];
00181     MPI_Irecv(coefs, size, MPI_DOUBLE, proc, tag2, comm_, &mpiReqs[i]);
00182   }
00183 
00184   //now pack and send the coefs that the other procs need from us.
00185   for(unsigned i=0; i<sendProcs.size(); ++i) {
00186     int proc = sendProcs[i];
00187 
00188     int num = send_sizes[i];
00189     int err = copyOutOfUnderlyingVector(num, &(send_ints[i][0]),
00190                                         &(send_doubles[i][0]), 0);
00191     if (err != 0) {
00192       FEI_COUT << "fei::Vector_core::scatterToOverlap ERROR getting data to send."<<FEI_ENDL;
00193       return(err);
00194     }
00195 
00196     MPI_Send(&(send_doubles[i][0]), num, MPI_DOUBLE, proc,
00197              tag2, comm_);
00198   }
00199 
00200   MPI_Waitall(recvProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
00201 
00202 #endif  //#ifndef FEI_SER
00203 
00204   return(0);
00205 }
00206 
00207 int fei::Vector_core::copyOut(int numValues,
00208           const int* indices,
00209           double* values,
00210           int vectorIndex) const
00211 {
00212   const std::vector<CSVec*>& remote = remotelyOwned();
00213 
00214   for(int i=0; i<numValues; ++i) {
00215     int ind = indices[i];
00216 
00217     int local = ind - firstLocalOffset_;
00218     if (local < 0 || local >= numLocal_) {
00219       if (ind < 0) {
00220   continue;
00221       }
00222 
00223       int proc = eqnComm_->getOwnerProc(ind);
00224 
00225       int insertPoint = -1;
00226       int idx = fei::binarySearch(ind, remote[proc]->indices(), insertPoint);
00227       if (idx < 0) {
00228   FEI_CERR << "fei::Vector_core::copyOut: proc " << fei::localProc(comm_)
00229        << ", index " << ind << " not in remotelyOwned_ vec object for proc "
00230             <<proc<<FEI_ENDL;
00231   ERReturn(-1);
00232       }
00233       else {
00234   values[i] = remote[proc]->coefs()[idx];
00235       }
00236     }
00237     else {
00238       CHK_ERR( copyOutOfUnderlyingVector(1, &ind, &(values[i]), vectorIndex) );
00239     }
00240   }
00241 
00242   return(0);
00243 }
00244 
00245 int fei::Vector_core::giveToVector(int numValues,
00246                const int* indices,
00247                const double* values,
00248                bool sumInto,
00249                int vectorIndex)
00250 {
00251   std::vector<CSVec*>& remote = remotelyOwned();
00252 
00253   for(int i=0; i<numValues; ++i) {
00254     int ind = indices[i];
00255     double val = values[i];
00256 
00257     if (ind < 0) {
00258 //      throw std::runtime_error("negative index not allowed");
00259       //preservation of existing behavior: certain Sierra scenarios
00260       //involve passing negative indices for positions that should be
00261       //ignored... so we'll continue rather than throwing.
00262       continue;
00263     }
00264     int proc = eqnComm_->getOwnerProc(ind);
00265     int local = ind - firstLocalOffset_;
00266     if (local < 0 || local >= numLocal_) {
00267       if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00268         FEI_OSTREAM& os = *output_stream_;
00269         os << dbgprefix_<<"giveToVector remote["<<proc<<"]("
00270          <<ind<<","<<val<<")"<<FEI_ENDL;
00271       }
00272 
00273       if (sumInto) {
00274   fei::add_entry( *remote[proc], ind, val);
00275       }
00276       else {
00277   fei::put_entry( *remote[proc], ind, val);
00278       }
00279     }
00280     else {
00281       int err = giveToUnderlyingVector(1, &ind, &val, sumInto, vectorIndex);
00282       if (err != 0) {
00283   FEI_CERR << "giveToVector sumIn ERROR, ind: " << ind
00284        << ", val: " << val << FEI_ENDL;
00285   ERReturn(-1);
00286       }
00287     }
00288   }
00289 
00290   return(0);
00291 }
00292 
00293 int fei::Vector_core::assembleFieldData(int fieldID,
00294               int idType,
00295               int numIDs,
00296               const int* IDs,
00297               const double* data,
00298               bool sumInto,
00299               int vectorIndex)
00300 {
00301   if (vecSpace_.get() == NULL) ERReturn(-1);
00302 
00303   int fieldSize = vecSpace_->getFieldSize(fieldID);
00304 
00305   work_indices_.resize(numIDs*fieldSize);
00306   int* indicesPtr = &work_indices_[0];
00307 
00308   CHK_ERR( vecSpace_->getGlobalIndices(numIDs, IDs, idType, fieldID,
00309           indicesPtr) );
00310 
00311   CHK_ERR( giveToVector(numIDs*fieldSize, indicesPtr, data, sumInto, vectorIndex) );
00312 
00313   return(0);
00314 }
00315 
00316 int fei::Vector_core::gatherFromOverlap(bool accumulate)
00317 {
00318   if (fei::numProcs(comm_) == 1 || haveFEVector()) {
00319     return(0);
00320   }
00321 
00322 #ifndef FEI_SER
00323   //first create the list of procs we'll be sending to.
00324   std::vector<int> sendProcs;
00325   for(unsigned i=0; i<remotelyOwned_.size(); ++i) {
00326     if ((int)i == fei::localProc(comm_)) continue;
00327     if (remotelyOwned_[i]->size() == 0) continue;
00328 
00329     sendProcs.push_back(i);
00330   }
00331 
00332   std::vector<int> recvProcs;
00333   fei::mirrorProcs(comm_, sendProcs, recvProcs);
00334 
00335   //declare arrays to hold the indices and coefs we'll be receiving.
00336   std::vector<std::vector<int> > recv_ints(recvProcs.size());
00337   std::vector<std::vector<double> > recv_doubles(recvProcs.size());
00338   std::vector<int> recv_sizes(recvProcs.size());
00339 
00340   std::vector<MPI_Request> mpiReqs(recvProcs.size()*2);
00341   std::vector<MPI_Status> mpiStatuses(recvProcs.size()*2);
00342   int tag1 = 11111;
00343   int tag2 = 11112;
00344 
00345   //post the recvs for the sizes.
00346   for(size_t i=0; i<recvProcs.size(); ++i) {
00347     int proc = recvProcs[i];
00348     MPI_Irecv(&recv_sizes[i], 1, MPI_INT, proc,
00349               tag1, comm_, &mpiReqs[i]);
00350   }
00351 
00352   //send the sizes of data we'll be sending.
00353   for(unsigned i=0; i<sendProcs.size(); ++i) {
00354     int proc = sendProcs[i];
00355     int size = remotelyOwned_[proc]->size();
00356     MPI_Send(&size, 1, MPI_INT, proc, tag1, comm_);
00357   }
00358 
00359   if (recvProcs.size() > 0) {
00360     MPI_Waitall(recvProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
00361   }
00362 
00363   //now post the recvs for the data.
00364   unsigned offset = 0;
00365   for(size_t i=0; i<recvProcs.size(); ++i) {
00366     int proc = recvProcs[i];
00367     int size = recv_sizes[i];
00368     std::vector<int>& recv_ints_i = recv_ints[i];
00369     std::vector<double>& recv_doubles_i = recv_doubles[i];
00370     recv_ints_i.resize(size);
00371     recv_doubles_i.resize(size);
00372     MPI_Irecv(&(recv_ints_i[0]), size, MPI_INT, proc,
00373               tag1, comm_, &mpiReqs[offset++]);
00374     MPI_Irecv(&(recv_doubles_i[0]), size, MPI_DOUBLE, proc,
00375               tag2, comm_, &mpiReqs[offset++]);
00376   }
00377 
00378   //now send the outgoing data.
00379   for(size_t i=0; i<sendProcs.size(); ++i) {
00380     int proc = sendProcs[i];
00381     int size = remotelyOwned_[proc]->size();
00382     int* indices = &(remotelyOwned_[proc]->indices())[0];
00383     MPI_Send(indices, size, MPI_INT, proc, tag1, comm_);
00384     double* coefs = &(remotelyOwned_[proc]->coefs())[0];
00385     MPI_Send(coefs, size, MPI_DOUBLE, proc, tag2, comm_);
00386 
00387     fei::set_values(*remotelyOwned_[proc], 0.0);
00388   }
00389 
00390   if (recvProcs.size() > 0) {
00391     MPI_Waitall(recvProcs.size()*2, &mpiReqs[0], &mpiStatuses[0]);
00392   }
00393 
00394   //now store the data we've received.
00395   for(size_t i=0; i<recvProcs.size(); ++i) {
00396     int num = recv_sizes[i];
00397     std::vector<int>& recv_ints_i = recv_ints[i];
00398     std::vector<double>& recv_doubles_i = recv_doubles[i];
00399     int err = giveToUnderlyingVector(num, &(recv_ints_i[0]),
00400                                      &(recv_doubles_i[0]), accumulate, 0);
00401     if (err != 0) {
00402       FEI_COUT << "fei::Vector_core::gatherFromOverlap ERROR storing recvd data" << FEI_ENDL;
00403       return(err);
00404     }
00405   }
00406 
00407 #endif  //#ifndef FEI_SER
00408 
00409   return(0);
00410 }
00411 
00412 int fei::Vector_core::copyOutFieldData(int fieldID,
00413              int idType,
00414              int numIDs,
00415              const int* IDs,
00416              double* data,
00417              int vectorIndex)
00418 {
00419   if (vecSpace_.get() == NULL) ERReturn(-1);
00420 
00421   int fieldSize = vecSpace_->getFieldSize(fieldID);
00422 
00423   if (haveFEVector_) {
00424     snl_fei::RecordCollection* collection = NULL;
00425     CHK_ERR( vecSpace_->getRecordCollection(idType, collection) );
00426     int nodeNumber;
00427     int dofOffset;
00428     int numInstances;
00429     int foffset;
00430     std::vector<int>& eqnNums = vecSpace_->getEqnNumbers();
00431     int* vspcEqnPtr = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
00432 
00433     int offset = 0;
00434     for(int i=0; i<numIDs; ++i) {
00435       fei::Record<int>* node = collection->getRecordWithID(IDs[i]);
00436       if (node == NULL) {
00437         ERReturn(-1);
00438       }
00439 
00440       nodeNumber = node->getNumber();
00441       int* eqnNumbers = vspcEqnPtr+node->getOffsetIntoEqnNumbers();
00442       node->getFieldMask()->getFieldEqnOffset(fieldID, foffset, numInstances);
00443       dofOffset = eqnNumbers[foffset] - eqnNumbers[0];
00444       for(int j=0; j<fieldSize; ++j) {
00445   CHK_ERR( copyOut_FE(nodeNumber, dofOffset+j, data[offset++]));
00446       }
00447     }
00448   }
00449   else {
00450     work_indices_.resize(numIDs*fieldSize*2);
00451     int* indicesPtr = &work_indices_[0];
00452 
00453     CHK_ERR( vecSpace_->getGlobalIndices(numIDs, IDs, idType,
00454            fieldID, indicesPtr) );
00455 
00456     CHK_ERR( copyOut(numIDs*fieldSize, indicesPtr, data) );
00457   }
00458 
00459   return(0);
00460 }
00461 
00462 int fei::Vector_core::writeToFile(const char* filename,
00463                                     bool matrixMarketFormat)
00464 {
00465   int numProcs = fei::numProcs(comm_);
00466   int localProc =fei::localProc(comm_);
00467 
00468   double coef;
00469 
00470   static char mmbanner[] = "%%MatrixMarket matrix array real general";
00471 
00472   for(int p=0; p<numProcs; ++p) {
00473     fei::Barrier(comm_);
00474     if (p != localProc) continue;
00475 
00476     FEI_OFSTREAM* outFile = NULL;
00477     if (p==0) {
00478       outFile = new FEI_OFSTREAM(filename, IOS_OUT);
00479       FEI_OFSTREAM& ofref = *outFile;
00480       if (matrixMarketFormat) {
00481         ofref << mmbanner << FEI_ENDL;
00482         ofref << eqnComm_->getGlobalOffsets()[numProcs] << " 1" << FEI_ENDL;
00483       }
00484       else {
00485         ofref << eqnComm_->getGlobalOffsets()[numProcs] << FEI_ENDL;
00486       }
00487     }
00488     else outFile = new FEI_OFSTREAM(filename, IOS_APP);
00489     FEI_OFSTREAM& ofref = *outFile;
00490     ofref.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00491     ofref.precision(13);
00492 
00493     for(int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
00494       CHK_ERR( copyOut(1, &i, &coef) );
00495       if (matrixMarketFormat) {
00496         ofref << " " << coef << FEI_ENDL;
00497       }
00498       else {
00499         ofref << i << " " << coef << FEI_ENDL;
00500       }
00501     }
00502     
00503     delete outFile;
00504   } 
00505     
00506   return(0);
00507 }
00508 
00509 int fei::Vector_core::writeToStream(FEI_OSTREAM& ostrm,
00510           bool matrixMarketFormat)
00511 {
00512   int numProcs = fei::numProcs(comm_);
00513   int local_proc =fei::localProc(comm_);
00514 
00515   double coef;
00516 
00517   static char mmbanner[] = "%%MatrixMarket matrix array real general";
00518 
00519   IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00520   ostrm.precision(13);
00521 
00522   for(int proc=0; proc<numProcs; ++proc) {
00523     fei::Barrier(comm_);
00524     if (proc != local_proc) continue;
00525 
00526     if (proc==0) {
00527       if (matrixMarketFormat) {
00528   ostrm << mmbanner << FEI_ENDL;
00529   ostrm << eqnComm_->getGlobalOffsets()[numProcs] << " 1" << FEI_ENDL;
00530       }
00531       else {
00532   ostrm << eqnComm_->getGlobalOffsets()[numProcs] << FEI_ENDL;
00533       }
00534     }
00535 
00536     for(int p=0; p<local_proc; ++p) {
00537       for(size_t ii=0; ii<remotelyOwned_[p]->size(); ++ii) {
00538   if (matrixMarketFormat) {
00539     ostrm << " " << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
00540   }
00541   else {
00542     ostrm << " " << remotelyOwned_[p]->indices()[ii] << " "
00543     << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
00544   }
00545       }
00546     }
00547 
00548     for(int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
00549       CHK_ERR( copyOut(1, &i, &coef) );
00550       if (matrixMarketFormat) {
00551   ostrm << " " << coef << FEI_ENDL;
00552       }
00553       else {
00554   ostrm << " " << i << " " << coef << FEI_ENDL;
00555       }
00556     }
00557 
00558     for(int p=local_proc+1; p<numProcs; ++p) {
00559       for(size_t ii=0; ii<remotelyOwned_[p]->size(); ++ii) {
00560   if (matrixMarketFormat) {
00561     ostrm << " " << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
00562   }
00563   else {
00564     ostrm << " " << remotelyOwned_[p]->indices()[ii] << " "
00565     << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
00566   }
00567       }
00568     }
00569   }
00570 
00571   ostrm.setf(oldf, IOS_FLOATFIELD);
00572 
00573   return(0);
00574 }
00575 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends