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