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