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