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