Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
Matrix_def.hpp
1#ifndef MATRIX_DEF_hpp
2#define MATRIX_DEF_hpp
3#include "Matrix_decl.hpp"
4
13
14namespace FEDD {
15//using namespace Teuchos;
16template <class SC, class LO, class GO, class NO>
17Matrix<SC,LO,GO,NO>::Matrix():
18 matrix_()
19{
20
21}
22
23template <class SC, class LO, class GO, class NO>
24Matrix<SC,LO,GO,NO>::Matrix( TpetraMatrixPtr_Type& tpetraMatPtrIn ):
25matrix_( tpetraMatPtrIn )
26{
27
28}
29
30template <class SC, class LO, class GO, class NO>
31Matrix<SC,LO,GO,NO>::Matrix( MapConstPtr_Type map , LO numEntries):
32matrix_( )
33{
34 matrix_ = Teuchos::RCP( new TpetraMatrix_Type(map->getTpetraMap(), numEntries)); //, plist) );
35
36}
37
38
39template <class SC, class LO, class GO, class NO>
40Matrix<SC,LO,GO,NO>::Matrix( MatrixPtr_Type matrixIn ):
41matrix_( )
42{
43 matrix_ =Teuchos::RCP( new TpetraMatrix_Type( matrixIn->getMap()->getTpetraMap(), matrixIn->getGlobalMaxNumRowEntries() ));
44 if(matrixIn->isLocallyIndexed())
45 {
46 Teuchos::ArrayView<const SC> values;
47 Teuchos::ArrayView<const LO> indices;
48
49 MapConstPtr_Type colMap = matrixIn->getMap("col");
50 MapConstPtr_Type rowMap = matrixIn->getMap("row");
51
52 for (UN i=0; i<matrixIn->getNodeNumRows(); i++)
53 {
54 matrixIn->getLocalRowView( i, indices, values );
55 Teuchos::Array<GO> indicesGlobal( indices.size() );
56 for (UN j=0; j<indices.size(); j++)
57 {
58 indicesGlobal[j] = colMap->getGlobalElement( indices[j] );
59
60 }
61 matrix_->insertGlobalValues( rowMap->getGlobalElement( i ), indicesGlobal(), values );
62 }
63 }
64
65 matrix_->fillComplete( matrixIn->getMap("domain")->getTpetraMap(), matrixIn->getMap("range")->getTpetraMap() );
66}
67
68template <class SC, class LO, class GO, class NO>
69Matrix<SC,LO,GO,NO>::~Matrix(){
70
71}
72
73template <class SC, class LO, class GO, class NO>
74void Matrix<SC,LO,GO,NO>::insertGlobalValues(GO globalRow, const Teuchos::ArrayView< const GO > &cols, const Teuchos::ArrayView< const SC > &vals){
75
76 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,"");
77 matrix_->insertGlobalValues( globalRow, cols, vals );
78}
79
80template <class SC, class LO, class GO, class NO>
82 return matrix_->getLocalNumRows();
83}
84
85
86template <class SC, class LO, class GO, class NO>
87typename Matrix<SC,LO,GO,NO>::MapConstPtr_Type Matrix<SC,LO,GO,NO>::getMap(std::string map_string){
88
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();
93 }
94 else if (!map_string.compare("col")) {
95 tpetraMap = matrix_->getColMap();
96 }
97 else if (!map_string.compare("domain")) {
98 tpetraMap = matrix_->getDomainMap();
99 }
100 else if (!map_string.compare("range")) {
101 tpetraMap = matrix_->getRangeMap();
102 }
103 else
104 tpetraMap = matrix_->getMap();
105
106 return Teuchos::rcp( new Map_Type(tpetraMap) );
108
109template <class SC, class LO, class GO, class NO>
110typename Matrix<SC,LO,GO,NO>::MapConstPtr_Type Matrix<SC,LO,GO,NO>::getMap(std::string map_string) const{
111
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();
116 }
117 else if (!map_string.compare("col")) {
118 tpetraMap = matrix_->getColMap();
119 }
120 else if (!map_string.compare("domain")) {
121 tpetraMap = matrix_->getDomainMap();
123 else if (!map_string.compare("range")) {
124 tpetraMap = matrix_->getRangeMap();
125 }
126 else
127 tpetraMap = matrix_->getMap();
128
129 return Teuchos::rcp( new Map_Type(tpetraMap) );
130}
131
133template <class SC, class LO, class GO, class NO>
134typename Matrix<SC,LO,GO,NO>::TpetraMapConstPtr_Type Matrix<SC,LO,GO,NO>::getMapTpetra(std::string map_string){
135
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();
140 }
141 else if (!map_string.compare("col")) {
142 return matrix_->getColMap();
143 }
144 else if (!map_string.compare("domain")) {
145 return matrix_->getDomainMap();
146 }
147 else if (!map_string.compare("range")) {
148 return matrix_->getRangeMap();
149 }
150 return matrix_->getMap();
151}
153template <class SC, class LO, class GO, class NO>
154Teuchos::RCP<const Thyra::LinearOpBase<SC> > Matrix<SC,LO,GO,NO>::getThyraLinOp() const{
155 //Xpetra::CrsMatrixWrap<SC,LO,GO,NO>& crsWrapMatrix = dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO>&>(*matrix_);
156 //return Xpetra::ThyraUtils<SC,LO,GO,NO>::toThyra(crsWrapMatrix.getCrsMatrix());
157
158 Teuchos::RCP<Thyra::LinearOpBase<SC> > thyraOp = Teuchos::null;
159
160 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpCrsMat = matrix_; // = xTpCrsMat->getTpetra_CrsMatrixNonConst();
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);
164
165 thyraOp = Thyra::createLinearOp(tpOperator);
166
167 return thyraOp;
168
169}
171template <class SC, class LO, class GO, class NO>
172Teuchos::RCP<Thyra::LinearOpBase<SC> > Matrix<SC,LO,GO,NO>::getThyraLinOpNonConst() {
173 //Xpetra::CrsMatrixWrap<SC,LO,GO,NO>& crsWrapMatrix = dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO>&>(*matrix_);
174 //return Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> > (Xpetra::ThyraUtils<SC,LO,GO,NO>::toThyra(crsWrapMatrix.getCrsMatrix()) );
175 Teuchos::RCP<Thyra::LinearOpBase<SC> > thyraOp = Teuchos::null;
176
177 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpCrsMat = matrix_; // = xTpCrsMat->getTpetra_CrsMatrixNonConst();
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);
181
182 thyraOp = Thyra::createLinearOp(tpOperator);
183
184 return Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> > (thyraOp); // is this now const or not??
186
187template <class SC, class LO, class GO, class NO>
188void Matrix<SC,LO,GO,NO>::print(Teuchos::EVerbosityLevel verbLevel){
189
190 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
191 matrix_->describe(*out,verbLevel);
192}
193
194template <class SC, class LO, class GO, class NO>
196 matrix_->resumeFill();
197}
198
199template <class SC, class LO, class GO, class NO>
201 matrix_->fillComplete();
202}
203
204template <class SC, class LO, class GO, class NO>
205void Matrix<SC,LO,GO,NO>::fillComplete(MapConstPtr_Type domainMap, MapConstPtr_Type rangeMap){
206 matrix_->fillComplete( domainMap->getTpetraMap(), rangeMap->getTpetraMap() );
207}
208
209template <class SC, class LO, class GO, class NO>
211 return matrix_->isFillComplete();
212}
213
214template <class SC, class LO, class GO, class NO>
215bool Matrix<SC,LO,GO,NO>::isLocallyIndexed(){
216 return matrix_->isLocallyIndexed();
217}
218
219template <class SC, class LO, class GO, class NO>
220void Matrix<SC,LO,GO,NO>::Multiply( const MatrixPtr_Type &tpA, bool transposeA , const MatrixPtr_Type &tpB, bool transposeB,bool fillComplete){
221
222 // Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = A;
223
224 // Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = B;
225
226 // Tpetra::CrsMatrix<SC,LO,GO,NO>& C;
227
228 Tpetra::MatrixMatrix::Multiply( *tpA->matrix_, transposeA , *tpB->matrix_, transposeB, *matrix_, fillComplete, std::string(), Teuchos::null);
230
231
232template <class SC,class LO,class GO,class NO>
233typename Matrix<SC,LO,GO,NO>::MatrixPtr_Type Matrix<SC,LO,GO,NO>::buildDiagonalInverse( std::string diagonalType){
234 MatrixPtr_Type matrix(new Matrix_Type( matrix_ )); // Diagonal matrix
235 MatrixPtr_Type diagInverse(new Matrix_Type( matrix->getMap("row"), 1) ); // Diagonal matrix
236 MapConstPtr_Type colMap = matrix->getMap("col");
237 MapConstPtr_Type rowMap = matrix->getMap("row");
238
239
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 );
248
249 Teuchos::Array<SC> values( 1, 0);
250 Teuchos::Array<GO> indicesGO( 1 , 0 );
251 bool setOne = false;
252
253 for (UN j=0; j<indices.size() && !setOne; j++) {
254 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
255 values[0] = 1./valuesOld[j]; // Diagonal Value
256 indicesGO[0] = colMap->getGlobalElement(indices[j]);
257 setOne=true;
258 }
259 }
260 GO row = GO ( rowMap->getGlobalElement( i) );
261 diagInverse->insertGlobalValues( row, indicesGO(), values() );
262 }
263 }
264 else if(diagonalType == "AbsRowSum")
265 {
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);
270
271 GO globalDof = rowMap->getGlobalElement( i );
272
273 double rowSum = 0.;
274 for (UN j=0; j<indices.size(); j++) {
275 rowSum += abs(valuesOld[j]);
276 }
277
278 Teuchos::Array<SC> values( 1, 0);
279 Teuchos::Array<GO> indicesGO( 1 , 0 );
280 bool setOne = false;
281
282 for (UN j=0; j<indices.size() && !setOne; j++) {
283 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
284 values[0] = 1./rowSum; // Diagonal Value
285 indicesGO[0] = colMap->getGlobalElement(indices[j]);
286 setOne = true;
287
288 }
289 }
290 GO row = GO ( rowMap->getGlobalElement( i) );
291 diagInverse->insertGlobalValues( row, indicesGO(), values() );
292 }
293 }
294
295 diagInverse->fillComplete();
296 // diagInverse->writeMM("Mu_Inverse_LSC");
297 return diagInverse;
298}
299
300//typename Matrix<SC,LO,GO,NO>::ThyraLinOpPtr_Type Matrix<SC,LO,GO,NO>::getThyraLinOp(){
301//
302// return ;
303//}
304
305template <class SC, class LO, class GO, class NO>
306void Matrix<SC,LO,GO,NO>::getGlobalRowView(GO globalRow, Teuchos::ArrayView< const GO > &indices, Teuchos::ArrayView< const SC > &values) const{
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.");
308
309 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::global_inds_host_view_type Indices; //ArrayView< const LO > indices
310 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type Values;
311 //typename Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::local_inds_host_view_type indices;
312 //typename Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::values_host_view_type values;
313
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));
317
318}
319
320template <class SC, class LO, class GO, class NO>
321void Matrix<SC,LO,GO,NO>::getLocalRowView(LO localRow, Teuchos::ArrayView< const LO > &indices, Teuchos::ArrayView< const SC > &values) const{
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; //ArrayView< const LO > indices
324 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type Values;
325 //typename Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::local_inds_host_view_type indices;
326 //typename Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::values_host_view_type values;
327
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));
331
332}
333
334
335template <class SC, class LO, class GO, class NO>
336void Matrix<SC,LO,GO,NO>::replaceGlobalValues(GO globalRow, const Teuchos::ArrayView< const GO > &indices, const Teuchos::ArrayView< const SC > &values){
337 matrix_->replaceGlobalValues(globalRow, indices, values);
338}
339
340template <class SC, class LO, class GO, class NO>
341void Matrix<SC,LO,GO,NO>::replaceLocalValues(LO localRow, const Teuchos::ArrayView< const LO > &indices, const Teuchos::ArrayView< const SC > &values){
342 matrix_->replaceLocalValues(localRow, indices, values);
343}
344
345template <class SC, class LO, class GO, class NO>
346typename Matrix<SC,LO,GO,NO>::TpetraMatrixConstPtr_Type Matrix<SC,LO,GO,NO>::getTpetraMatrix() const{
347 return matrix_;
348}
349
350template <class SC, class LO, class GO, class NO>
351void Matrix<SC,LO,GO,NO>::apply(const MultiVector_Type& X,
352 MultiVector_Type& Y,
353 Teuchos::ETransp mode,
354 SC alpha,
355 SC beta) const{
356
357 TpetraMVPtr_Type yTpetra = Y.getTpetraMultiVectorNonConst();
358
359 matrix_->apply( *X.getTpetraMultiVector(), *yTpetra, mode, alpha, beta );
360}
361
362template <class SC, class LO, class GO, class NO>
363void Matrix<SC,LO,GO,NO>::scale(const SC& alpha) {
364
365 matrix_->scale( alpha );
366
367}
368
369template <class SC, class LO, class GO, class NO>
370void Matrix<SC,LO,GO,NO>::writeMM(std::string fileName) const{
371 TEUCHOS_TEST_FOR_EXCEPTION( matrix_.is_null(), std::runtime_error,"Matrix in writeMM is null.");
372 //typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMatrix;
373 //typedef Teuchos::RCP<TpetraCrsMatrix> TpetraCrsMatrixPtr;
374
375 //Xpetra::CrsMatrixWrap<SC,LO,GO,NO>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO>&>(*matrix_);
376 //Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>& xTpetraMat = dynamic_cast<Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>&>(*crsOp.getCrsMatrix());
377
378 TpetraMatrixPtr_Type tpetraMat = matrix_;
379
380 Tpetra::MatrixMarket::Writer< TpetraMatrix_Type > tpetraWriter;
381
382 tpetraWriter.writeSparseFile(fileName, tpetraMat, "matrix", "");
383}
384
385template <class SC, class LO, class GO, class NO>
386void Matrix<SC,LO,GO,NO>::addMatrix(SC alpha, const MatrixPtr_Type &B, SC beta){
387 //B = alpha*A + beta*B.
388 if (B->isFillComplete())
389 B->resumeFill();
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.");
391
392 //const Tpetra::CrsMatrix<SC,LO,GO,NO>& tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
393 //Tpetra::CrsMatrix<SC,LO,GO,NO>& tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(B);
394
395 Tpetra::MatrixMatrix::Add(*matrix_, false, alpha, *B->matrix_, beta);
396
397 //Xpetra::MatrixMatrix<SC,LO,GO,NO>::TwoMatrixAdd( *matrix_, false, alpha, *B->matrix_, beta );
398}
399
400template <class SC, class LO, class GO, class NO>
401void Matrix<SC,LO,GO,NO>::toMV( MultiVectorPtr_Type& mv ){
402
403 MapConstPtr_Type map = this->getMap("row");
404 MapConstPtr_Type mapCol = this->getMap("col");
405 // number of matrix columns will be the number of multivectors
406 int numMV = mapCol->getMaxAllGlobalIndex() + 1;
407
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;
412 GO globalCol = -1;
413 for (int i=0; i<this->getNodeNumRows(); i++) {
414 this->getLocalRowView( i, indices, values );
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 ];
419 }
420 }
421}
422template <class SC, class LO, class GO, class NO>
424 return matrix_->getGlobalMaxNumRowEntries();
425}
426
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){
429
430 TEUCHOS_TEST_FOR_EXCEPTION(matrix_.is_null(),std::runtime_error,"");
431 matrix_->insertLocalValues( localRow, cols, vals );
432}
433
434
435template <class SC, class LO, class GO, class NO>
436void Matrix<SC,LO,GO,NO>::importFromVector( MatrixPtr_Type mvIn, bool reuseImport, std::string combineMode, std::string type) {
437
438 //TEUCHOS_TEST_FOR_EXCEPTION( getNodeNumCols() != mvIn->getNodeNumCols(), std::logic_error,"MultiVectors for fillFromVector have different number of vectors.");
439
440 if ( importer_.is_null() || !reuseImport) {
441 if (type=="Forward")
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() ));
445 else
446 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown type for import. Choose Forward or Reverse");
447 }
448 else{
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.");
451 }
452
453
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);
459 else
460 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown combine mode.");
461 }
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);
467 else
468 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown combine mode.");
469 }
470 else
471 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown type for import. Choose Forward or Reverse");
472}
473
474template <class SC, class LO, class GO, class NO>
475void Matrix<SC,LO,GO,NO>::exportFromVector( MatrixPtr_Type mvIn, bool reuseExport, std::string combineMode, std::string type) {
476
477 //TEUCHOS_TEST_FOR_EXCEPTION( getNodeNumCols() != mvIn->getNodeNumCols(), std::logic_error,"MultiVectors for fillFromVector have different number of vectors.");
478
479 if ( exporter_.is_null() || !reuseExport) {
480 if (type=="Forward")
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() ));
484 else
485 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown type for import. Choose Forward or Reverse");
486 }
487 else{
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.");
490 }
491
492
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);
498 else
499 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown combine mode.");
500 }
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);
506 else
507 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown combine mode.");
508 }
509 else
510 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Unknown type for import. Choose Forward or Reverse");
511}
512
513}
514#endif
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