3#include "Matrix_decl.hpp"
16template <
class SC,
class LO,
class GO,
class NO>
17Matrix<SC,LO,GO,NO>::Matrix():
23template <
class SC,
class LO,
class GO,
class NO>
24Matrix<SC,LO,GO,NO>::Matrix( TpetraMatrixPtr_Type& tpetraMatPtrIn ):
25matrix_( tpetraMatPtrIn )
30template <
class SC,
class LO,
class GO,
class NO>
31Matrix<SC,LO,GO,NO>::Matrix( MapConstPtr_Type map , LO numEntries):
34 matrix_ = Teuchos::RCP(
new TpetraMatrix_Type(map->getTpetraMap(), numEntries));
39template <
class SC,
class LO,
class GO,
class NO>
40Matrix<SC,LO,GO,NO>::Matrix( MatrixPtr_Type matrixIn ):
43 matrix_ =Teuchos::RCP(
new TpetraMatrix_Type( matrixIn->getMap()->getTpetraMap(), matrixIn->getGlobalMaxNumRowEntries() ));
44 if(matrixIn->isLocallyIndexed())
46 Teuchos::ArrayView<const SC> values;
47 Teuchos::ArrayView<const LO> indices;
49 MapConstPtr_Type colMap = matrixIn->getMap(
"col");
50 MapConstPtr_Type rowMap = matrixIn->getMap(
"row");
52 for (UN i=0; i<matrixIn->getNodeNumRows(); i++)
54 matrixIn->getLocalRowView( i, indices, values );
55 Teuchos::Array<GO> indicesGlobal( indices.size() );
56 for (UN j=0; j<indices.size(); j++)
58 indicesGlobal[j] = colMap->getGlobalElement( indices[j] );
61 matrix_->insertGlobalValues( rowMap->getGlobalElement( i ), indicesGlobal(), values );
65 matrix_->fillComplete( matrixIn->getMap(
"domain")->getTpetraMap(), matrixIn->getMap(
"range")->getTpetraMap() );
68template <
class SC,
class LO,
class GO,
class NO>
69Matrix<SC,LO,GO,NO>::~Matrix(){
73template <
class SC,
class LO,
class GO,
class NO>
76 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,
"");
77 matrix_->insertGlobalValues( globalRow, cols, vals );
80template <
class SC,
class LO,
class GO,
class NO>
82 return matrix_->getLocalNumRows();
86template <
class SC,
class LO,
class GO,
class NO>
89 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,
"RCP<Matrix> is null.");
90 TpetraMapConstPtr_Type tpetraMap;
91 if (!map_string.compare(
"row")) {
92 tpetraMap = matrix_->getRowMap();
94 else if (!map_string.compare(
"col")) {
95 tpetraMap = matrix_->getColMap();
97 else if (!map_string.compare(
"domain")) {
98 tpetraMap = matrix_->getDomainMap();
100 else if (!map_string.compare(
"range")) {
101 tpetraMap = matrix_->getRangeMap();
104 tpetraMap = matrix_->getMap();
106 return Teuchos::rcp(
new Map_Type(tpetraMap) );
109template <
class SC,
class LO,
class GO,
class NO>
112 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,
"RCP<Matrix> is null.");
113 TpetraMapConstPtr_Type tpetraMap;
114 if (!map_string.compare(
"row")) {
115 tpetraMap = matrix_->getRowMap();
117 else if (!map_string.compare(
"col")) {
118 tpetraMap = matrix_->getColMap();
120 else if (!map_string.compare(
"domain")) {
121 tpetraMap = matrix_->getDomainMap();
123 else if (!map_string.compare(
"range")) {
124 tpetraMap = matrix_->getRangeMap();
127 tpetraMap = matrix_->getMap();
129 return Teuchos::rcp(
new Map_Type(tpetraMap) );
133template <
class SC,
class LO,
class GO,
class NO>
136 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,
"RCP<Matrix> is null.");
138 if (!map_string.compare(
"row")) {
139 return matrix_->getRowMap();
141 else if (!map_string.compare(
"col")) {
142 return matrix_->getColMap();
144 else if (!map_string.compare(
"domain")) {
145 return matrix_->getDomainMap();
147 else if (!map_string.compare(
"range")) {
148 return matrix_->getRangeMap();
150 return matrix_->getMap();
153template <
class SC,
class LO,
class GO,
class NO>
158 Teuchos::RCP<Thyra::LinearOpBase<SC> > thyraOp = Teuchos::null;
160 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpCrsMat = matrix_;
161 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
162 Teuchos::RCP<Tpetra::RowMatrix<SC,LO,GO,NO> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<SC,LO,GO,NO> >(tpCrsMat,
true);
163 Teuchos::RCP<Tpetra::Operator <SC,LO,GO,NO> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<SC,LO,GO,NO> >(tpRowMat,
true);
165 thyraOp = Thyra::createLinearOp(tpOperator);
171template <
class SC,
class LO,
class GO,
class NO>
175 Teuchos::RCP<Thyra::LinearOpBase<SC> > thyraOp = Teuchos::null;
177 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpCrsMat = matrix_;
178 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
179 Teuchos::RCP<Tpetra::RowMatrix<SC,LO,GO,NO> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<SC,LO,GO,NO> >(tpCrsMat,
true);
180 Teuchos::RCP<Tpetra::Operator <SC,LO,GO,NO> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<SC,LO,GO,NO> >(tpRowMat,
true);
182 thyraOp = Thyra::createLinearOp(tpOperator);
184 return Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> > (thyraOp);
187template <
class SC,
class LO,
class GO,
class NO>
190 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
191 matrix_->describe(*out,verbLevel);
194template <
class SC,
class LO,
class GO,
class NO>
196 matrix_->resumeFill();
199template <
class SC,
class LO,
class GO,
class NO>
201 matrix_->fillComplete();
204template <
class SC,
class LO,
class GO,
class NO>
206 matrix_->fillComplete( domainMap->getTpetraMap(), rangeMap->getTpetraMap() );
209template <
class SC,
class LO,
class GO,
class NO>
211 return matrix_->isFillComplete();
214template <
class SC,
class LO,
class GO,
class NO>
215bool Matrix<SC,LO,GO,NO>::isLocallyIndexed(){
216 return matrix_->isLocallyIndexed();
219template <
class SC,
class LO,
class GO,
class NO>
228 Tpetra::MatrixMatrix::Multiply( *tpA->matrix_, transposeA , *tpB->matrix_, transposeB, *matrix_,
fillComplete, std::string(), Teuchos::null);
232template <
class SC,
class LO,
class GO,
class NO>
234 MatrixPtr_Type matrix(
new Matrix_Type( matrix_ ));
235 MatrixPtr_Type diagInverse(
new Matrix_Type( matrix->getMap(
"row"), 1) );
236 MapConstPtr_Type colMap = matrix->getMap(
"col");
237 MapConstPtr_Type rowMap = matrix->getMap(
"row");
240 if(diagonalType ==
"Diagonal")
242 for(
int i =0; i< rowMap->getNodeNumElements(); i ++){
243 Teuchos::ArrayView<const SC> valuesOld;
244 Teuchos::ArrayView<const LO> indices;
245 matrix->getLocalRowView(i, indices, valuesOld);
247 GO globalDof = rowMap->getGlobalElement( i );
249 Teuchos::Array<SC> values( 1, 0);
250 Teuchos::Array<GO> indicesGO( 1 , 0 );
253 for (UN j=0; j<indices.size() && !setOne; j++) {
254 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
255 values[0] = 1./valuesOld[j];
256 indicesGO[0] = colMap->getGlobalElement(indices[j]);
260 GO row = GO ( rowMap->getGlobalElement( i) );
261 diagInverse->insertGlobalValues( row, indicesGO(), values() );
264 else if(diagonalType ==
"AbsRowSum")
266 for(
int i =0; i< rowMap->getNodeNumElements(); i ++){
267 Teuchos::ArrayView<const SC> valuesOld;
268 Teuchos::ArrayView<const LO> indices;
269 matrix->getLocalRowView(i, indices, valuesOld);
271 GO globalDof = rowMap->getGlobalElement( i );
274 for (UN j=0; j<indices.size(); j++) {
275 rowSum += abs(valuesOld[j]);
278 Teuchos::Array<SC> values( 1, 0);
279 Teuchos::Array<GO> indicesGO( 1 , 0 );
282 for (UN j=0; j<indices.size() && !setOne; j++) {
283 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
284 values[0] = 1./rowSum;
285 indicesGO[0] = colMap->getGlobalElement(indices[j]);
290 GO row = GO ( rowMap->getGlobalElement( i) );
291 diagInverse->insertGlobalValues( row, indicesGO(), values() );
295 diagInverse->fillComplete();
305template <
class SC,
class LO,
class GO,
class NO>
307 TEUCHOS_TEST_FOR_EXCEPTION(matrix_->isLocallyIndexed(),std::logic_error,
"Underlying matrix is locally indexed and we can not use a global row view. Global row copy is valid here and can be implemented.");
309 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::global_inds_host_view_type Indices;
310 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type Values;
314 matrix_->getGlobalRowView(globalRow, Indices, Values);
315 indices = Teuchos::ArrayView<const GO> (Indices.data(), Indices.extent(0));
316 values = Teuchos::ArrayView<const SC> (
reinterpret_cast<const SC*
>(Values.data()), Values.extent(0));
320template <
class SC,
class LO,
class GO,
class NO>
322 TEUCHOS_TEST_FOR_EXCEPTION(matrix_->isGloballyIndexed(),std::logic_error,
"Underlying matrix is globally indexed and we can not use a local row view.");
323 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type Indices;
324 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type Values;
328 matrix_->getLocalRowView(localRow, Indices, Values);
329 indices = Teuchos::ArrayView<const LO> (Indices.data(), Indices.extent(0));
330 values = Teuchos::ArrayView<const SC> (
reinterpret_cast<const SC*
>(Values.data()), Values.extent(0));
335template <
class SC,
class LO,
class GO,
class NO>
337 matrix_->replaceGlobalValues(globalRow, indices, values);
340template <
class SC,
class LO,
class GO,
class NO>
342 matrix_->replaceLocalValues(localRow, indices, values);
345template <
class SC,
class LO,
class GO,
class NO>
350template <
class SC,
class LO,
class GO,
class NO>
353 Teuchos::ETransp mode,
357 TpetraMVPtr_Type yTpetra = Y.getTpetraMultiVectorNonConst();
359 matrix_->apply( *X.getTpetraMultiVector(), *yTpetra, mode, alpha, beta );
362template <
class SC,
class LO,
class GO,
class NO>
365 matrix_->scale( alpha );
369template <
class SC,
class LO,
class GO,
class NO>
371 TEUCHOS_TEST_FOR_EXCEPTION( matrix_.is_null(), std::runtime_error,
"Matrix in writeMM is null.");
378 TpetraMatrixPtr_Type tpetraMat = matrix_;
380 Tpetra::MatrixMarket::Writer< TpetraMatrix_Type > tpetraWriter;
382 tpetraWriter.writeSparseFile(fileName, tpetraMat,
"matrix",
"");
385template <
class SC,
class LO,
class GO,
class NO>
388 if (B->isFillComplete())
390 TEUCHOS_TEST_FOR_EXCEPTION( B->isLocallyIndexed(), std::runtime_error,
"Matrix in is locally index but Trilinos Tpetra cannot add to a matrix at this stage.");
395 Tpetra::MatrixMatrix::Add(*matrix_,
false, alpha, *B->matrix_, beta);
400template <
class SC,
class LO,
class GO,
class NO>
403 MapConstPtr_Type map = this->
getMap(
"row");
404 MapConstPtr_Type mapCol = this->
getMap(
"col");
406 int numMV = mapCol->getMaxAllGlobalIndex() + 1;
408 mv = Teuchos::rcp(
new MultiVector_Type ( map, numMV ) );
409 Teuchos::ArrayView< const LO > indices;
410 Teuchos::ArrayView< const SC > values;
411 Teuchos::ArrayRCP< SC > valuesMV;
415 for (
int j=0; j<indices.size(); j++) {
416 globalCol = mapCol->getGlobalElement( indices[j] );
417 valuesMV = mv->getDataNonConst( globalCol );
418 valuesMV[ i ] = values[ j ];
422template <
class SC,
class LO,
class GO,
class NO>
424 return matrix_->getGlobalMaxNumRowEntries();
427template <
class SC,
class LO,
class GO,
class NO>
428void Matrix<SC,LO,GO,NO>::insertLocalValues(LO localRow,
const Teuchos::ArrayView< const LO > &cols,
const Teuchos::ArrayView< const SC > &vals){
430 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,
"");
431 matrix_->insertLocalValues( localRow, cols, vals );
435template <
class SC,
class LO,
class GO,
class NO>
440 if ( importer_.is_null() || !reuseImport) {
442 importer_ = Teuchos::RCP(
new Tpetra::Import<LO, GO, NO>( mvIn->getMapTpetra(), this->getMapTpetra() ));
443 else if(type==
"Reverse")
444 importer_ = Teuchos::RCP(
new Tpetra::Import<LO, GO, NO>( this->
getMapTpetra(), mvIn->getMapTpetra() ));
446 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown type for import. Choose Forward or Reverse");
449 TEUCHOS_TEST_FOR_EXCEPTION( !importer_->getSourceMap()->isSameAs( *mvIn->getMap()->getTpetraMap() ), std::logic_error,
"Source maps of Importer and Matrix are not the same.");
450 TEUCHOS_TEST_FOR_EXCEPTION( !importer_->getTargetMap()->isSameAs( *this->getMap()->getTpetraMap() ), std::logic_error,
"Target maps of Importer and Matrix are not the same.");
454 if (type==
"Forward") {
455 if ( !combineMode.compare(
"Insert") )
456 matrix_->doImport ( *mvIn->getTpetraMatrix(), *importer_, Tpetra::INSERT);
457 else if ( !combineMode.compare(
"Add") )
458 matrix_->doImport ( *mvIn->getTpetraMatrix(), *importer_, Tpetra::ADD);
460 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown combine mode.");
462 else if(type==
"Reverse"){
463 if ( !combineMode.compare(
"Insert") )
464 matrix_->doExport ( *mvIn->getTpetraMatrix(), *importer_, Tpetra::INSERT);
465 else if ( !combineMode.compare(
"Add") )
466 matrix_->doExport ( *mvIn->getTpetraMatrix(), *importer_, Tpetra::ADD);
468 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown combine mode.");
471 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown type for import. Choose Forward or Reverse");
474template <
class SC,
class LO,
class GO,
class NO>
479 if ( exporter_.is_null() || !reuseExport) {
481 exporter_ = Teuchos::RCP(
new Tpetra::Export<LO, GO, NO>( mvIn->getMapTpetra(), this->getMapTpetra() ));
482 else if(type==
"Reverse")
483 exporter_ = Teuchos::RCP(
new Tpetra::Export<LO, GO, NO>( this->
getMapTpetra(), mvIn->getMapTpetra() ));
485 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown type for import. Choose Forward or Reverse");
488 TEUCHOS_TEST_FOR_EXCEPTION( !exporter_->getSourceMap()->isSameAs( *this->getMap()->getTpetraMap() ), std::logic_error,
"Source maps of Exporter and Multivector are not the same.");
489 TEUCHOS_TEST_FOR_EXCEPTION( !exporter_->getTargetMap()->isSameAs( *mvIn->getMap()->getTpetraMap() ), std::logic_error,
"Target maps of Exporter and Multivector are not the same.");
493 if (type==
"Forward") {
494 if ( !combineMode.compare(
"Insert") )
495 matrix_->doExport ( *mvIn->getTpetraMatrix(), *exporter_, Tpetra::INSERT);
496 else if ( !combineMode.compare(
"Add") )
497 matrix_->doExport ( *mvIn->getTpetraMatrix(), *exporter_, Tpetra::ADD);
499 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown combine mode.");
501 else if(type==
"Reverse"){
502 if ( !combineMode.compare(
"Insert") )
503 matrix_->doImport ( *mvIn->getTpetraMatrix(), *exporter_, Tpetra::INSERT);
504 else if ( !combineMode.compare(
"Add") )
505 matrix_->doImport ( *mvIn->getTpetraMatrix(), *exporter_, Tpetra::ADD);
507 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown combine mode.");
510 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Unknown type for import. Choose Forward or Reverse");
void getLocalRowView(LO localRow, Teuchos::ArrayView< const LO > &indices, Teuchos::ArrayView< const SC > &values) const
Extracting single rows of Matrix with local row ID. Indices returns local indices of entries stored i...
Definition Matrix_def.hpp:321
MatrixPtr_Type buildDiagonalInverse(std::string diagonalType)
Definition Matrix_def.hpp:233
Teuchos::RCP< const Thyra::LinearOpBase< SC > > getThyraLinOp() const
i.e. for NOX
Definition Matrix_def.hpp:154
LO getNodeNumRows() const
Returns the local number of rows.
Definition Matrix_def.hpp:81
void print(Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_EXTREME)
printing matrix
Definition Matrix_def.hpp:188
void replaceLocalValues(LO localRow, const Teuchos::ArrayView< const LO > &indices, const Teuchos::ArrayView< const SC > &values)
Replacing single rows of Matrix with local row ID. Indices returns local indices of entries stored in...
Definition Matrix_def.hpp:341
void replaceGlobalValues(GO globalRow, const Teuchos::ArrayView< const GO > &indices, const Teuchos::ArrayView< const SC > &values)
Replacing single rows of Matrix with global row ID. Indices returns global indices of entries stored ...
Definition Matrix_def.hpp:336
void exportFromVector(MatrixPtr_Type mvIn, bool reuseExport=false, std::string combineMode="Insert", std::string type="Forward")
Definition Matrix_def.hpp:475
void getGlobalRowView(GO globalRow, Teuchos::ArrayView< const GO > &indices, Teuchos::ArrayView< const SC > &values) const
Extracting single rows of Matrix with global row ID. Indices returns global indices of entries stored...
Definition Matrix_def.hpp:306
void toMV(MultiVectorPtr_Type &mv)
Turning Matrix into MultiVector Format.
Definition Matrix_def.hpp:401
void scale(const SC &alpha)
Scaling this with constant alpha.
Definition Matrix_def.hpp:363
TpetraMapConstPtr_Type getMapTpetra(std::string map_string="")
Return map in Xpetra Format of type " ".
Definition Matrix_def.hpp:134
MapConstPtr_Type getMap(std::string map_string="")
Returns map of type " ". i.e. row or column map.
Definition Matrix_def.hpp:87
Teuchos::RCP< Thyra::LinearOpBase< SC > > getThyraLinOpNonConst()
i.e. for NOX
Definition Matrix_def.hpp:172
void insertGlobalValues(GO globalRow, const Teuchos::ArrayView< const GO > &cols, const Teuchos::ArrayView< const SC > &vals)
Intertion of values in global row 'globalRow'. Matrix is distributed nodewise. If values are added to...
Definition Matrix_def.hpp:74
bool isFillComplete()
Check if matrix is already filled complete.
Definition Matrix_def.hpp:210
void writeMM(std::string fileName="matrix.mm") const
Writing Matrix in file.
Definition Matrix_def.hpp:370
void importFromVector(MatrixPtr_Type mvIn, bool reuseImport=false, std::string combineMode="Insert", std::string type="Forward")
Definition Matrix_def.hpp:436
void resumeFill()
Resuming filling process. But only limited options i.e. scaling remain.
Definition Matrix_def.hpp:195
void addMatrix(SC alpha, const MatrixPtr_Type &B, SC beta)
B = alpha*this + beta*B.
Definition Matrix_def.hpp:386
void apply(const MultiVector_Type &X, MultiVector_Type &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Matrix Vector Operation. Applying MultiVector X to this. Y = alpha * (this)^mode * X + beta * Y....
Definition Matrix_def.hpp:351
void fillComplete()
after inserting global values into matrix. After this step the column map is fixed....
Definition Matrix_def.hpp:200
void Multiply(const MatrixPtr_Type &tpA, bool transposeA, const MatrixPtr_Type &tpB, bool transposeB, bool fillComplete=true)
Definition Matrix_def.hpp:220
LO getGlobalMaxNumRowEntries() const
Maximum number of entries in any row of the matrix, over all processes.
Definition Matrix_def.hpp:423
TpetraMatrixConstPtr_Type getTpetraMatrix() const
Return matrix in Xpetra Format of type " ".
Definition Matrix_def.hpp:346
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:33