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