PlayaEpetraMatrix.cpp

00001 /* @HEADER@ */
00002 //   
00003  /* @HEADER@ */
00004 
00005 #include "PlayaEpetraMatrix.hpp"
00006 #include "PlayaEpetraVector.hpp"
00007 #include "PlayaVectorSpaceDecl.hpp"  // changed from Impl
00008 #include "PlayaVectorDecl.hpp"
00009 #include "PlayaLinearOperatorDecl.hpp"  // changed from Impl
00010 #include "PlayaIfpackOperator.hpp"
00011 #include "PlayaPreconditioner.hpp"
00012 #include "PlayaGenericLeftPreconditioner.hpp"
00013 #include "PlayaGenericRightPreconditioner.hpp"
00014 
00015 
00016 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00017 #include "PlayaVectorImpl.hpp"
00018 #include "PlayaLinearOperatorImpl.hpp"
00019 #endif
00020 
00021 using namespace Playa;
00022 using namespace Teuchos;
00023 
00024 
00025 EpetraMatrix::EpetraMatrix(const Epetra_CrsGraph& graph,
00026   const VectorSpace<double>& domain,
00027   const VectorSpace<double>& range)
00028   : LinearOpWithSpaces<double>(domain, range),
00029     matrix_(rcp(new Epetra_CrsMatrix(Copy, graph)))
00030 {}
00031 
00032 EpetraMatrix::EpetraMatrix(const RCP<Epetra_CrsMatrix>& mat,
00033   const VectorSpace<double>& domain,
00034   const VectorSpace<double>& range)
00035   : LinearOpWithSpaces<double>(domain, range),
00036     matrix_(mat)
00037 {}
00038 
00039 
00040 void EpetraMatrix::apply(
00041   Teuchos::ETransp applyType,
00042   const Vector<double>& in,
00043   Vector<double> out) const
00044 {
00045   const Epetra_Vector& epIn = EpetraVector::getConcrete(in);
00046   Epetra_Vector& epOut = EpetraVector::getConcrete(out);
00047   using Teuchos::rcp_dynamic_cast;
00048 
00049   bool trans = applyType == TRANS;
00050   
00051   int ierr = matrix_->Multiply(trans, epIn, epOut);
00052   TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error, 
00053     "EpetraMatrix::generalApply() detected ierr="
00054     << ierr << " in matrix_->Multiply()");
00055 }
00056 
00057 
00058 
00059 void EpetraMatrix::addToRow(int globalRowIndex,
00060   int nElemsToInsert,
00061   const int* globalColumnIndices,
00062   const double* elementValues)
00063 {
00064   int ierr = crsMatrix()->SumIntoGlobalValues(globalRowIndex,
00065     nElemsToInsert,
00066     (double*) elementValues,
00067     (int*) globalColumnIndices);
00068 
00069   TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error, 
00070     "failed to add to row " << globalRowIndex
00071     << " in EpetraMatrix::addToRow() with nnz="
00072     << nElemsToInsert 
00073     << ". Error code was " << ierr);
00074 }
00075 
00076 
00077 
00078 void EpetraMatrix::addToElementBatch(int numRows, 
00079   int rowBlockSize,
00080   const int* globalRowIndices,
00081   int numColumnsPerRow,
00082   const int* globalColumnIndices,
00083   const double* values,
00084   const int* skipRow)
00085 {
00086   Epetra_CrsMatrix* crs = crsMatrix();
00087 
00088   int numRowBlocks = numRows/rowBlockSize;
00089   int row = 0;
00090 
00091   for (int rb=0; rb<numRowBlocks; rb++)
00092   {
00093     const int* cols = globalColumnIndices + rb*numColumnsPerRow;
00094     for (int r=0; r<rowBlockSize; r++, row++)
00095     {
00096       if (skipRow[row]) continue;
00097       const double* rowVals = values + row*numColumnsPerRow;
00098       int ierr=crs->SumIntoGlobalValues(globalRowIndices[row], 
00099         numColumnsPerRow,
00100         (double*) rowVals,
00101         (int*) cols);
00102       TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error, 
00103         "failed to add to row " << globalRowIndices[row]
00104         << " in EpetraMatrix::addToRow() with nnz="
00105         << numColumnsPerRow
00106         << ". Error code was " << ierr);
00107     }
00108   }
00109 }
00110 
00111 
00112 void EpetraMatrix::zero()
00113 {
00114   crsMatrix()->PutScalar(0.0);
00115 }
00116 
00117 
00118 
00119 void EpetraMatrix::getILUKPreconditioner(int fillLevels,
00120   int overlapFill,
00121   double relaxationValue,
00122   double relativeThreshold,
00123   double absoluteThreshold,
00124   LeftOrRight leftOrRight,
00125   Preconditioner<double>& rtn) const
00126 {
00127   RCP<LinearOperatorBase<double> > a = rcp(new IfpackOperator(this, 
00128       fillLevels,
00129       overlapFill,
00130       relaxationValue,
00131       relativeThreshold,
00132       absoluteThreshold));
00133 
00134   LinearOperator<double> ilu = a;
00135 
00136 
00137   if (leftOrRight == Left)
00138   {
00139     rtn = new GenericLeftPreconditioner<double>(ilu);
00140   }
00141   else
00142   {
00143     rtn = new GenericRightPreconditioner<double>(ilu);
00144   }
00145 }
00146 
00147 
00148 void EpetraMatrix::print(std::ostream& os) const 
00149 {
00150   crsMatrix()->Print(os);
00151 }
00152 
00153 
00154 string EpetraMatrix::description() const 
00155 {
00156   std::string rtn = "EpetraMatrix[nRow=" 
00157     + Teuchos::toString(crsMatrix()->NumGlobalRows())
00158     + ", nCol=" + Teuchos::toString(crsMatrix()->NumGlobalCols())
00159     + "]";
00160   return rtn;
00161 }
00162 
00163 
00164 Epetra_CrsMatrix* EpetraMatrix::crsMatrix()
00165 {
00166   return matrix_.get();
00167 }
00168 
00169 
00170 const Epetra_CrsMatrix* EpetraMatrix::crsMatrix() const 
00171 {
00172   return matrix_.get();
00173 }
00174 
00175 
00176 Epetra_CrsMatrix& EpetraMatrix::getConcrete(const LinearOperator<double>& A)
00177 {
00178   return *Teuchos::dyn_cast<EpetraMatrix>(*A.ptr()).crsMatrix();
00179 }
00180 
00181 
00182 RCP<const Epetra_CrsMatrix>
00183 EpetraMatrix::getConcretePtr(const LinearOperator<double>& A)
00184 {
00185   return Teuchos::rcp_dynamic_cast<EpetraMatrix>(A.ptr())->matrix_;
00186 }
00187 
00188 
00189 void EpetraMatrix::getRow(const int& row, 
00190   Teuchos::Array<int>& indices, 
00191   Teuchos::Array<double>& values) const
00192 {
00193   const Epetra_CrsMatrix* crs = crsMatrix();
00194 
00195   int numEntries;
00196   int* epIndices;
00197   double* epValues;
00198 
00199   int info = crs->ExtractGlobalRowView(row, numEntries, epValues, epIndices);
00200   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
00201     "call to ExtractGlobalRowView not successful");
00202 
00203   indices.resize(numEntries);
00204   values.resize(numEntries);
00205   for (int i = 0; i < numEntries; i++)
00206   {
00207     indices[i] = *epIndices;
00208     values[i] = *epValues;
00209     epIndices++;
00210     epValues++;
00211   }
00212 }
00213 
00214 
00215 

doxygen