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