00001
00002
00003
00004
00005 #include "PlayaEpetraMatrix.hpp"
00006 #include "PlayaEpetraVector.hpp"
00007 #include "PlayaVectorSpaceDecl.hpp"
00008 #include "PlayaVectorDecl.hpp"
00009 #include "PlayaLinearOperatorDecl.hpp"
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