47 #include "Teko_Config.h" 51 #include "Thyra_MultiVectorStdOps.hpp" 52 #include "Thyra_ZeroLinearOpBase.hpp" 53 #include "Thyra_DefaultDiagonalLinearOp.hpp" 54 #include "Thyra_DefaultAddedLinearOp.hpp" 55 #include "Thyra_EpetraExtDiagScaledMatProdTransformer.hpp" 56 #include "Thyra_EpetraExtDiagScalingTransformer.hpp" 57 #include "Thyra_EpetraExtAddTransformer.hpp" 58 #include "Thyra_DefaultScaledAdjointLinearOp.hpp" 59 #include "Thyra_DefaultMultipliedLinearOp.hpp" 60 #include "Thyra_DefaultZeroLinearOp.hpp" 61 #include "Thyra_DefaultProductMultiVector.hpp" 62 #include "Thyra_DefaultProductVectorSpace.hpp" 63 #include "Thyra_MultiVectorStdOps.hpp" 64 #include "Thyra_VectorStdOps.hpp" 65 #include "Thyra_SpmdVectorBase.hpp" 66 #include "Thyra_get_Epetra_Operator.hpp" 67 #include "Thyra_EpetraThyraWrappers.hpp" 68 #include "Thyra_EpetraOperatorWrapper.hpp" 69 #include "Thyra_EpetraLinearOp.hpp" 72 #include "Teuchos_Array.hpp" 75 #include "Epetra_Operator.h" 76 #include "Epetra_CrsGraph.h" 77 #include "Epetra_CrsMatrix.h" 78 #include "Epetra_Vector.h" 79 #include "Epetra_Map.h" 81 #include "EpetraExt_Transpose_RowMatrix.h" 82 #include "EpetraExt_MatrixMatrix.h" 85 #include "AnasaziBasicEigenproblem.hpp" 86 #include "AnasaziThyraAdapter.hpp" 87 #include "AnasaziBlockKrylovSchurSolMgr.hpp" 88 #include "AnasaziBlockKrylovSchur.hpp" 89 #include "AnasaziStatusTestMaxIters.hpp" 92 #ifdef Teko_ENABLE_Isorropia 93 #include "Isorropia_EpetraProber.hpp" 97 #include "Teko_EpetraHelpers.hpp" 98 #include "Teko_EpetraOperatorWrapper.hpp" 99 #include "Teko_TpetraHelpers.hpp" 100 #include "Teko_TpetraOperatorWrapper.hpp" 103 #include "Thyra_TpetraLinearOp.hpp" 104 #include "Tpetra_CrsMatrix.hpp" 105 #include "Tpetra_Vector.hpp" 106 #include "Thyra_TpetraThyraWrappers.hpp" 107 #include "TpetraExt_MatrixMatrix.hpp" 108 #include "Tpetra_RowMatrixTransposer.hpp" 115 using Teuchos::rcp_dynamic_cast;
117 #ifdef Teko_ENABLE_Isorropia 118 using Isorropia::Epetra::Prober;
121 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
123 Teuchos::RCP<Teuchos::FancyOStream> os =
124 Teuchos::VerboseObjectBase::getDefaultOStream();
132 inline double dist(
int dim,
double * coords,
int row,
int col)
135 for(
int i=0;i<dim;i++)
136 value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
139 return std::sqrt(value);
143 inline double dist(
double * x,
double * y,
double * z,
int stride,
int row,
int col)
146 if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
147 if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
148 if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
151 return std::sqrt(value);
172 RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double * coords,
const Epetra_CrsMatrix & stencil)
175 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
176 stencil.MaxNumEntries()+1),
true);
179 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
180 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
183 for(
int j=0;j<gl->NumMyRows();j++) {
184 int row = gl->GRID(j);
185 double diagValue = 0.0;
190 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
193 for(
int i=0;i<rowSz;i++) {
197 double value = rowData[i];
198 if(value==0)
continue;
202 double d = dist(dim,coords,row,col);
204 diagValue += rowData[i];
212 rowData[rowSz] = -diagValue;
217 rowData[diagInd] = -diagValue;
218 rowInd[diagInd] = row;
222 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
230 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(
int dim,ST * coords,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
233 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
234 stencil.getGlobalMaxNumRowEntries()+1));
237 std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
238 std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
241 for(LO j=0;j< (LO) gl->getNodeNumRows();j++) {
242 GO row = gl->getRowMap()->getGlobalElement(j);
248 stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
251 for(
size_t i=0;i<rowSz;i++) {
255 ST value = rowData[i];
256 if(value==0)
continue;
260 ST d = dist(dim,coords,row,col);
262 diagValue += rowData[i];
270 rowData[rowSz] = -diagValue;
275 rowData[diagInd] = -diagValue;
276 rowInd[diagInd] = row;
280 gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
310 RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double * x,
double * y,
double * z,
int stride,
const Epetra_CrsMatrix & stencil)
313 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
314 stencil.MaxNumEntries()+1),
true);
317 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
318 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
321 for(
int j=0;j<gl->NumMyRows();j++) {
322 int row = gl->GRID(j);
323 double diagValue = 0.0;
328 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
331 for(
int i=0;i<rowSz;i++) {
335 double value = rowData[i];
336 if(value==0)
continue;
340 double d = dist(x,y,z,stride,row,col);
342 diagValue += rowData[i];
350 rowData[rowSz] = -diagValue;
355 rowData[diagInd] = -diagValue;
356 rowInd[diagInd] = row;
360 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
368 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
371 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
372 stencil.getGlobalMaxNumRowEntries()+1),
true);
375 std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
376 std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
379 for(LO j=0;j<(LO) gl->getNodeNumRows();j++) {
380 GO row = gl->getRowMap()->getGlobalElement(j);
386 stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
389 for(
size_t i=0;i<rowSz;i++) {
393 ST value = rowData[i];
394 if(value==0)
continue;
398 ST d = dist(x,y,z,stride,row,col);
400 diagValue += rowData[i];
408 rowData[rowSz] = -diagValue;
413 rowData[diagInd] = -diagValue;
414 rowInd[diagInd] = row;
418 gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
441 void applyOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha,
double beta)
443 Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
461 void applyTransposeOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha,
double beta)
463 Thyra::apply(*A,Thyra::TRANS,*x,y.ptr(),alpha,beta);
468 void update(
double alpha,
const MultiVector & x,
double beta,MultiVector & y)
470 Teuchos::Array<double>
scale;
471 Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
474 scale.push_back(alpha);
475 vec.push_back(x.ptr());
478 Thyra::linear_combination<double>(
scale,vec,beta,y.ptr());
482 BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
488 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
489 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
495 upper->beginBlockFill(rows,rows);
497 for(
int i=0;i<rows;i++) {
501 RCP<const Thyra::LinearOpBase<double> > zed
502 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
503 upper->setBlock(i,i,zed);
505 for(
int j=i+1;j<rows;j++) {
507 LinearOp uij = blo->getBlock(i,j);
510 if(uij!=Teuchos::null)
511 upper->setBlock(i,j,uij);
515 upper->endBlockFill();
521 BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
527 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
528 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
534 lower->beginBlockFill(rows,rows);
536 for(
int i=0;i<rows;i++) {
540 RCP<const Thyra::LinearOpBase<double> > zed
541 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
542 lower->setBlock(i,i,zed);
544 for(
int j=0;j<i;j++) {
546 LinearOp uij = blo->getBlock(i,j);
549 if(uij!=Teuchos::null)
550 lower->setBlock(i,j,uij);
554 lower->endBlockFill();
578 BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp & blo)
584 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
585 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
591 zeroOp->beginBlockFill(rows,rows);
593 for(
int i=0;i<rows;i++) {
597 RCP<const Thyra::LinearOpBase<double> > zed
598 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
599 zeroOp->setBlock(i,i,zed);
606 bool isZeroOp(
const LinearOp op)
609 if(op==Teuchos::null)
return true;
612 const LinearOp test = rcp_dynamic_cast<
const Thyra::ZeroLinearOpBase<double> >(op);
615 return test!=Teuchos::null;
626 ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp & op)
628 bool isTpetra =
false;
629 RCP<const Epetra_CrsMatrix> eCrsOp;
630 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
634 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
635 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
638 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
640 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
642 else if (!tOp.is_null()){
643 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
647 throw std::logic_error(
"Neither Epetra nor Tpetra");
649 catch (std::exception & e) {
650 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
652 *out <<
"Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
653 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
655 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
657 *out <<
"*** THROWN EXCEPTION ***\n";
658 *out << e.what() << std::endl;
659 *out <<
"************************\n";
666 const RCP<Epetra_Vector> ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
667 Epetra_Vector & diag = *ptrDiag;
671 for(
int i=0;i<eCrsOp->NumMyRows();i++) {
674 eCrsOp->ExtractMyRowView(i,numEntries,values);
677 for(
int j=0;j<numEntries;j++)
678 diag[i] += std::abs(values[j]);
682 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),
"absRowSum( " + op->getObjectLabel() +
" )");
686 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
687 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
691 for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
692 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
693 std::vector<LO> indices(numEntries);
694 std::vector<ST> values(numEntries);
695 Teuchos::ArrayView<const LO> indices_av(indices);
696 Teuchos::ArrayView<const ST> values_av(values);
697 tCrsOp->getLocalRowView(i,indices_av,values_av);
700 for(LO j=0;j<numEntries;j++)
701 diag.sumIntoLocalValue(i,std::abs(values_av[j]));
705 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),
"absRowSum( " + op->getObjectLabel() +
" ))");
719 ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp & op)
721 if(Teko::TpetraHelpers::isTpetraLinearOp(op)) {
724 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
727 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
728 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
732 for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
733 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
734 std::vector<LO> indices(numEntries);
735 std::vector<ST> values(numEntries);
736 Teuchos::ArrayView<const LO> indices_av(indices);
737 Teuchos::ArrayView<const ST> values_av(values);
738 tCrsOp->getLocalRowView(i,indices_av,values_av);
741 for(LO j=0;j<numEntries;j++)
742 diag.sumIntoLocalValue(i,std::abs(values_av[j]));
745 diag.reciprocal(diag);
748 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),
"absRowSum( " + op->getObjectLabel() +
" ))");
752 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op,
true);
753 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
756 const RCP<Epetra_Vector> ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
757 Epetra_Vector & diag = *ptrDiag;
761 for(
int i=0;i<eCrsOp->NumMyRows();i++) {
764 eCrsOp->ExtractMyRowView(i,numEntries,values);
767 for(
int j=0;j<numEntries;j++)
768 diag[i] += std::abs(values[j]);
770 diag.Reciprocal(diag);
773 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),
"absRowSum( " + op->getObjectLabel() +
" )");
785 ModifiableLinearOp getLumpedMatrix(
const LinearOp & op)
787 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
788 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
791 Thyra::assign(ones.ptr(),1.0);
795 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
797 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
808 ModifiableLinearOp getInvLumpedMatrix(
const LinearOp & op)
810 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
811 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
814 Thyra::assign(ones.ptr(),1.0);
817 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
818 Thyra::reciprocal(*diag,diag.ptr());
820 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
834 const ModifiableLinearOp getDiagonalOp(
const LinearOp & op)
836 bool isTpetra =
false;
837 RCP<const Epetra_CrsMatrix> eCrsOp;
838 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
842 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
843 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
846 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
848 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
850 else if (!tOp.is_null()){
851 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
855 throw std::logic_error(
"Neither Epetra nor Tpetra");
857 catch (std::exception & e) {
858 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
860 *out <<
"Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
861 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
863 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
865 *out <<
"*** THROWN EXCEPTION ***\n";
866 *out << e.what() << std::endl;
867 *out <<
"************************\n";
874 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
875 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
878 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
882 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
883 tCrsOp->getLocalDiagCopy(*diag);
886 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
891 const MultiVector getDiagonal(
const LinearOp & op)
893 bool isTpetra =
false;
894 RCP<const Epetra_CrsMatrix> eCrsOp;
895 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
899 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
900 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
903 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
905 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
907 else if (!tOp.is_null()){
908 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
912 throw std::logic_error(
"Neither Epetra nor Tpetra");
914 catch (std::exception & e) {
915 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
917 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
918 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
922 *out <<
"Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
923 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
925 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
927 *out <<
"*** THROWN EXCEPTION ***\n";
928 *out << e.what() << std::endl;
929 *out <<
"************************\n";
936 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
937 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
939 return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
943 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
944 tCrsOp->getLocalDiagCopy(*diag);
947 return Thyra::createVector<ST,LO,GO,NT>(diag,Thyra::createVectorSpace<ST,LO,GO,NT>(tCrsOp->getRowMap()));
952 const MultiVector getDiagonal(
const Teko::LinearOp & A,
const DiagonalType & dt)
954 LinearOp diagOp = Teko::getDiagonalOp(A,dt);
956 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
957 Teuchos::rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
958 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
972 const ModifiableLinearOp getInvDiagonalOp(
const LinearOp & op)
974 if (Teko::TpetraHelpers::isTpetraLinearOp(op)){
977 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
980 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
982 tCrsOp->getLocalDiagCopy(*diag);
983 diag->reciprocal(*diag);
986 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
990 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op,
true);
991 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
994 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
995 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
996 diag->Reciprocal(*diag);
999 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
1015 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr)
1017 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1018 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1019 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1021 if(isTpetral && isTpetram && isTpetrar){
1025 bool transpl =
false;
1026 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1028 bool transpm =
false;
1029 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1031 bool transpr =
false;
1032 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1035 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1036 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1039 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1040 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1041 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1042 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1043 explicitCrsOp->resumeFill();
1044 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1045 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1046 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1049 }
else if (isTpetral && !isTpetram && isTpetrar){
1053 bool transpl =
false;
1054 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1056 bool transpr =
false;
1057 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1060 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opm,
true);
1061 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1062 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1063 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1066 tCrsOplm->rightScale(*diagPtr);
1069 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1070 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1073 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1074 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1075 explicitCrsOp->resumeFill();
1076 explicitCrsOp->scale(scalarl*scalarr);
1077 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1078 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1084 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1087 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1088 Thyra::epetraExtDiagScaledMatProdTransformer();
1091 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1092 prodTrans->transform(*implicitOp,explicitOp.ptr());
1093 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1094 " * " + opm->getObjectLabel() +
1095 " * " + opr->getObjectLabel() +
" )");
1116 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr,
1117 const ModifiableLinearOp & destOp)
1119 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1120 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1121 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1123 if(isTpetral && isTpetram && isTpetrar){
1127 bool transpl =
false;
1128 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1130 bool transpm =
false;
1131 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1133 bool transpr =
false;
1134 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1137 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1138 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1141 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1142 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1143 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1144 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1145 explicitCrsOp->resumeFill();
1146 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1147 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1148 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1151 }
else if (isTpetral && !isTpetram && isTpetrar){
1155 bool transpl =
false;
1156 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1158 bool transpr =
false;
1159 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1162 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opm,
true);
1163 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1164 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1165 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1168 tCrsOplm->rightScale(*diagPtr);
1171 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1172 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1175 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1176 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1177 explicitCrsOp->resumeFill();
1178 explicitCrsOp->scale(scalarl*scalarr);
1179 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1180 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1186 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1189 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1190 Thyra::epetraExtDiagScaledMatProdTransformer();
1193 ModifiableLinearOp explicitOp;
1196 if(destOp==Teuchos::null)
1197 explicitOp = prodTrans->createOutputOp();
1199 explicitOp = destOp;
1202 prodTrans->transform(*implicitOp,explicitOp.ptr());
1205 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1206 " * " + opm->getObjectLabel() +
1207 " * " + opr->getObjectLabel() +
" )");
1224 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr)
1226 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1227 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1229 if(isTpetral && isTpetrar){
1232 bool transpl =
false;
1233 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1235 bool transpr =
false;
1236 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1239 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1240 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1243 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1244 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1245 explicitCrsOp->resumeFill();
1246 explicitCrsOp->scale(scalarl*scalarr);
1247 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1248 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1251 }
else if (isTpetral && !isTpetrar){
1255 bool transpl =
false;
1256 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1259 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opr);
1260 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1261 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1262 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1264 explicitCrsOp->rightScale(*diagPtr);
1265 explicitCrsOp->resumeFill();
1266 explicitCrsOp->scale(scalarl);
1267 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1269 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1271 }
else if (!isTpetral && isTpetrar){
1275 bool transpr =
false;
1276 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1279 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opl);
1280 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1281 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1282 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpr, Tpetra::Import<LO,GO,NT>(tCrsOpr->getRowMap(),tCrsOpr->getRowMap()));
1284 explicitCrsOp->leftScale(*diagPtr);
1285 explicitCrsOp->resumeFill();
1286 explicitCrsOp->scale(scalarr);
1287 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1289 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1294 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1297 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1298 = Thyra::epetraExtDiagScalingTransformer();
1302 if(not prodTrans->isCompatible(*implicitOp))
1303 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1306 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1307 prodTrans->transform(*implicitOp,explicitOp.ptr());
1308 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1309 " * " + opr->getObjectLabel() +
" )");
1328 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr,
1329 const ModifiableLinearOp & destOp)
1331 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1332 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1334 if(isTpetral && isTpetrar){
1338 bool transpl =
false;
1339 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1341 bool transpr =
false;
1342 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1345 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1346 if(destOp!=Teuchos::null)
1347 explicitOp = destOp;
1349 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1350 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1353 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1354 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1355 explicitCrsOp->resumeFill();
1356 explicitCrsOp->scale(scalarl*scalarr);
1357 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1358 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1361 }
else if (isTpetral && !isTpetrar){
1365 bool transpl =
false;
1366 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1369 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opr);
1370 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1371 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1374 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1375 explicitCrsOp->rightScale(*diagPtr);
1376 explicitCrsOp->resumeFill();
1377 explicitCrsOp->scale(scalarl);
1378 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1379 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1381 }
else if (!isTpetral && isTpetrar){
1385 bool transpr =
false;
1386 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1389 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opl,
true);
1390 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1391 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1394 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpr, Tpetra::Import<LO,GO,NT>(tCrsOpr->getRowMap(),tCrsOpr->getRowMap()));
1395 explicitCrsOp->leftScale(*diagPtr);
1396 explicitCrsOp->resumeFill();
1397 explicitCrsOp->scale(scalarr);
1398 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1399 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1404 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1408 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1409 = Thyra::epetraExtDiagScalingTransformer();
1413 if(not prodTrans->isCompatible(*implicitOp))
1414 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1417 ModifiableLinearOp explicitOp;
1420 if(destOp==Teuchos::null)
1421 explicitOp = prodTrans->createOutputOp();
1423 explicitOp = destOp;
1426 prodTrans->transform(*implicitOp,explicitOp.ptr());
1429 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1430 " * " + opr->getObjectLabel() +
" )");
1446 const LinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr)
1448 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1449 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1451 if(isTpetral && isTpetrar){
1455 bool transpl =
false;
1456 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1458 bool transpr =
false;
1459 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1462 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1463 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1466 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1467 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1472 const LinearOp implicitOp = Thyra::add(opl,opr);
1475 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1476 Thyra::epetraExtAddTransformer();
1479 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1480 prodTrans->transform(*implicitOp,explicitOp.ptr());
1481 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1482 " + " + opr->getObjectLabel() +
" )");
1500 const ModifiableLinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr,
1501 const ModifiableLinearOp & destOp)
1503 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1504 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1506 if(isTpetral && isTpetrar){
1510 bool transpl =
false;
1511 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1513 bool transpr =
false;
1514 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1517 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1518 if(destOp!=Teuchos::null)
1519 explicitOp = destOp;
1521 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1522 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1525 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1526 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1532 const LinearOp implicitOp = Thyra::add(opl,opr);
1535 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1536 Thyra::epetraExtAddTransformer();
1539 RCP<Thyra::LinearOpBase<double> > explicitOp;
1540 if(destOp!=Teuchos::null)
1541 explicitOp = destOp;
1543 explicitOp = prodTrans->createOutputOp();
1546 prodTrans->transform(*implicitOp,explicitOp.ptr());
1547 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1548 " + " + opr->getObjectLabel() +
" )");
1558 const ModifiableLinearOp explicitSum(
const LinearOp & op,
1559 const ModifiableLinearOp & destOp)
1562 const RCP<const Epetra_CrsMatrix> epetraOp =
1563 rcp_dynamic_cast<
const Epetra_CrsMatrix>(get_Epetra_Operator(*op),
true);
1565 if(destOp==Teuchos::null) {
1566 Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(
new Epetra_CrsMatrix(*epetraOp));
1568 return Thyra::nonconstEpetraLinearOp(epetraDest);
1571 const RCP<Epetra_CrsMatrix> epetraDest =
1572 rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp),
true);
1574 EpetraExt::MatrixMatrix::Add(*epetraOp,
false,1.0,*epetraDest,1.0);
1579 const LinearOp explicitTranspose(
const LinearOp & op)
1581 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
1583 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op,
true);
1584 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
1586 Tpetra::RowMatrixTransposer<ST,LO,GO,NT> transposer(tCrsOp);
1587 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > transOp = transposer.createTranspose();
1589 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getDomainMap()),transOp);
1593 RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1594 TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
1595 "Teko::explicitTranspose Not an Epetra_Operator");
1596 RCP<const Epetra_RowMatrix> eRowMatrixOp
1597 = Teuchos::rcp_dynamic_cast<
const Epetra_RowMatrix>(eOp);
1598 TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
1599 "Teko::explicitTranspose Not an Epetra_RowMatrix");
1602 EpetraExt::RowMatrix_Transpose tranposeOp;
1603 Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
1607 Teuchos::RCP<Epetra_CrsMatrix> crsMat
1608 = Teuchos::rcp(
new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
1610 return Thyra::epetraLinearOp(crsMat);
1614 const double frobeniusNorm(
const LinearOp & op)
1616 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
1617 const RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
1618 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
1619 return crsOp->getFrobeniusNorm();
1621 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
1622 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(epOp,
true);
1623 return crsOp->NormFrobenius();
1627 const LinearOp buildDiagonal(
const MultiVector & src,
const std::string & lbl)
1629 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
1630 Thyra::copy(*src->col(0),dst.ptr());
1632 return Thyra::diagonal<double>(dst,lbl);
1635 const LinearOp buildInvDiagonal(
const MultiVector & src,
const std::string & lbl)
1637 const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
1638 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
1639 Thyra::reciprocal<double>(*srcV,dst.ptr());
1641 return Thyra::diagonal<double>(dst,lbl);
1645 BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> & mvv)
1647 Teuchos::Array<MultiVector> mvA;
1648 Teuchos::Array<VectorSpace> vsA;
1651 std::vector<MultiVector>::const_iterator itr;
1652 for(itr=mvv.begin();itr!=mvv.end();++itr) {
1653 mvA.push_back(*itr);
1654 vsA.push_back((*itr)->range());
1658 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
1659 = Thyra::productVectorSpace<double>(vsA);
1661 return Thyra::defaultProductMultiVector<double>(vs,mvA);
1674 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
1675 const std::vector<int> & indices,
1676 const VectorSpace & vs,
1677 double onValue,
double offValue)
1683 RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
1684 Thyra::put_scalar<double>(offValue,v.ptr());
1687 for(std::size_t i=0;i<indices.size();i++)
1688 Thyra::set_ele<double>(indices[i],onValue,v.ptr());
1717 double computeSpectralRad(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
1718 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
1720 typedef Thyra::LinearOpBase<double> OP;
1721 typedef Thyra::MultiVectorBase<double> MV;
1723 int startVectors = 1;
1726 const RCP<MV> ivec = Thyra::createMember(A->domain());
1727 Thyra::randomize(-1.0,1.0,ivec.ptr());
1729 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1730 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1732 eigProb->setHermitian(isHermitian);
1735 if(not eigProb->setProblem()) {
1737 return Teuchos::ScalarTraits<double>::nan();
1741 std::string which(
"LM");
1745 Teuchos::ParameterList MyPL;
1746 MyPL.set(
"Verbosity", verbosity );
1747 MyPL.set(
"Which", which );
1748 MyPL.set(
"Block Size", startVectors );
1749 MyPL.set(
"Num Blocks", numBlocks );
1750 MyPL.set(
"Maximum Restarts", restart );
1751 MyPL.set(
"Convergence Tolerance", tol );
1759 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1762 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1764 if(solverreturn==Anasazi::Unconverged) {
1765 double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
1766 double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
1768 return -std::abs(std::sqrt(real*real+comp*comp));
1774 double real = eigProb->getSolution().Evals.begin()->realpart;
1775 double comp = eigProb->getSolution().Evals.begin()->imagpart;
1777 return std::abs(std::sqrt(real*real+comp*comp));
1807 double computeSmallestMagEig(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
1808 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
1810 typedef Thyra::LinearOpBase<double> OP;
1811 typedef Thyra::MultiVectorBase<double> MV;
1813 int startVectors = 1;
1816 const RCP<MV> ivec = Thyra::createMember(A->domain());
1817 Thyra::randomize(-1.0,1.0,ivec.ptr());
1819 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1820 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1822 eigProb->setHermitian(isHermitian);
1825 if(not eigProb->setProblem()) {
1827 return Teuchos::ScalarTraits<double>::nan();
1831 std::string which(
"SM");
1834 Teuchos::ParameterList MyPL;
1835 MyPL.set(
"Verbosity", verbosity );
1836 MyPL.set(
"Which", which );
1837 MyPL.set(
"Block Size", startVectors );
1838 MyPL.set(
"Num Blocks", numBlocks );
1839 MyPL.set(
"Maximum Restarts", restart );
1840 MyPL.set(
"Convergence Tolerance", tol );
1848 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1851 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1853 if(solverreturn==Anasazi::Unconverged) {
1855 return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1859 return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1871 ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt)
1875 return getDiagonalOp(A);
1877 return getLumpedMatrix(A);
1879 return getAbsRowSumMatrix(A);
1882 TEUCHOS_TEST_FOR_EXCEPT(
true);
1885 return Teuchos::null;
1896 ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp & A,
const Teko::DiagonalType & dt)
1900 return getInvDiagonalOp(A);
1902 return getInvLumpedMatrix(A);
1904 return getAbsRowSumInvMatrix(A);
1907 TEUCHOS_TEST_FOR_EXCEPT(
true);
1910 return Teuchos::null;
1919 std::string getDiagonalName(
const DiagonalType & dt)
1947 if(name==
"Diagonal")
1951 if(name==
"AbsRowSum")
1959 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp & Op){
1960 #ifdef Teko_ENABLE_Isorropia 1961 Teuchos::ParameterList probeList;
1962 Prober prober(G,probeList,
true);
1963 Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(
new Epetra_CrsMatrix(Copy,*G));
1965 prober.probe(Mwrap,*Mat);
1966 return Thyra::epetraLinearOp(Mat);
1968 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
"Probe requires Isorropia");
1972 double norm_1(
const MultiVector & v,std::size_t col)
1974 Teuchos::Array<double> n(v->domain()->dim());
1975 Thyra::norms_1<double>(*v,n);
1980 double norm_2(
const MultiVector & v,std::size_t col)
1982 Teuchos::Array<double> n(v->domain()->dim());
1983 Thyra::norms_2<double>(*v,n);
1988 void putScalar(
const ModifiableLinearOp & op,
double scalar)
1992 RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1995 RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
1997 eCrsOp->PutScalar(scalar);
1999 catch (std::exception & e) {
2000 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2002 *out <<
"Teko: putScalar requires an Epetra_CrsMatrix\n";
2003 *out <<
" Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2005 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2007 *out <<
"*** THROWN EXCEPTION ***\n";
2008 *out << e.what() << std::endl;
2009 *out <<
"************************\n";
2015 void clipLower(MultiVector & v,
double lowerBound)
2018 using Teuchos::rcp_dynamic_cast;
2024 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2025 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2026 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2028 Teuchos::ArrayRCP<double> values;
2030 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2031 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2032 if(values[j]<lowerBound)
2033 values[j] = lowerBound;
2038 void clipUpper(MultiVector & v,
double upperBound)
2041 using Teuchos::rcp_dynamic_cast;
2046 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2047 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2048 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2050 Teuchos::ArrayRCP<double> values;
2052 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2053 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2054 if(values[j]>upperBound)
2055 values[j] = upperBound;
2060 void replaceValue(MultiVector & v,
double currentValue,
double newValue)
2063 using Teuchos::rcp_dynamic_cast;
2068 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2069 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2070 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2072 Teuchos::ArrayRCP<double> values;
2074 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2075 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2076 if(values[j]==currentValue)
2077 values[j] = newValue;
2082 void columnAverages(
const MultiVector & v,std::vector<double> & averages)
2084 averages.resize(v->domain()->dim());
2087 Thyra::sums<double>(*v,averages);
2090 Thyra::Ordinal rows = v->range()->dim();
2091 for(std::size_t i=0;i<averages.size();i++)
2092 averages[i] = averages[i]/rows;
2095 double average(
const MultiVector & v)
2097 Thyra::Ordinal rows = v->range()->dim();
2098 Thyra::Ordinal cols = v->domain()->dim();
2100 std::vector<double> averages;
2101 columnAverages(v,averages);
2104 for(std::size_t i=0;i<averages.size();i++)
2105 sum += averages[i] * rows;
2107 return sum/(rows*cols);
Specifies that a block diagonal approximation is to be used.
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
Specifies that row sum is used to form a diagonal.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...
For user convenience, if Teko recieves this value, exceptions will be thrown.
DiagonalType
Type describing the type of diagonal to construct.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
Specifies that the diagonal entry is .