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