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