FEI Version of the Day
fei_Matrix_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 #include <fei_macros.hpp>
00009 
00010 #include <fei_utils.hpp>
00011 
00012 #include <fei_ParameterSet.hpp>
00013 #include <fei_LogManager.hpp>
00014 #include <fei_Matrix_core.hpp>
00015 #include <fei_CSRMat.hpp>
00016 #include <fei_TemplateUtils.hpp>
00017 #include <fei_CommUtils.hpp>
00018 #include <fei_impl_utils.hpp>
00019 
00020 #include <fei_VectorSpace.hpp>
00021 #include <snl_fei_PointBlockMap.hpp>
00022 #include <fei_MatrixGraph.hpp>
00023 #include <snl_fei_Utils.hpp>
00024 
00025 #include <fei_MatrixTraits.hpp>
00026 #include <fei_MatrixTraits_FillableMat.hpp>
00027 
00028 #undef fei_file
00029 #define fei_file "fei_Matrix_core.cpp"
00030 #include <fei_ErrMacros.hpp>
00031 
00032 fei::Matrix_core::Matrix_core(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00033                               int numLocalEqns)
00034   : name_(),
00035     work_indices_(),
00036     work_indices2_(),
00037     work_ints_(),
00038     work_data1D_(),
00039     work_data2D_(),
00040     eqnComm_(),
00041     rhsVector_(),
00042     comm_(matrixGraph->getRowSpace()->getCommunicator()),
00043     localProc_(0),
00044     numProcs_(1),
00045     vecSpace_(),
00046     matrixGraph_(matrixGraph),
00047     remotelyOwned_(),
00048     remotelyOwned_last_requested_(NULL),
00049     sendProcs_(),
00050     recvProcs_(),
00051     sendRecvProcsNeedUpdated_(true),
00052     proc_last_requested_(-1),
00053     haveBlockMatrix_(false),
00054     haveFEMatrix_(false),
00055     globalOffsets_(),
00056     firstLocalOffset_(0),
00057     lastLocalOffset_(0)
00058 {
00059   if (matrixGraph.get() == NULL) {
00060     throw std::runtime_error("fei::Matrix_core constructed with NULL fei::MatrixGraph");
00061   }
00062 
00063   eqnComm_.reset(new fei::EqnComm(comm_, numLocalEqns));
00064 
00065   localProc_ = fei::localProc(comm_);
00066   numProcs_ = fei::numProcs(comm_);
00067 
00068   vecSpace_ = matrixGraph->getRowSpace();
00069 
00070   setName("dbg");
00071 
00072   globalOffsets_ = eqnComm_->getGlobalOffsets();
00073 
00074   firstLocalOffset_ = globalOffsets_[localProc_];
00075   lastLocalOffset_ = globalOffsets_[localProc_+1]-1;
00076 }
00077 
00078 std::map<int,fei::FillableMat*>&
00079 fei::Matrix_core::getRemotelyOwnedMatrices()
00080 {
00081   return remotelyOwned_;
00082 }
00083 
00084 void
00085 fei::Matrix_core::putScalar_remotelyOwned(double scalar)
00086 {
00087   std::map<int,FillableMat*>::iterator
00088     it = remotelyOwned_.begin(),
00089     it_end = remotelyOwned_.end();
00090   for(; it!=it_end; ++it) {
00091     fei::MatrixTraits<FillableMat>::setValues(it->second, scalar);
00092   }
00093 }
00094 
00095 void
00096 fei::Matrix_core::setEqnComm(fei::SharedPtr<fei::EqnComm> eqnComm)
00097 {
00098   eqnComm_ = eqnComm;
00099   globalOffsets_ = eqnComm_->getGlobalOffsets();
00100   firstLocalOffset_ = globalOffsets_[localProc_];
00101   lastLocalOffset_ = globalOffsets_[localProc_+1]-1;
00102 }
00103 
00104 fei::Matrix_core::~Matrix_core()
00105 {
00106   std::map<int,FillableMat*>::iterator
00107     it = remotelyOwned_.begin(),
00108     it_end = remotelyOwned_.end();
00109   for(; it!=it_end; ++it) {
00110     delete it->second;
00111   }
00112 }
00113 
00114 void
00115 fei::Matrix_core::parameters(const fei::ParameterSet& paramset)
00116 {
00117   const fei::Param* param = paramset.get("name");
00118   fei::Param::ParamType ptype = param != NULL ?
00119     param->getType() : fei::Param::BAD_TYPE;
00120   if (ptype == fei::Param::STRING) {
00121     setName(param->getStringValue().c_str());
00122   }
00123 
00124   param = paramset.get("FEI_OUTPUT_LEVEL");
00125   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00126   if (ptype == fei::Param::STRING) {
00127     fei::LogManager& log_manager = fei::LogManager::getLogManager();
00128     log_manager.setOutputLevel(param->getStringValue().c_str());
00129     setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
00130   }
00131 
00132 }
00133 
00134 void fei::Matrix_core::setName(const char* name)
00135 {
00136   if (name == NULL) return;
00137 
00138   if (name_ == name) return;
00139 
00140   name_ = name;
00141 }
00142 
00143 int fei::Matrix_core::getOwnerProc(int globalEqn) const
00144 {
00145   int len = globalOffsets_.size();
00146   if (globalEqn > globalOffsets_[len-1]) return(-1);
00147 
00148   for(int p=len-2; p>=0; --p) {
00149     if (globalEqn >= globalOffsets_[p]) {
00150       return(p);
00151     }
00152   }
00153 
00154   return(-1);
00155 }
00156 
00157 void fei::Matrix_core::setRHS(fei::SharedPtr<fei::Vector> rhsvector)
00158 {
00159   rhsVector_ = rhsvector;
00160 }
00161 
00162 int fei::Matrix_core::gatherFromOverlap(bool accumulate)
00163 {
00164   if (numProcs() == 1) return(0);
00165 
00166 #ifndef FEI_SER
00167   //this function gathers shared-but-not-owned data onto the
00168   //owning processors.
00169 
00170   if (sendRecvProcsNeedUpdated_) {
00171     //iterate the remotelyOwned_ map and create a list of processors
00172     //that we will be sending data to. (processors which own matrix rows
00173     //that we share.)
00174     sendProcs_.clear();
00175     std::map<int,FillableMat*>::iterator
00176       it = remotelyOwned_.begin(),
00177       it_end = remotelyOwned_.end();
00178     for(; it!=it_end; ++it) {
00179       if (it->first == localProc_) continue;
00180       if (it->second != NULL) {
00181         if (it->second->getNumRows() == 0) {
00182           continue;
00183         }
00184   
00185         sendProcs_.push_back(it->first);
00186       }
00187     }
00188     sendRecvProcsNeedUpdated_ = false;
00189 
00190     //now we can find out which procs we'll be receiving from.
00191     recvProcs_.clear();
00192     fei::mirrorProcs(comm_, sendProcs_, recvProcs_);
00193   }
00194 
00195   //next we'll declare arrays to receive into.
00196   std::vector<std::vector<char> > recv_chars(recvProcs_.size());
00197 
00198   //...and an array for the sizes of the recv buffers:
00199   std::vector<int> recv_sizes(recvProcs_.size());
00200   std::vector<MPI_Request> mpiReqs(recvProcs_.size());
00201   std::vector<MPI_Request> mpiReqs2(recvProcs_.size());
00202   std::vector<MPI_Status> mpiStatuses(recvProcs_.size());
00203 
00204   int tag1 = 11111;
00205 
00206   for(size_t i=0; i<recvProcs_.size(); ++i) {
00207     MPI_Irecv(&recv_sizes[i], 1, MPI_INT, recvProcs_[i],
00208               tag1, comm_, &mpiReqs[i]);
00209   }
00210 
00211   //now we'll pack our to-be-sent data into buffers, and send the
00212   //sizes to the receiving procs:
00213   std::vector<std::vector<char> > send_chars(sendProcs_.size());
00214 
00215   for(size_t i=0; i<sendProcs_.size(); ++i) {
00216     int proc = sendProcs_[i];
00217 
00218     fei::impl_utils::pack_FillableMat(*(remotelyOwned_[proc]),
00219                                       send_chars[i]);
00220 
00221     int bsize = send_chars[i].size();
00222 
00223     MPI_Send(&bsize, 1, MPI_INT, proc, tag1, comm_);
00224 
00225     remotelyOwned_[proc]->clear();//setValues(0.0);
00226   }
00227 
00228   int numRecvProcs = recvProcs_.size();
00229 
00230   //now resize our recv buffers, and post the recvs.
00231   for(size_t i=0; i<recvProcs_.size(); ++i) {
00232     int index;
00233     MPI_Status status;
00234     MPI_Waitany(numRecvProcs, &mpiReqs[0], &index, &status);
00235 
00236     int bsize = recv_sizes[index];
00237 
00238     recv_chars[index].resize(bsize);
00239 
00240     MPI_Irecv(&(recv_chars[index][0]), bsize, MPI_CHAR, recvProcs_[index],
00241               tag1, comm_, &mpiReqs2[index]);
00242   }
00243 
00244   //now send our packed buffers.
00245   for(size_t i=0; i<sendProcs_.size(); ++i) {
00246     int proc = sendProcs_[i];
00247 
00248     MPI_Send(&(send_chars[i][0]), send_chars[i].size(), MPI_CHAR, proc, tag1, comm_);
00249   }
00250 
00251   for(size_t ir=0; ir<recvProcs_.size(); ++ir) {
00252     int index;
00253     MPI_Status status;
00254     MPI_Waitany(numRecvProcs, &mpiReqs2[0], &index, &status);
00255   }
00256 
00257   //and finally, unpack and store the received buffers.
00258   CSRMat recvMat;
00259   for(size_t ir=0; ir<recvProcs_.size(); ++ir) {
00260     fei::impl_utils::unpack_CSRMat(recv_chars[ir], recvMat);
00261 
00262     int nrows = recvMat.getNumRows();
00263 
00264     for(int i=0; i<nrows; ++i) {
00265       int row = recvMat.getGraph().rowNumbers[i];
00266       int rowOffset = recvMat.getGraph().rowOffsets[i];
00267       int rowLen = recvMat.getGraph().rowOffsets[i+1] - rowOffset;
00268 
00269       int* indices = &(recvMat.getGraph().packedColumnIndices[rowOffset]);
00270       double* vals = &(recvMat.getPackedCoefs()[rowOffset]);
00271       int err = 0;
00272       if (haveBlockMatrix()) {
00273         err = giveToBlockMatrix(1, &row, rowLen, indices, &vals, accumulate);
00274       }
00275       else {
00276         err = giveToUnderlyingMatrix(1, &row, rowLen, indices,
00277                                      &vals, accumulate, FEI_DENSE_ROW);
00278       }
00279       if (err != 0) {
00280         FEI_COUT << "fei::Matrix_core::gatherFromOverlap ERROR storing "
00281          << "values for row=" << row<<", recv'd from proc " << recvProcs_[i]
00282          << FEI_ENDL;
00283         
00284         return(err);
00285       }
00286     }
00287   }
00288 #endif
00289 
00290   return(0);
00291 }
00292 
00293 void
00294 fei::Matrix_core::copyTransposeToWorkArrays(int numRows, int numCols,
00295                                             const double*const* values,
00296                                             std::vector<double>& work_1D,
00297                                             std::vector<const double*>& work_2D)
00298 {
00299   //copy the transpose of 'values' into work-arrays.
00300 
00301   int arrayLen = numCols*numRows;
00302   if ((int)work_1D.size() != arrayLen) {
00303     work_1D.resize(arrayLen);
00304   }
00305   if ((int)work_2D.size() != numCols) {
00306     work_2D.resize(numCols);
00307   }
00308 
00309   const double** dataPtr = &work_2D[0];
00310   double* data1DPtr = &work_1D[0];
00311 
00312   for(int i=0; i<numCols; ++i) {
00313     dataPtr[i] = data1DPtr;
00314 
00315     for(int j=0; j<numRows; ++j) {
00316       data1DPtr[j] = values[j][i];
00317     }
00318 
00319     data1DPtr += numRows;
00320   }
00321 }
00322 
00323 
00324 int fei::Matrix_core::convertPtToBlk(int numRows,
00325                const int* rows,
00326                int numCols,
00327                const int* cols,
00328                int* blkRows,
00329                int* blkRowOffsets,
00330                int* blkCols,
00331                int* blkColOffsets)
00332 {
00333   snl_fei::PointBlockMap* pointBlockMap = vecSpace_->getPointBlockMap();
00334 
00335   int i;
00336   for(i=0; i<numRows; ++i) {
00337     if (pointBlockMap->getPtEqnInfo(rows[i], blkRows[i], blkRowOffsets[i])!=0){
00338       ERReturn(-1);
00339     }
00340   }
00341   for(i=0; i<numCols; ++i) {
00342     if(pointBlockMap->getPtEqnInfo(cols[i], blkCols[i], blkColOffsets[i])!=0){
00343       ERReturn(-1);
00344     }
00345   }
00346 
00347   return(0);
00348 }
00349 
00350 int fei::Matrix_core::copyPointRowsToBlockRow(int numPtRows,
00351             int numPtCols,
00352             const double*const* ptValues,
00353             int numBlkCols,
00354             const int* blkColDims,
00355             double** blkValues)
00356 {
00357   int ptColOffset = 0;
00358   for(int blki=0; blki<numBlkCols; ++blki) {
00359     int blkvalueOffset = 0;
00360     double* blkvalues_i = blkValues[blki];
00361     for(int j=0; j<blkColDims[blki]; ++j) {
00362       int loc = ptColOffset+j;
00363       for(int i=0; i<numPtRows; ++i) {
00364         blkvalues_i[blkvalueOffset++] = ptValues[i][loc];
00365       }
00366     }
00367     ptColOffset += blkColDims[blki];
00368   }
00369 
00370   return(0);
00371 }
00372 
00373 void fei::Matrix_core::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00374 {
00375   matrixGraph_ = matrixGraph;
00376   vecSpace_ = matrixGraph->getRowSpace();
00377 }
00378 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends