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_DefaultScaledAdjointLinearOp.hpp"
56#include "Thyra_DefaultMultipliedLinearOp.hpp"
57#include "Thyra_DefaultZeroLinearOp.hpp"
58#include "Thyra_DefaultProductMultiVector.hpp"
59#include "Thyra_DefaultProductVectorSpace.hpp"
60#include "Thyra_MultiVectorStdOps.hpp"
61#include "Thyra_VectorStdOps.hpp"
62#include "Thyra_SpmdVectorBase.hpp"
65#ifdef TEKO_HAVE_EPETRA
66#include "Thyra_EpetraExtDiagScaledMatProdTransformer.hpp"
67#include "Thyra_EpetraExtDiagScalingTransformer.hpp"
68#include "Thyra_EpetraExtAddTransformer.hpp"
69#include "Thyra_get_Epetra_Operator.hpp"
70#include "Thyra_EpetraThyraWrappers.hpp"
71#include "Thyra_EpetraOperatorWrapper.hpp"
72#include "Thyra_EpetraLinearOp.hpp"
76#include "Teuchos_Array.hpp"
79#ifdef TEKO_HAVE_EPETRA
80#include "Epetra_Operator.h"
81#include "Epetra_CrsGraph.h"
82#include "Epetra_CrsMatrix.h"
83#include "Epetra_Vector.h"
84#include "Epetra_Map.h"
86#include "EpetraExt_Transpose_RowMatrix.h"
87#include "EpetraExt_MatrixMatrix.h"
89#include "Teko_EpetraHelpers.hpp"
90#include "Teko_EpetraOperatorWrapper.hpp"
94#include "AnasaziBasicEigenproblem.hpp"
95#include "AnasaziThyraAdapter.hpp"
96#include "AnasaziBlockKrylovSchurSolMgr.hpp"
97#include "AnasaziBlockKrylovSchur.hpp"
98#include "AnasaziStatusTestMaxIters.hpp"
101#ifdef Teko_ENABLE_Isorropia
102#include "Isorropia_EpetraProber.hpp"
106#include "Teko_TpetraHelpers.hpp"
107#include "Teko_TpetraOperatorWrapper.hpp"
110#include "Thyra_TpetraLinearOp.hpp"
111#include "Tpetra_CrsMatrix.hpp"
112#include "Tpetra_Vector.hpp"
113#include "Thyra_TpetraThyraWrappers.hpp"
114#include "TpetraExt_MatrixMatrix.hpp"
115#include "Tpetra_RowMatrixTransposer.hpp"
122using Teuchos::rcp_dynamic_cast;
124#ifdef Teko_ENABLE_Isorropia
125using Isorropia::Epetra::Prober;
128const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
130 Teuchos::RCP<Teuchos::FancyOStream> os =
131 Teuchos::VerboseObjectBase::getDefaultOStream();
139inline double dist(
int dim,
double * coords,
int row,
int col)
142 for(
int i=0;i<dim;i++)
143 value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
146 return std::sqrt(value);
150inline double dist(
double * x,
double * y,
double * z,
int stride,
int row,
int col)
153 if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
154 if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
155 if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
158 return std::sqrt(value);
179#ifdef TEKO_HAVE_EPETRA
180RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double * coords,
const Epetra_CrsMatrix & stencil)
183 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
184 stencil.MaxNumEntries()+1),
true);
187 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
188 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
191 for(
int j=0;j<gl->NumMyRows();j++) {
192 int row = gl->GRID(j);
193 double diagValue = 0.0;
198 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
201 for(
int i=0;i<rowSz;i++) {
205 double value = rowData[i];
206 if(value==0)
continue;
210 double d = dist(dim,coords,row,col);
212 diagValue += rowData[i];
220 rowData[rowSz] = -diagValue;
225 rowData[diagInd] = -diagValue;
226 rowInd[diagInd] = row;
230 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
239RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(
int dim,ST * coords,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
242 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
243 stencil.getGlobalMaxNumRowEntries()+1));
246 auto rowInd =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
247 auto rowData =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
250 for(LO j=0;j< (LO) gl->getLocalNumRows();j++) {
251 GO row = gl->getRowMap()->getGlobalElement(j);
257 stencil.getGlobalRowCopy(row,rowInd,rowData,rowSz);
260 for(
size_t i=0;i<rowSz;i++) {
264 ST value = rowData[i];
265 if(value==0)
continue;
269 ST d = dist(dim,coords,row,col);
271 diagValue += rowData(i);
279 rowData(rowSz) = -diagValue;
284 rowData(diagInd) = -diagValue;
285 rowInd(diagInd) = row;
289 gl->replaceGlobalValues(row,rowInd,rowData);
319#ifdef TEKO_HAVE_EPETRA
320RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double * x,
double * y,
double * z,
int stride,
const Epetra_CrsMatrix & stencil)
323 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
324 stencil.MaxNumEntries()+1),
true);
327 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
328 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
331 for(
int j=0;j<gl->NumMyRows();j++) {
332 int row = gl->GRID(j);
333 double diagValue = 0.0;
338 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
341 for(
int i=0;i<rowSz;i++) {
345 double value = rowData[i];
346 if(value==0)
continue;
350 double d = dist(x,y,z,stride,row,col);
352 diagValue += rowData[i];
360 rowData[rowSz] = -diagValue;
365 rowData[diagInd] = -diagValue;
366 rowInd[diagInd] = row;
370 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
379RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
382 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
383 stencil.getGlobalMaxNumRowEntries()+1),
true);
386 auto rowInd =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
387 auto rowData =
typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing(
"rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
390 for(LO j=0;j<(LO) gl->getLocalNumRows();j++) {
391 GO row = gl->getRowMap()->getGlobalElement(j);
397 stencil.getGlobalRowCopy(row,rowInd,rowData,rowSz);
400 for(
size_t i=0;i<rowSz;i++) {
404 ST value = rowData[i];
405 if(value==0)
continue;
409 ST d = dist(x,y,z,stride,row,col);
411 diagValue += rowData(i);
419 rowData(rowSz) = -diagValue;
424 rowData(diagInd) = -diagValue;
425 rowInd(diagInd) = row;
429 gl->replaceGlobalValues(row,rowInd,rowData);
452void applyOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha,
double beta)
454 Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
472void applyTransposeOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha,
double beta)
474 Thyra::apply(*A,Thyra::TRANS,*x,y.ptr(),alpha,beta);
479void update(
double alpha,
const MultiVector & x,
double beta,MultiVector & y)
481 Teuchos::Array<double> scale;
482 Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
485 scale.push_back(alpha);
486 vec.push_back(x.ptr());
489 Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
493BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
499 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
500 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
506 upper->beginBlockFill(rows,rows);
508 for(
int i=0;i<rows;i++) {
512 RCP<const Thyra::LinearOpBase<double> > zed
513 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
514 upper->setBlock(i,i,zed);
516 for(
int j=i+1;j<rows;j++) {
518 LinearOp uij = blo->getBlock(i,j);
521 if(uij!=Teuchos::null)
522 upper->setBlock(i,j,uij);
526 upper->endBlockFill();
532BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
538 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
539 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
545 lower->beginBlockFill(rows,rows);
547 for(
int i=0;i<rows;i++) {
551 RCP<const Thyra::LinearOpBase<double> > zed
552 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
553 lower->setBlock(i,i,zed);
555 for(
int j=0;j<i;j++) {
557 LinearOp uij = blo->getBlock(i,j);
560 if(uij!=Teuchos::null)
561 lower->setBlock(i,j,uij);
565 lower->endBlockFill();
589BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp & blo)
595 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
596 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
602 zeroOp->beginBlockFill(rows,rows);
604 for(
int i=0;i<rows;i++) {
608 RCP<const Thyra::LinearOpBase<double> > zed
609 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
610 zeroOp->setBlock(i,i,zed);
617bool isZeroOp(
const LinearOp op)
620 if(op==Teuchos::null)
return true;
623 LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
626 if(test!=Teuchos::null)
return true;
630 Thyra::EOpTransp transp = Thyra::NOTRANS;
631 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
632 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
633 test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(wrapped_op);
634 return test!=Teuchos::null;
637std::pair<ModifiableLinearOp, bool> getAbsRowSumMatrixEpetra(
const LinearOp &op) {
638#ifndef TEKO_HAVE_EPETRA
639 return std::make_pair(ModifiableLinearOp{},
false);
641 RCP<const Epetra_CrsMatrix> eCrsOp;
643 const auto eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op);
646 return std::make_pair(ModifiableLinearOp{},
false);
649 eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
653 const auto ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
654 Epetra_Vector &diag = *ptrDiag;
658 for (
int i = 0; i < eCrsOp->NumMyRows(); i++) {
661 eCrsOp->ExtractMyRowView(i, numEntries, values);
664 for (
int j = 0; j < numEntries; j++)
665 diag[i] += std::abs(values[j]);
669 return std::make_pair(
670 Teko::Epetra::thyraDiagOp(ptrDiag, eCrsOp->RowMap(),
671 "absRowSum( " + op->getObjectLabel() +
" )"),
676std::pair<ModifiableLinearOp, bool> getAbsRowSumMatrixTpetra(
const LinearOp &op) {
677 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT>> tCrsOp;
680 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT>>(op);
682 tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT>>(
683 tOp->getConstTpetraOperator(),
true);
687 Tpetra::createVector<ST, LO, GO, NT>(tCrsOp->getRowMap());
688 auto &diag = *ptrDiag;
692 for (LO i = 0; i < (LO)tCrsOp->getLocalNumRows(); i++) {
693 auto numEntries = tCrsOp->getNumEntriesInLocalRow(i);
694 typename Tpetra::CrsMatrix<ST, LO, GO, NT>::local_inds_host_view_type
696 typename Tpetra::CrsMatrix<ST, LO, GO, NT>::values_host_view_type values;
697 tCrsOp->getLocalRowView(i, indices, values);
700 for (
size_t j = 0; j < numEntries; j++)
701 diag.sumIntoLocalValue(i, std::abs(values(j)));
705 return std::make_pair(Teko::TpetraHelpers::thyraDiagOp(
706 ptrDiag, *tCrsOp->getRowMap(),
707 "absRowSum( " + op->getObjectLabel() +
" ))"),
719ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp & op)
722 auto eResult = getAbsRowSumMatrixEpetra(op);
723 if (eResult.second) {
724 return eResult.first;
727 auto tResult = getAbsRowSumMatrixTpetra(op);
728 if (tResult.second) {
729 return tResult.first;
731 throw std::logic_error(
"Neither Epetra nor Tpetra");
733 }
catch (std::exception &e) {
734 auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
736 *out <<
"Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a "
737 "Tpetra::CrsMatrix\n";
738 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator "
740 << op->description() << std::endl;
742 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or "
743 "a Tpetra_Operator to a Tpetra::CrsMatrix\n";
745 *out <<
"*** THROWN EXCEPTION ***\n";
746 *out << e.what() << std::endl;
747 *out <<
"************************\n";
761ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp & op)
765 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
766 if(blocked_op != Teuchos::null){
767 int numRows = blocked_op->productRange()->numBlocks();
768 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == numRows);
769 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_diag = Thyra::defaultBlockedLinearOp<double>();
770 blocked_diag->beginBlockFill(numRows,numRows);
771 for(
int r = 0; r < numRows; ++r){
772 for(
int c = 0; c < numRows; ++c){
774 blocked_diag->setNonconstBlock(r,c,getAbsRowSumInvMatrix(blocked_op->getBlock(r,c)));
776 blocked_diag->setBlock(r,c,Thyra::zero<double>(blocked_op->getBlock(r,c)->range(),blocked_op->getBlock(r,c)->domain()));
779 blocked_diag->endBlockFill();
783 if(Teko::TpetraHelpers::isTpetraLinearOp(op)) {
786 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
789 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
790 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
794 for(LO i=0;i<(LO) tCrsOp->getLocalNumRows();i++) {
795 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
796 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::local_inds_host_view_type indices;
797 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::values_host_view_type values;
798 tCrsOp->getLocalRowView(i,indices,values);
801 for(LO j=0;j<numEntries;j++)
802 diag.sumIntoLocalValue(i,std::abs(values(j)));
805 diag.reciprocal(diag);
808 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),
"absRowSum( " + op->getObjectLabel() +
" ))");
812#ifdef TEKO_HAVE_EPETRA
813 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,
true);
814 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
817 const RCP<Epetra_Vector> ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
818 Epetra_Vector & diag = *ptrDiag;
822 for(
int i=0;i<eCrsOp->NumMyRows();i++) {
825 eCrsOp->ExtractMyRowView(i,numEntries,values);
828 for(
int j=0;j<numEntries;j++)
829 diag[i] += std::abs(values[j]);
831 diag.Reciprocal(diag);
834 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),
"absRowSum( " + op->getObjectLabel() +
" )");
836 throw std::logic_error(
"getAbsRowSumInvMatrix is trying to use Epetra "
837 "code, but TEKO_HAVE_EPETRA is disabled!");
850ModifiableLinearOp getLumpedMatrix(
const LinearOp & op)
852 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
853 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
856 Thyra::assign(ones.ptr(),1.0);
860 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
862 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
873ModifiableLinearOp getInvLumpedMatrix(
const LinearOp & op)
875 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
876 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
879 Thyra::assign(ones.ptr(),1.0);
882 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
883 Thyra::reciprocal(*diag,diag.ptr());
885 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
888const std::pair<ModifiableLinearOp, bool> getDiagonalOpEpetra(
const LinearOp &op) {
889#ifndef TEKO_HAVE_EPETRA
890 return std::make_pair(ModifiableLinearOp{},
false);
892 RCP<const Epetra_CrsMatrix> eCrsOp;
894 const auto eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op);
896 return std::make_pair(ModifiableLinearOp{},
false);
899 eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
902 const auto diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
903 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
906 return std::make_pair(
907 Teko::Epetra::thyraDiagOp(diag, eCrsOp->RowMap(),
908 "inv(diag( " + op->getObjectLabel() +
" ))"),
913const std::pair<ModifiableLinearOp, bool> getDiagonalOpTpetra(
const LinearOp &op) {
914 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT>> tCrsOp;
917 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT>>(op);
919 return std::make_pair(ModifiableLinearOp{},
false);
922 tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT>>(
923 tOp->getConstTpetraOperator(),
true);
926 const auto diag = Tpetra::createVector<ST, LO, GO, NT>(tCrsOp->getRowMap());
927 tCrsOp->getLocalDiagCopy(*diag);
930 return std::make_pair(Teko::TpetraHelpers::thyraDiagOp(
931 diag, *tCrsOp->getRowMap(),
932 "inv(diag( " + op->getObjectLabel() +
" ))"),
947const ModifiableLinearOp getDiagonalOp(
const LinearOp & op)
951 const auto eDiagOp = getDiagonalOpEpetra(op);
953 if (eDiagOp.second) {
954 return eDiagOp.first;
957 const auto tDiagOp = getDiagonalOpTpetra(op);
958 if (tDiagOp.second) {
959 return tDiagOp.first;
962 throw std::logic_error(
"Neither Epetra nor Tpetra");
964 catch (std::exception & e) {
965 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
967 *out <<
"Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
968 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
970 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
972 *out <<
"*** THROWN EXCEPTION ***\n";
973 *out << e.what() << std::endl;
974 *out <<
"************************\n";
980const MultiVector getDiagonal(
const LinearOp & op)
984 auto diagOp = getDiagonalOpEpetra(op);
986 if (!diagOp.second) {
987 diagOp = getDiagonalOpTpetra(op);
990 throw std::logic_error(
"Neither Epetra nor Tpetra");
994 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
995 Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp.first)->getDiag();
996 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
998 catch (std::exception & e) {
999 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1001 *out <<
"Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
1002 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
1004 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
1006 *out <<
"*** THROWN EXCEPTION ***\n";
1007 *out << e.what() << std::endl;
1008 *out <<
"************************\n";
1014const MultiVector getDiagonal(
const Teko::LinearOp & A,
const DiagonalType & dt)
1016 LinearOp diagOp = Teko::getDiagonalOp(A,dt);
1018 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
1019 Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
1020 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
1034const ModifiableLinearOp getInvDiagonalOp(
const LinearOp & op)
1037 auto diagonal_op = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double>>(op);
1038 if(diagonal_op != Teuchos::null){
1039 auto diag = diagonal_op->getDiag();
1040 auto inv_diag = diag->clone_v();
1041 Thyra::reciprocal(*diag,inv_diag.ptr());
1042 return rcp(
new Thyra::DefaultDiagonalLinearOp<double>(inv_diag));
1046 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
1047 if(blocked_op != Teuchos::null){
1048 int numRows = blocked_op->productRange()->numBlocks();
1049 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == numRows);
1050 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_diag = Thyra::defaultBlockedLinearOp<double>();
1051 blocked_diag->beginBlockFill(numRows,numRows);
1052 for(
int r = 0; r < numRows; ++r){
1053 for(
int c = 0; c < numRows; ++c){
1055 blocked_diag->setNonconstBlock(r,c,getInvDiagonalOp(blocked_op->getBlock(r,c)));
1057 blocked_diag->setBlock(r,c,Thyra::zero<double>(blocked_op->getBlock(r,c)->range(),blocked_op->getBlock(r,c)->domain()));
1060 blocked_diag->endBlockFill();
1061 return blocked_diag;
1064 if (Teko::TpetraHelpers::isTpetraLinearOp(op)){
1066 bool transp =
false;
1067 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
1070 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
1071 diag->scale(scalar);
1072 tCrsOp->getLocalDiagCopy(*diag);
1073 diag->reciprocal(*diag);
1076 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
1080#ifdef TEKO_HAVE_EPETRA
1081 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,
true);
1082 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
1085 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
1086 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
1087 diag->Reciprocal(*diag);
1090 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
1092 throw std::logic_error(
"getInvDiagonalOp is trying to use Epetra "
1093 "code, but TEKO_HAVE_EPETRA is disabled!");
1110const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr)
1115 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1116 bool isBlockedM = isPhysicallyBlockedLinearOp(opm);
1117 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1120 if((isBlockedL && isBlockedM && isBlockedR)){
1122 double scalarl = 0.0;
1123 bool transpl =
false;
1124 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1125 double scalarm = 0.0;
1126 bool transpm =
false;
1127 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opm = getPhysicallyBlockedLinearOp(opm,&scalarm,&transpm);
1128 double scalarr = 0.0;
1129 bool transpr =
false;
1130 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1131 double scalar = scalarl*scalarm*scalarr;
1133 int numRows = blocked_opl->productRange()->numBlocks();
1134 int numCols = blocked_opr->productDomain()->numBlocks();
1135 int numMiddle = blocked_opm->productRange()->numBlocks();
1138 TEUCHOS_ASSERT(blocked_opm->productDomain()->numBlocks() == numMiddle);
1139 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1140 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1142 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1143 blocked_product->beginBlockFill(numRows,numCols);
1144 for(
int r = 0; r < numRows; ++r){
1145 for(
int c = 0; c < numCols; ++c){
1146 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opm->getBlock(0,0),blocked_opr->getBlock(0,c));
1147 for(
int m = 1; m < numMiddle; ++m){
1148 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opm->getBlock(m,m),blocked_opr->getBlock(m,c));
1149 product_rc = explicitAdd(product_rc,product_m);
1151 blocked_product->setBlock(r,c,product_rc);
1154 blocked_product->endBlockFill();
1155 return Thyra::scale<double>(scalar,blocked_product.getConst());
1159 if(isBlockedL && !isBlockedM && isBlockedR){
1160 double scalarl = 0.0;
1161 bool transpl =
false;
1162 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1163 double scalarr = 0.0;
1164 bool transpr =
false;
1165 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1166 double scalar = scalarl*scalarr;
1168 int numRows = blocked_opl->productRange()->numBlocks();
1169 int numCols = blocked_opr->productDomain()->numBlocks();
1173 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1174 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1176 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1177 blocked_product->beginBlockFill(numRows,numCols);
1178 for(
int r = 0; r < numRows; ++r){
1179 for(
int c = 0; c < numCols; ++c){
1180 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),opm,blocked_opr->getBlock(0,c));
1181 blocked_product->setBlock(r,c,product_rc);
1184 blocked_product->endBlockFill();
1185 return Thyra::scale<double>(scalar,blocked_product.getConst());
1189 if(!isBlockedL && !isBlockedM && isBlockedR){
1190 double scalarr = 0.0;
1191 bool transpr =
false;
1192 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1193 double scalar = scalarr;
1196 int numCols = blocked_opr->productDomain()->numBlocks();
1200 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1202 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1203 blocked_product->beginBlockFill(numRows,numCols);
1204 for(
int c = 0; c < numCols; ++c){
1205 LinearOp product_c = explicitMultiply(opl,opm,blocked_opr->getBlock(0,c));
1206 blocked_product->setBlock(0,c,product_c);
1208 blocked_product->endBlockFill();
1209 return Thyra::scale<double>(scalar,blocked_product.getConst());
1214 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1215 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1216 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1218 if(isTpetral && isTpetram && isTpetrar){
1222 bool transpl =
false;
1223 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1225 bool transpm =
false;
1226 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1228 bool transpr =
false;
1229 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1232 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1233 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1236 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1237 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1238 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1239 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1240 explicitCrsOp->resumeFill();
1241 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1242 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1243 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1246 }
else if (isTpetral && !isTpetram && isTpetrar){
1250 bool transpl =
false;
1251 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1253 bool transpr =
false;
1254 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1256 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr;
1259 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm);
1260 if(dOpm != Teuchos::null){
1261 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1262 diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1265 else if(rcp_dynamic_cast<
const Thyra::ZeroLinearOpBase<ST> >(opm) != Teuchos::null){
1266 diagPtr = rcp(
new Tpetra::Vector<ST,LO,GO,NT>(tCrsOpl->getDomainMap()));
1269 TEUCHOS_ASSERT(
false);
1271 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()));
1274 tCrsOplm->rightScale(*diagPtr);
1277 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1278 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1281 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1282 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1283 explicitCrsOp->resumeFill();
1284 explicitCrsOp->scale(scalarl*scalarr);
1285 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1286 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1290#ifdef TEKO_HAVE_EPETRA
1292 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1295 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1296 Thyra::epetraExtDiagScaledMatProdTransformer();
1299 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1300 prodTrans->transform(*implicitOp,explicitOp.ptr());
1301 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1302 " * " + opm->getObjectLabel() +
1303 " * " + opr->getObjectLabel() +
" )");
1307 throw std::logic_error(
"explicitMultiply is trying to use Epetra "
1308 "code, but TEKO_HAVE_EPETRA is disabled!");
1327const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr,
1328 const ModifiableLinearOp & destOp)
1330 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1331 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1332 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1334 if(isTpetral && isTpetram && isTpetrar){
1338 bool transpl =
false;
1339 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1341 bool transpm =
false;
1342 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opm, &scalarm, &transpm);
1344 bool transpr =
false;
1345 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1348 auto tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(destOp);
1349 if(tExplicitOp.is_null())
1350 tExplicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1353 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1354 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1355 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1356 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1357 explicitCrsOp->resumeFill();
1358 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1359 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1360 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1363 }
else if (isTpetral && !isTpetram && isTpetrar){
1367 bool transpl =
false;
1368 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1370 bool transpr =
false;
1371 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1374 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,
true);
1375 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1376 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1377 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()));
1380 tCrsOplm->rightScale(*diagPtr);
1383 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1384 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1387 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1388 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1389 explicitCrsOp->resumeFill();
1390 explicitCrsOp->scale(scalarl*scalarr);
1391 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1392 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1396#ifdef TEKO_HAVE_EPETRA
1398 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1401 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1402 Thyra::epetraExtDiagScaledMatProdTransformer();
1405 ModifiableLinearOp explicitOp;
1408 if(destOp==Teuchos::null)
1409 explicitOp = prodTrans->createOutputOp();
1411 explicitOp = destOp;
1414 prodTrans->transform(*implicitOp,explicitOp.ptr());
1417 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1418 " * " + opm->getObjectLabel() +
1419 " * " + opr->getObjectLabel() +
" )");
1423 throw std::logic_error(
"explicitMultiply is trying to use Epetra "
1424 "code, but TEKO_HAVE_EPETRA is disabled!");
1439const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr)
1444 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1445 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1448 if((isBlockedL && isBlockedR)){
1449 double scalarl = 0.0;
1450 bool transpl =
false;
1451 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1452 double scalarr = 0.0;
1453 bool transpr =
false;
1454 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1455 double scalar = scalarl*scalarr;
1457 int numRows = blocked_opl->productRange()->numBlocks();
1458 int numCols = blocked_opr->productDomain()->numBlocks();
1459 int numMiddle = blocked_opl->productDomain()->numBlocks();
1461 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1463 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1464 blocked_product->beginBlockFill(numRows,numCols);
1465 for(
int r = 0; r < numRows; ++r){
1466 for(
int c = 0; c < numCols; ++c){
1467 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opr->getBlock(0,c));
1468 for(
int m = 1; m < numMiddle; ++m){
1469 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opr->getBlock(m,c));
1470 product_rc = explicitAdd(product_rc,product_m);
1472 blocked_product->setBlock(r,c,Thyra::scale(scalar,product_rc));
1475 blocked_product->endBlockFill();
1476 return blocked_product;
1480 if((isBlockedL && !isBlockedR)){
1481 double scalarl = 0.0;
1482 bool transpl =
false;
1483 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1484 double scalar = scalarl;
1486 int numRows = blocked_opl->productRange()->numBlocks();
1490 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1492 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1493 blocked_product->beginBlockFill(numRows,numCols);
1494 for(
int r = 0; r < numRows; ++r){
1495 LinearOp product_r = explicitMultiply(blocked_opl->getBlock(r,0),opr);
1496 blocked_product->setBlock(r,0,Thyra::scale(scalar,product_r));
1498 blocked_product->endBlockFill();
1499 return blocked_product;
1503 if((!isBlockedL && isBlockedR)){
1504 double scalarr = 0.0;
1505 bool transpr =
false;
1506 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1507 double scalar = scalarr;
1510 int numCols = blocked_opr->productDomain()->numBlocks();
1513 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1515 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1516 blocked_product->beginBlockFill(numRows,numCols);
1517 for(
int c = 0; c < numCols; ++c){
1518 LinearOp product_c = explicitMultiply(opl,blocked_opr->getBlock(0,c));
1519 blocked_product->setBlock(0,c,Thyra::scale(scalar,product_c));
1521 blocked_product->endBlockFill();
1522 return blocked_product;
1525 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1526 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1528 if(isTpetral && isTpetrar){
1531 bool transpl =
false;
1532 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1534 bool transpr =
false;
1535 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1538 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1539 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1542 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1543 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1544 explicitCrsOp->resumeFill();
1545 explicitCrsOp->scale(scalarl*scalarr);
1546 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1547 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1550 }
else if (isTpetral && !isTpetrar){
1554 bool transpl =
false;
1555 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1558 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr,
true);
1559 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1560 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1561 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()));
1563 explicitCrsOp->rightScale(*diagPtr);
1564 explicitCrsOp->resumeFill();
1565 explicitCrsOp->scale(scalarl);
1566 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1568 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1570 }
else if (!isTpetral && isTpetrar){
1574 bool transpr =
false;
1575 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1577 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr;
1580 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl);
1581 if(dOpl != Teuchos::null){
1582 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1583 diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1586 else if(rcp_dynamic_cast<
const Thyra::ZeroLinearOpBase<ST> >(opl) != Teuchos::null){
1587 diagPtr = rcp(
new Tpetra::Vector<ST,LO,GO,NT>(tCrsOpr->getRangeMap()));
1590 TEUCHOS_ASSERT(
false);
1592 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()));
1594 explicitCrsOp->leftScale(*diagPtr);
1595 explicitCrsOp->resumeFill();
1596 explicitCrsOp->scale(scalarr);
1597 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1599 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1602#ifdef TEKO_HAVE_EPETRA
1604 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1607 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1608 = Thyra::epetraExtDiagScalingTransformer();
1612 if(not prodTrans->isCompatible(*implicitOp))
1613 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1616 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1617 prodTrans->transform(*implicitOp,explicitOp.ptr());
1618 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1619 " * " + opr->getObjectLabel() +
" )");
1623 throw std::logic_error(
"explicitMultiply is trying to use Epetra "
1624 "code, but TEKO_HAVE_EPETRA is disabled!");
1642const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr,
1643 const ModifiableLinearOp & destOp)
1648 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1649 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1652 if((isBlockedL && isBlockedR)){
1653 double scalarl = 0.0;
1654 bool transpl =
false;
1655 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1656 double scalarr = 0.0;
1657 bool transpr =
false;
1658 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1659 double scalar = scalarl*scalarr;
1661 int numRows = blocked_opl->productRange()->numBlocks();
1662 int numCols = blocked_opr->productDomain()->numBlocks();
1663 int numMiddle = blocked_opl->productDomain()->numBlocks();
1665 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1667 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1668 blocked_product->beginBlockFill(numRows,numCols);
1669 for(
int r = 0; r < numRows; ++r){
1670 for(
int c = 0; c < numCols; ++c){
1672 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opr->getBlock(0,c));
1673 for(
int m = 1; m < numMiddle; ++m){
1674 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opr->getBlock(m,c));
1675 product_rc = explicitAdd(product_rc,product_m);
1677 blocked_product->setBlock(r,c,Thyra::scale(scalar,product_rc));
1680 blocked_product->endBlockFill();
1681 return blocked_product;
1685 if((isBlockedL && !isBlockedR)){
1686 double scalarl = 0.0;
1687 bool transpl =
false;
1688 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1689 double scalar = scalarl;
1691 int numRows = blocked_opl->productRange()->numBlocks();
1695 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1697 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1698 blocked_product->beginBlockFill(numRows,numCols);
1699 for(
int r = 0; r < numRows; ++r){
1700 LinearOp product_r = explicitMultiply(blocked_opl->getBlock(r,0),opr);
1701 blocked_product->setBlock(r,0,Thyra::scale(scalar,product_r));
1703 blocked_product->endBlockFill();
1704 return blocked_product;
1708 if((!isBlockedL && isBlockedR)){
1709 double scalarr = 0.0;
1710 bool transpr =
false;
1711 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1712 double scalar = scalarr;
1715 int numCols = blocked_opr->productDomain()->numBlocks();
1718 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1720 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1721 blocked_product->beginBlockFill(numRows,numCols);
1722 for(
int c = 0; c < numCols; ++c){
1723 LinearOp product_c = explicitMultiply(opl,blocked_opr->getBlock(0,c));
1724 blocked_product->setBlock(0,c,Thyra::scale(scalar,product_c));
1726 blocked_product->endBlockFill();
1727 return blocked_product;
1730 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1731 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1733 if(isTpetral && isTpetrar){
1737 bool transpl =
false;
1738 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1740 bool transpr =
false;
1741 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1744 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1745 if(destOp!=Teuchos::null)
1746 explicitOp = destOp;
1748 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1749 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1752 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1753 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1754 explicitCrsOp->resumeFill();
1755 explicitCrsOp->scale(scalarl*scalarr);
1756 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1757 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1760 }
else if (isTpetral && !isTpetrar){
1764 bool transpl =
false;
1765 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1768 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1769 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1770 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1773 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()));
1774 explicitCrsOp->rightScale(*diagPtr);
1775 explicitCrsOp->resumeFill();
1776 explicitCrsOp->scale(scalarl);
1777 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1778 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1780 }
else if (!isTpetral && isTpetrar){
1784 bool transpr =
false;
1785 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1788 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl,
true);
1789 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1790 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1793 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()));
1794 explicitCrsOp->leftScale(*diagPtr);
1795 explicitCrsOp->resumeFill();
1796 explicitCrsOp->scale(scalarr);
1797 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1798 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1801#ifdef TEKO_HAVE_EPETRA
1803 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1807 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1808 = Thyra::epetraExtDiagScalingTransformer();
1812 if(not prodTrans->isCompatible(*implicitOp))
1813 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1816 ModifiableLinearOp explicitOp;
1819 if(destOp==Teuchos::null)
1820 explicitOp = prodTrans->createOutputOp();
1822 explicitOp = destOp;
1825 prodTrans->transform(*implicitOp,explicitOp.ptr());
1828 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1829 " * " + opr->getObjectLabel() +
" )");
1833 throw std::logic_error(
"explicitMultiply is trying to use Epetra "
1834 "code, but TEKO_HAVE_EPETRA is disabled!");
1849const LinearOp explicitAdd(
const LinearOp & opl_in,
const LinearOp & opr_in)
1852 if(isPhysicallyBlockedLinearOp(opl_in) && isPhysicallyBlockedLinearOp(opr_in)){
1853 double scalarl = 0.0;
1854 bool transpl =
false;
1855 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl_in, &scalarl, &transpl);
1857 double scalarr = 0.0;
1858 bool transpr =
false;
1859 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr_in, &scalarr, &transpr);
1861 int numRows = blocked_opl->productRange()->numBlocks();
1862 int numCols = blocked_opl->productDomain()->numBlocks();
1863 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numRows);
1864 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == numCols);
1866 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_sum = Thyra::defaultBlockedLinearOp<double>();
1867 blocked_sum->beginBlockFill(numRows,numCols);
1868 for(
int r = 0; r < numRows; ++r)
1869 for(
int c = 0; c < numCols; ++c)
1870 blocked_sum->setBlock(r,c,explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c))));
1871 blocked_sum->endBlockFill();
1876 LinearOp opl = opl_in;
1877 LinearOp opr = opr_in;
1878 if(isPhysicallyBlockedLinearOp(opl_in)){
1879 double scalarl = 0.0;
1880 bool transpl =
false;
1881 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl_in, &scalarl, &transpl);
1882 TEUCHOS_ASSERT(blocked_opl->productRange()->numBlocks() == 1);
1883 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == 1);
1884 opl = Thyra::scale(scalarl,blocked_opl->getBlock(0,0));
1886 if(isPhysicallyBlockedLinearOp(opr_in)){
1887 double scalarr = 0.0;
1888 bool transpr =
false;
1889 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr_in, &scalarr, &transpr);
1890 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == 1);
1891 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == 1);
1892 opr = Thyra::scale(scalarr,blocked_opr->getBlock(0,0));
1895 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1896 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1904 bool transp =
false;
1905 auto crs_op = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalar, &transp);
1906 auto zero_crs = Tpetra::createCrsMatrix<ST,LO,GO,NT>(crs_op->getRowMap());
1907 zero_crs->fillComplete();
1908 opl = Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getDomainMap()),zero_crs);
1911 return opr->clone();
1916 bool transp =
false;
1917 auto crs_op = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalar, &transp);
1918 auto zero_crs = Tpetra::createCrsMatrix<ST,LO,GO,NT>(crs_op->getRowMap());
1919 zero_crs->fillComplete();
1920 opr = Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getDomainMap()),zero_crs);
1923 return opl->clone();
1926 if(isTpetral && isTpetrar){
1930 bool transpl =
false;
1931 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1933 bool transpr =
false;
1934 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1937 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1938 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1941 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1942 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1946#ifdef TEKO_HAVE_EPETRA
1948 const LinearOp implicitOp = Thyra::add(opl,opr);
1951 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1952 Thyra::epetraExtAddTransformer();
1955 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1956 prodTrans->transform(*implicitOp,explicitOp.ptr());
1957 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1958 " + " + opr->getObjectLabel() +
" )");
1962 throw std::logic_error(
"explicitAdd is trying to use Epetra "
1963 "code, but TEKO_HAVE_EPETRA is disabled!");
1980const ModifiableLinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr,
1981 const ModifiableLinearOp & destOp)
1984 if(isPhysicallyBlockedLinearOp(opl)){
1985 TEUCHOS_ASSERT(isPhysicallyBlockedLinearOp(opr));
1987 double scalarl = 0.0;
1988 bool transpl =
false;
1989 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl, &scalarl, &transpl);
1991 double scalarr = 0.0;
1992 bool transpr =
false;
1993 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr, &scalarr, &transpr);
1995 int numRows = blocked_opl->productRange()->numBlocks();
1996 int numCols = blocked_opl->productDomain()->numBlocks();
1997 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numRows);
1998 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == numCols);
2000 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_sum = Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<double>>(destOp);
2001 if(blocked_sum.is_null()) {
2003 blocked_sum = Thyra::defaultBlockedLinearOp<double>();
2005 blocked_sum->beginBlockFill(numRows,numCols);
2006 for(
int r = 0; r < numRows; ++r) {
2007 for(
int c = 0; c < numCols; ++c) {
2008 auto block = explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c)),Teuchos::null);
2009 blocked_sum->setNonconstBlock(r,c,block);
2012 blocked_sum->endBlockFill();
2017 for(
int r = 0; r < numRows; ++r)
2018 for(
int c = 0; c < numCols; ++c)
2019 explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c)),blocked_sum->getNonconstBlock(r,c));
2025 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
2026 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
2028 if(isTpetral && isTpetrar){
2032 bool transpl =
false;
2033 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
2035 bool transpr =
false;
2036 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
2039 RCP<Thyra::LinearOpBase<ST> > explicitOp;
2040 if(destOp!=Teuchos::null)
2041 explicitOp = destOp;
2043 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
2044 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
2047 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
2048 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
2052#ifdef TEKO_HAVE_EPETRA
2054 const LinearOp implicitOp = Thyra::add(opl,opr);
2057 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
2058 Thyra::epetraExtAddTransformer();
2061 RCP<Thyra::LinearOpBase<double> > explicitOp;
2062 if(destOp!=Teuchos::null)
2063 explicitOp = destOp;
2065 explicitOp = prodTrans->createOutputOp();
2068 prodTrans->transform(*implicitOp,explicitOp.ptr());
2069 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
2070 " + " + opr->getObjectLabel() +
" )");
2074 throw std::logic_error(
"explicitAdd is trying to use Epetra "
2075 "code, but TEKO_HAVE_EPETRA is disabled!");
2084const ModifiableLinearOp explicitSum(
const LinearOp & op,
2085 const ModifiableLinearOp & destOp)
2087#ifdef TEKO_HAVE_EPETRA
2089 const RCP<const Epetra_CrsMatrix> epetraOp =
2090 rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*op),
true);
2092 if(destOp==Teuchos::null) {
2093 Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(
new Epetra_CrsMatrix(*epetraOp));
2095 return Thyra::nonconstEpetraLinearOp(epetraDest);
2098 const RCP<Epetra_CrsMatrix> epetraDest =
2099 rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp),
true);
2101 EpetraExt::MatrixMatrix::Add(*epetraOp,
false,1.0,*epetraDest,1.0);
2105 throw std::logic_error(
"explicitSum is trying to use Epetra "
2106 "code, but TEKO_HAVE_EPETRA is disabled!");
2110const LinearOp explicitTranspose(
const LinearOp & op)
2112 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2114 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op,
true);
2115 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
2117 Tpetra::RowMatrixTransposer<ST,LO,GO,NT> transposer(tCrsOp);
2118 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > transOp = transposer.createTranspose();
2120 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getRangeMap()),
2121 Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getDomainMap()),
2125#ifdef TEKO_HAVE_EPETRA
2126 RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2127 TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
2128 "Teko::explicitTranspose Not an Epetra_Operator");
2129 RCP<const Epetra_RowMatrix> eRowMatrixOp
2130 = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
2131 TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
2132 "Teko::explicitTranspose Not an Epetra_RowMatrix");
2135 EpetraExt::RowMatrix_Transpose tranposeOp;
2136 Epetra_RowMatrix & eMat = tranposeOp(
const_cast<Epetra_RowMatrix &
>(*eRowMatrixOp));
2140 Teuchos::RCP<Epetra_CrsMatrix> crsMat
2141 = Teuchos::rcp(
new Epetra_CrsMatrix(
dynamic_cast<Epetra_CrsMatrix &
>(eMat)));
2143 return Thyra::epetraLinearOp(crsMat);
2145 throw std::logic_error(
"explicitTranspose is trying to use Epetra "
2146 "code, but TEKO_HAVE_EPETRA is disabled!");
2151double frobeniusNorm(
const LinearOp & op_in)
2154 double scalar = 1.0;
2157 if(isPhysicallyBlockedLinearOp(op_in)){
2158 bool transp =
false;
2159 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = getPhysicallyBlockedLinearOp(op_in,&scalar,&transp);
2160 TEUCHOS_ASSERT(blocked_op->productRange()->numBlocks() == 1);
2161 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == 1);
2162 op = blocked_op->getBlock(0,0);
2166 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2167 const RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
2168 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
2169 return crsOp->getFrobeniusNorm();
2171#ifdef TEKO_HAVE_EPETRA
2172 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2173 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,
true);
2174 return crsOp->NormFrobenius();
2176 throw std::logic_error(
"frobeniusNorm is trying to use Epetra "
2177 "code, but TEKO_HAVE_EPETRA is disabled!");
2182double oneNorm(
const LinearOp & op)
2184 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2185 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"One norm not currently implemented for Tpetra matrices");
2188#ifdef TEKO_HAVE_EPETRA
2189 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2190 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,
true);
2191 return crsOp->NormOne();
2193 throw std::logic_error(
"oneNorm is trying to use Epetra "
2194 "code, but TEKO_HAVE_EPETRA is disabled!");
2199double infNorm(
const LinearOp & op)
2201 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2203 bool transp =
false;
2204 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
2207 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
2208 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
2211 diag.putScalar(0.0);
2212 for(LO i=0;i<(LO) tCrsOp->getLocalNumRows();i++) {
2213 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
2214 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::local_inds_host_view_type indices;
2215 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::values_host_view_type values;
2216 tCrsOp->getLocalRowView(i,indices,values);
2219 for(LO j=0;j<numEntries;j++)
2220 diag.sumIntoLocalValue(i,std::abs(values(j)));
2222 return diag.normInf()*scalar;
2225#ifdef TEKO_HAVE_EPETRA
2226 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2227 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,
true);
2228 return crsOp->NormInf();
2230 throw std::logic_error(
"infNorm is trying to use Epetra "
2231 "code, but TEKO_HAVE_EPETRA is disabled!");
2236const LinearOp buildDiagonal(
const MultiVector & src,
const std::string & lbl)
2238 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
2239 Thyra::copy(*src->col(0),dst.ptr());
2241 return Thyra::diagonal<double>(dst,lbl);
2244const LinearOp buildInvDiagonal(
const MultiVector & src,
const std::string & lbl)
2246 const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
2247 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
2248 Thyra::reciprocal<double>(*srcV,dst.ptr());
2250 return Thyra::diagonal<double>(dst,lbl);
2254BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> & mvv)
2256 Teuchos::Array<MultiVector> mvA;
2257 Teuchos::Array<VectorSpace> vsA;
2260 std::vector<MultiVector>::const_iterator itr;
2261 for(itr=mvv.begin();itr!=mvv.end();++itr) {
2262 mvA.push_back(*itr);
2263 vsA.push_back((*itr)->range());
2267 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
2268 = Thyra::productVectorSpace<double>(vsA);
2270 return Thyra::defaultProductMultiVector<double>(vs,mvA);
2283Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
2284 const std::vector<int> & indices,
2285 const VectorSpace & vs,
2286 double onValue,
double offValue)
2292 RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
2293 Thyra::put_scalar<double>(offValue,v.ptr());
2296 for(std::size_t i=0;i<indices.size();i++)
2297 Thyra::set_ele<double>(indices[i],onValue,v.ptr());
2326double computeSpectralRad(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
2327 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
2329 typedef Thyra::LinearOpBase<double> OP;
2330 typedef Thyra::MultiVectorBase<double> MV;
2332 int startVectors = 1;
2335 const RCP<MV> ivec = Thyra::createMember(A->domain());
2336 Thyra::randomize(-1.0,1.0,ivec.ptr());
2338 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
2339 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
2341 eigProb->setHermitian(isHermitian);
2344 if(not eigProb->setProblem()) {
2346 return Teuchos::ScalarTraits<double>::nan();
2350 std::string which(
"LM");
2354 Teuchos::ParameterList MyPL;
2355 MyPL.set(
"Verbosity", verbosity );
2356 MyPL.set(
"Which", which );
2357 MyPL.set(
"Block Size", startVectors );
2358 MyPL.set(
"Num Blocks", numBlocks );
2359 MyPL.set(
"Maximum Restarts", restart );
2360 MyPL.set(
"Convergence Tolerance", tol );
2368 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
2371 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
2373 if(solverreturn==Anasazi::Unconverged) {
2374 double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
2375 double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
2377 return -std::abs(std::sqrt(real*real+comp*comp));
2383 double real = eigProb->getSolution().Evals.begin()->realpart;
2384 double comp = eigProb->getSolution().Evals.begin()->imagpart;
2386 return std::abs(std::sqrt(real*real+comp*comp));
2416double computeSmallestMagEig(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
2417 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
2419 typedef Thyra::LinearOpBase<double> OP;
2420 typedef Thyra::MultiVectorBase<double> MV;
2422 int startVectors = 1;
2425 const RCP<MV> ivec = Thyra::createMember(A->domain());
2426 Thyra::randomize(-1.0,1.0,ivec.ptr());
2428 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
2429 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
2431 eigProb->setHermitian(isHermitian);
2434 if(not eigProb->setProblem()) {
2436 return Teuchos::ScalarTraits<double>::nan();
2440 std::string which(
"SM");
2443 Teuchos::ParameterList MyPL;
2444 MyPL.set(
"Verbosity", verbosity );
2445 MyPL.set(
"Which", which );
2446 MyPL.set(
"Block Size", startVectors );
2447 MyPL.set(
"Num Blocks", numBlocks );
2448 MyPL.set(
"Maximum Restarts", restart );
2449 MyPL.set(
"Convergence Tolerance", tol );
2457 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
2460 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
2462 if(solverreturn==Anasazi::Unconverged) {
2464 return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
2468 return std::abs(eigProb->getSolution().Evals.begin()->realpart);
2480ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt)
2484 return getDiagonalOp(A);
2486 return getLumpedMatrix(A);
2488 return getAbsRowSumMatrix(A);
2491 TEUCHOS_TEST_FOR_EXCEPT(
true);
2494 return Teuchos::null;
2505ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp & A,
const Teko::DiagonalType & dt)
2509 return getInvDiagonalOp(A);
2511 return getInvLumpedMatrix(A);
2513 return getAbsRowSumInvMatrix(A);
2516 TEUCHOS_TEST_FOR_EXCEPT(
true);
2519 return Teuchos::null;
2528std::string getDiagonalName(
const DiagonalType & dt)
2556 if(name==
"Diagonal")
2560 if(name==
"AbsRowSum")
2568#ifdef TEKO_HAVE_EPETRA
2569LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp & Op){
2570#ifdef Teko_ENABLE_Isorropia
2571 Teuchos::ParameterList probeList;
2572 Prober prober(G,probeList,
true);
2573 Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(
new Epetra_CrsMatrix(Copy,*G));
2575 prober.probe(Mwrap,*Mat);
2576 return Thyra::epetraLinearOp(Mat);
2580 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
"Probe requires Isorropia");
2585double norm_1(
const MultiVector & v,std::size_t col)
2587 Teuchos::Array<double> n(v->domain()->dim());
2588 Thyra::norms_1<double>(*v,n);
2593double norm_2(
const MultiVector & v,std::size_t col)
2595 Teuchos::Array<double> n(v->domain()->dim());
2596 Thyra::norms_2<double>(*v,n);
2601#ifdef TEKO_HAVE_EPETRA
2602void putScalar(
const ModifiableLinearOp & op,
double scalar)
2606 RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2609 RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
2611 eCrsOp->PutScalar(scalar);
2613 catch (std::exception & e) {
2614 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2616 *out <<
"Teko: putScalar requires an Epetra_CrsMatrix\n";
2617 *out <<
" Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2619 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2621 *out <<
"*** THROWN EXCEPTION ***\n";
2622 *out << e.what() << std::endl;
2623 *out <<
"************************\n";
2630void clipLower(MultiVector & v,
double lowerBound)
2633 using Teuchos::rcp_dynamic_cast;
2639 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2640 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2641 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2643 Teuchos::ArrayRCP<double> values;
2645 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2646 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2647 if(values[j]<lowerBound)
2648 values[j] = lowerBound;
2653void clipUpper(MultiVector & v,
double upperBound)
2656 using Teuchos::rcp_dynamic_cast;
2661 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2662 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2663 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2665 Teuchos::ArrayRCP<double> values;
2667 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2668 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2669 if(values[j]>upperBound)
2670 values[j] = upperBound;
2675void replaceValue(MultiVector & v,
double currentValue,
double newValue)
2678 using Teuchos::rcp_dynamic_cast;
2683 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2684 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2685 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2687 Teuchos::ArrayRCP<double> values;
2689 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2690 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2691 if(values[j]==currentValue)
2692 values[j] = newValue;
2697void columnAverages(
const MultiVector & v,std::vector<double> & averages)
2699 averages.resize(v->domain()->dim());
2702 Thyra::sums<double>(*v,averages);
2705 Thyra::Ordinal rows = v->range()->dim();
2706 for(std::size_t i=0;i<averages.size();i++)
2707 averages[i] = averages[i]/rows;
2710double average(
const MultiVector & v)
2712 Thyra::Ordinal rows = v->range()->dim();
2713 Thyra::Ordinal cols = v->domain()->dim();
2715 std::vector<double> averages;
2716 columnAverages(v,averages);
2719 for(std::size_t i=0;i<averages.size();i++)
2720 sum += averages[i] * rows;
2722 return sum/(rows*cols);
2725bool isPhysicallyBlockedLinearOp(
const LinearOp & op)
2728 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
2729 if (!pblo.is_null())
2734 Thyra::EOpTransp transp = Thyra::NOTRANS;
2735 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
2736 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
2737 pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(wrapped_op);
2738 if (!pblo.is_null())
2744RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(
const LinearOp & op, ST *scalar,
bool *transp)
2747 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
2748 if(!pblo.is_null()){
2755 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
2756 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
2757 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
2758 pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(wrapped_op,
true);
2759 if(!pblo.is_null()){
2761 if(eTransp == Thyra::NOTRANS)
2766 return Teuchos::null;
DiagonalType
Type describing the type of diagonal to construct.
@ BlkDiag
Specifies that a block diagonal approximation is to be used.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ AbsRowSum
Specifies that the diagonal entry is .
@ Diagonal
Specifies that just the diagonal is used.
@ Lumped
Specifies that row sum is used to form a diagonal.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...