FEI Version of the Day
fei_MatrixReducer.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2007 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 
00009 #include <fei_ParameterSet.hpp>
00010 #include "fei_MatrixReducer.hpp"
00011 #include "fei_EqnComm.hpp"
00012 #include "fei_Matrix_core.hpp"
00013 #include "fei_sstream.hpp"
00014 #include "fei_fstream.hpp"
00015 #include <fei_CommUtils.hpp>
00016 
00017 namespace fei {
00018 
00019 MatrixReducer::MatrixReducer(fei::SharedPtr<fei::Reducer> reducer,
00020                              fei::SharedPtr<fei::Matrix> target)
00021  : reducer_(reducer),
00022    target_(target),
00023    globalAssembleCalled_(false),
00024    changedSinceMark_(false)
00025 {
00026   fei::SharedPtr<fei::VectorSpace> vspace =
00027     target->getMatrixGraph()->getRowSpace();
00028   MPI_Comm comm = vspace->getCommunicator();
00029   int numLocal = reducer_->getLocalReducedEqns().size();
00030   fei::SharedPtr<fei::EqnComm> eqnComm(new fei::EqnComm(comm, numLocal));
00031   fei::Matrix_core* target_core =
00032     dynamic_cast<fei::Matrix_core*>(target_.get());
00033     if (target_core == NULL) {
00034     throw std::runtime_error("fei::MatrixReducer ERROR, target matrix not dynamic_cast-able to fei::Matrix_core.");
00035   }
00036 
00037   target_core->setEqnComm(eqnComm);
00038 }
00039 
00040 MatrixReducer::~MatrixReducer()
00041 {
00042 }
00043 
00044 int
00045 MatrixReducer::parameters(const fei::ParameterSet& paramset)
00046 {
00047   return(target_->parameters(paramset));
00048 }
00049 
00050 void
00051 MatrixReducer::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00052 {
00053   target_->setMatrixGraph(matrixGraph);
00054 }
00055 
00056 int
00057 MatrixReducer::getGlobalNumRows() const
00058 {
00059   return(target_->getMatrixGraph()->getRowSpace()->getGlobalNumIndices());
00060 }
00061 
00062 int
00063 MatrixReducer::getLocalNumRows() const
00064 {
00065   return(target_->getMatrixGraph()->getRowSpace()->getNumIndices_Owned());
00066 }
00067 
00068 int MatrixReducer::putScalar(double scalar)
00069 { return(target_->putScalar(scalar)); }
00070 
00071 int
00072 MatrixReducer::getRowLength(int row, int& length) const
00073 {
00074   if (reducer_->isSlaveEqn(row)) {
00075     FEI_OSTRINGSTREAM osstr;
00076     osstr << "fei::MatrixReducer::getRowLength ERROR, row="<<row<<" is a slave eqn. You can't get a slave row from the reduced matrix.";
00077     throw std::runtime_error(osstr.str());
00078   }
00079 
00080   int reducedrow = reducer_->translateToReducedEqn(row);
00081   return(target_->getRowLength(reducedrow, length));
00082 }
00083 
00084 int
00085 MatrixReducer::copyOutRow(int row, int len, double* coefs, int* indices) const
00086 {
00087   if (reducer_->isSlaveEqn(row)) {
00088     FEI_OSTRINGSTREAM osstr;
00089     osstr << "fei::MatrixReducer::copyOutRow ERROR, requested row ("<<row
00090       <<") is a slave eqn. You can't get a slave row from the reduced matrix.";
00091     throw std::runtime_error(osstr.str());
00092   }
00093 
00094   int reducedrow = reducer_->translateToReducedEqn(row);
00095   int err = target_->copyOutRow(reducedrow, len, coefs, indices);
00096   for(int i=0; i<len; ++i) {
00097     indices[i] = reducer_->translateFromReducedEqn(indices[i]);
00098   }
00099   return(err);
00100 }
00101 
00102 int
00103 MatrixReducer::sumIn(int numRows, const int* rows,
00104                      int numCols, const int* cols,
00105                      const double* const* values,
00106                      int format)
00107 {
00108   int err = reducer_->addMatrixValues(numRows, rows, numCols, cols,
00109                                       values, true, *target_, format);
00110   return(err);
00111 }
00112 
00113 int
00114 MatrixReducer::copyIn(int numRows, const int* rows,
00115                       int numCols, const int* cols,
00116                       const double* const* values,
00117                       int format)
00118 {
00119   int err = reducer_->addMatrixValues(numRows, rows, numCols, cols,
00120                                       values, false, *target_, format);
00121   return(err);
00122 }
00123 
00124 int
00125 MatrixReducer::sumInFieldData(int fieldID,
00126                               int idType,
00127                               int rowID,
00128                               int colID,
00129                               const double* const* data,
00130                               int format)
00131 {
00132   fei::SharedPtr<fei::VectorSpace> rowSpace =
00133     target_->getMatrixGraph()->getRowSpace();
00134   fei::SharedPtr<fei::VectorSpace> colSpace =
00135     target_->getMatrixGraph()->getColSpace();
00136 
00137   unsigned fieldSize = rowSpace->getFieldSize(fieldID);
00138   std::vector<int> indices(fieldSize*2);
00139   int* rowIndices = &indices[0];
00140   int* colIndices = rowIndices+fieldSize;
00141 
00142   rowSpace->getGlobalIndices(1, &rowID, idType, fieldID, rowIndices);
00143   colSpace->getGlobalIndices(1, &colID, idType, fieldID, colIndices);
00144 
00145   if (format != FEI_DENSE_ROW) {
00146     throw std::runtime_error("MatrixReducer: bad format");
00147   }
00148 
00149   int err = reducer_->addMatrixValues(fieldSize, rowIndices,
00150                                       fieldSize, colIndices,
00151                                       data, true, *target_, format);
00152   return(err);
00153 }
00154 
00155 int
00156 MatrixReducer::sumInFieldData(int fieldID,
00157              int idType,
00158              int rowID,
00159              int colID,
00160              const double* data,
00161              int format)
00162 {
00163   fei::SharedPtr<fei::VectorSpace> rowSpace =
00164     target_->getMatrixGraph()->getRowSpace();
00165 
00166   unsigned fieldSize = rowSpace->getFieldSize(fieldID);
00167 
00168   std::vector<const double*> data_2d(fieldSize);
00169   for(unsigned i=0; i<fieldSize; ++i) {
00170     data_2d[i] = &data[i*fieldSize];
00171   }
00172 
00173   return(sumInFieldData(fieldID, idType, rowID, colID, &data_2d[0], format));
00174 }
00175 
00176 int
00177 MatrixReducer::sumIn(int blockID, int connectivityID,
00178           const double* const* values,
00179           int format)
00180 {
00181   fei::SharedPtr<fei::MatrixGraph> matGraph = getMatrixGraph();
00182   int numRowIndices, numColIndices, dummy;
00183   matGraph->getConnectivityNumIndices(blockID, numRowIndices, numColIndices);
00184 
00185   std::vector<int> indices(numRowIndices+numColIndices);
00186   int* rowIndices = &indices[0];
00187   int* colIndices = rowIndices+numRowIndices;
00188 
00189   matGraph->getConnectivityIndices(blockID, connectivityID,
00190                                    numRowIndices, rowIndices, dummy,
00191                                    numColIndices, colIndices, dummy);
00192 
00193   return(sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
00194                values, format));
00195 }
00196 
00197 int
00198 MatrixReducer::globalAssemble()
00199 {
00200   reducer_->assembleReducedMatrix(*target_);
00201   return(target_->globalAssemble());
00202 }
00203 
00204 int MatrixReducer::multiply(fei::Vector* x, fei::Vector* y)
00205 { return(target_->multiply(x, y)); }
00206 
00207 int MatrixReducer::gatherFromOverlap(bool accumulate)
00208 {
00209   reducer_->assembleReducedMatrix(*target_);
00210   return(target_->gatherFromOverlap(accumulate));
00211 }
00212 
00213 int MatrixReducer::writeToFile(const char* filename,
00214           bool matrixMarketFormat)
00215 {
00216   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
00217   std::vector<int>& localrows = reducer_->getLocalReducedEqns();
00218   int localNumRows = localrows.size();
00219 
00220   int globalNNZ = 0;
00221   int localNNZ = 0;
00222 
00223   for(int i=0; i<localNumRows; ++i) {
00224     int len;
00225     CHK_ERR( target_->getRowLength(localrows[i], len) );
00226     localNNZ += len;
00227   }
00228 
00229   MPI_Comm comm = getMatrixGraph()->getRowSpace()->getCommunicator();
00230 
00231   CHK_MPI( fei::GlobalSum(comm, localNNZ, globalNNZ) );
00232   int globalNumRows = 0;
00233   CHK_MPI( fei::GlobalSum(comm, localNumRows, globalNumRows) );
00234 
00235   int globalNumCols = globalNumRows;
00236 
00237   for(int p=0; p<fei::numProcs(comm); ++p) {
00238     fei::Barrier(comm);
00239     if (p != fei::localProc(comm)) continue;
00240 
00241     FEI_OFSTREAM* outFile = NULL;
00242     if (p==0) {
00243       outFile = new FEI_OFSTREAM(filename, IOS_OUT);
00244       FEI_OFSTREAM& ofs = *outFile;
00245       if (matrixMarketFormat) {
00246         ofs << mmbanner << FEI_ENDL;
00247         ofs <<globalNumRows<< " " <<globalNumCols<< " " <<globalNNZ<<FEI_ENDL;
00248       }
00249       else {
00250         ofs <<globalNumRows<< " " <<globalNumCols<<FEI_ENDL;
00251       }
00252     }
00253     else outFile = new FEI_OFSTREAM(filename, IOS_APP);
00254 
00255     outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00256     outFile->precision(13);
00257     FEI_OFSTREAM& ofs = *outFile;
00258 
00259     int rowLength;
00260     std::vector<int> work_indices;
00261     std::vector<double> work_data1D;
00262 
00263     for(int i=0; i<localNumRows; ++i) {
00264       int row = localrows[i];
00265       CHK_ERR( target_->getRowLength(row, rowLength) );
00266 
00267       work_indices.resize(rowLength);
00268       work_data1D.resize(rowLength);
00269 
00270       int* indPtr = &work_indices[0];
00271       double* coefPtr = &work_data1D[0];
00272 
00273       CHK_ERR( target_->copyOutRow(row, rowLength, coefPtr, indPtr) );
00274 
00275       for(int j=0; j<rowLength; ++j) {
00276         if (matrixMarketFormat) {
00277           ofs << row+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
00278         }
00279         else {
00280           ofs << row <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
00281         }
00282       }
00283     }
00284 
00285     delete outFile;
00286   }
00287 
00288   return(0);
00289 }
00290 
00291 int MatrixReducer::writeToStream(FEI_OSTREAM& ostrm,
00292             bool matrixMarketFormat)
00293 {
00294   return(target_->writeToStream(ostrm, matrixMarketFormat));
00295 }
00296 
00297 void MatrixReducer::markState()
00298 { target_->markState(); }
00299 
00300 bool MatrixReducer::changedSinceMark()
00301 { return(target_->changedSinceMark()); }
00302 
00303 }//namespace fei
00304 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends