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 "Epetra/Teko_EpetraHelpers.hpp" 98 #include "Epetra/Teko_EpetraOperatorWrapper.hpp" 99 #include "Tpetra/Teko_TpetraHelpers.hpp" 100 #include "Tpetra/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" 114 using Teuchos::rcp_dynamic_cast;
116 #ifdef Teko_ENABLE_Isorropia 117 using Isorropia::Epetra::Prober;
120 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
122 Teuchos::RCP<Teuchos::FancyOStream> os =
123 Teuchos::VerboseObjectBase::getDefaultOStream();
131 inline double dist(
int dim,
double * coords,
int row,
int col)
134 for(
int i=0;i<dim;i++)
135 value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
138 return std::sqrt(value);
142 inline double dist(
double * x,
double * y,
double * z,
int stride,
int row,
int col)
145 if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
146 if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
147 if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
150 return std::sqrt(value);
171 RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double * coords,
const Epetra_CrsMatrix & stencil)
174 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
175 stencil.MaxNumEntries()+1),
true);
178 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
179 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
182 for(
int j=0;j<gl->NumMyRows();j++) {
183 int row = gl->GRID(j);
184 double diagValue = 0.0;
189 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
192 for(
int i=0;i<rowSz;i++) {
196 double value = rowData[i];
197 if(value==0)
continue;
201 double d = dist(dim,coords,row,col);
203 diagValue += rowData[i];
211 rowData[rowSz] = -diagValue;
216 rowData[diagInd] = -diagValue;
217 rowInd[diagInd] = row;
221 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
229 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(
int dim,ST * coords,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
232 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
233 stencil.getGlobalMaxNumRowEntries()+1));
236 std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
237 std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
240 for(LO j=0;j< (LO) gl->getNodeNumRows();j++) {
241 GO row = gl->getRowMap()->getGlobalElement(j);
247 stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
250 for(
size_t i=0;i<rowSz;i++) {
254 ST value = rowData[i];
255 if(value==0)
continue;
259 ST d = dist(dim,coords,row,col);
261 diagValue += rowData[i];
269 rowData[rowSz] = -diagValue;
274 rowData[diagInd] = -diagValue;
275 rowInd[diagInd] = row;
279 gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
309 RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double * x,
double * y,
double * z,
int stride,
const Epetra_CrsMatrix & stencil)
312 RCP<Epetra_CrsMatrix> gl = rcp(
new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
313 stencil.MaxNumEntries()+1),
true);
316 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
317 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
320 for(
int j=0;j<gl->NumMyRows();j++) {
321 int row = gl->GRID(j);
322 double diagValue = 0.0;
327 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
330 for(
int i=0;i<rowSz;i++) {
334 double value = rowData[i];
335 if(value==0)
continue;
339 double d = dist(x,y,z,stride,row,col);
341 diagValue += rowData[i];
349 rowData[rowSz] = -diagValue;
354 rowData[diagInd] = -diagValue;
355 rowInd[diagInd] = row;
359 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
367 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
370 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(
new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
371 stencil.getGlobalMaxNumRowEntries()+1),
true);
374 std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
375 std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
378 for(LO j=0;j<(LO) gl->getNodeNumRows();j++) {
379 GO row = gl->getRowMap()->getGlobalElement(j);
385 stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
388 for(
size_t i=0;i<rowSz;i++) {
392 ST value = rowData[i];
393 if(value==0)
continue;
397 ST d = dist(x,y,z,stride,row,col);
399 diagValue += rowData[i];
407 rowData[rowSz] = -diagValue;
412 rowData[diagInd] = -diagValue;
413 rowInd[diagInd] = row;
417 gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
440 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);
448 void update(
double alpha,
const MultiVector & x,
double beta,MultiVector & y)
450 Teuchos::Array<double>
scale;
451 Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
454 scale.push_back(alpha);
455 vec.push_back(x.ptr());
458 Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
462 BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
468 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
469 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
475 upper->beginBlockFill(rows,rows);
477 for(
int i=0;i<rows;i++) {
481 RCP<const Thyra::LinearOpBase<double> > zed
482 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
483 upper->setBlock(i,i,zed);
485 for(
int j=i+1;j<rows;j++) {
487 LinearOp uij = blo->getBlock(i,j);
490 if(uij!=Teuchos::null)
491 upper->setBlock(i,j,uij);
495 upper->endBlockFill();
501 BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill)
507 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
508 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
514 lower->beginBlockFill(rows,rows);
516 for(
int i=0;i<rows;i++) {
520 RCP<const Thyra::LinearOpBase<double> > zed
521 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
522 lower->setBlock(i,i,zed);
524 for(
int j=0;j<i;j++) {
526 LinearOp uij = blo->getBlock(i,j);
529 if(uij!=Teuchos::null)
530 lower->setBlock(i,j,uij);
534 lower->endBlockFill();
558 BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp & blo)
564 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
565 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
571 zeroOp->beginBlockFill(rows,rows);
573 for(
int i=0;i<rows;i++) {
577 RCP<const Thyra::LinearOpBase<double> > zed
578 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
579 zeroOp->setBlock(i,i,zed);
586 bool isZeroOp(
const LinearOp op)
589 if(op==Teuchos::null)
return true;
592 const LinearOp test = rcp_dynamic_cast<
const Thyra::ZeroLinearOpBase<double> >(op);
595 return test!=Teuchos::null;
606 ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp & op)
608 bool isTpetra =
false;
609 RCP<const Epetra_CrsMatrix> eCrsOp;
610 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
614 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
615 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
618 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
620 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
622 else if (!tOp.is_null()){
623 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
627 throw std::logic_error(
"Neither Epetra nor Tpetra");
629 catch (std::exception & e) {
630 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
632 *out <<
"Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
633 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
635 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
637 *out <<
"*** THROWN EXCEPTION ***\n";
638 *out << e.what() << std::endl;
639 *out <<
"************************\n";
646 const RCP<Epetra_Vector> ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
647 Epetra_Vector & diag = *ptrDiag;
651 for(
int i=0;i<eCrsOp->NumMyRows();i++) {
654 eCrsOp->ExtractMyRowView(i,numEntries,values);
657 for(
int j=0;j<numEntries;j++)
658 diag[i] += std::abs(values[j]);
662 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),
"absRowSum( " + op->getObjectLabel() +
" )");
666 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
667 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
671 for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
672 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
673 std::vector<LO> indices(numEntries);
674 std::vector<ST> values(numEntries);
675 Teuchos::ArrayView<const LO> indices_av(indices);
676 Teuchos::ArrayView<const ST> values_av(values);
677 tCrsOp->getLocalRowView(i,indices_av,values_av);
680 for(LO j=0;j<numEntries;j++)
681 diag.sumIntoLocalValue(i,std::abs(values_av[j]));
685 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),
"absRowSum( " + op->getObjectLabel() +
" ))");
699 ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp & op)
701 bool isTpetra =
false;
702 RCP<const Epetra_CrsMatrix> eCrsOp;
703 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
707 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
708 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
711 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
713 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
715 else if (!tOp.is_null()){
716 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
720 throw std::logic_error(
"Neither Epetra nor Tpetra");
722 catch (std::exception & e) {
723 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
725 *out <<
"Teko: getAbsRowSumInvMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
726 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
728 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
730 *out <<
"*** THROWN EXCEPTION ***\n";
731 *out << e.what() << std::endl;
732 *out <<
"************************\n";
739 const RCP<Epetra_Vector> ptrDiag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
740 Epetra_Vector & diag = *ptrDiag;
744 for(
int i=0;i<eCrsOp->NumMyRows();i++) {
747 eCrsOp->ExtractMyRowView(i,numEntries,values);
750 for(
int j=0;j<numEntries;j++)
751 diag[i] += std::abs(values[j]);
753 diag.Reciprocal(diag);
756 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),
"absRowSum( " + op->getObjectLabel() +
" )");
760 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
761 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
765 for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
766 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
767 std::vector<LO> indices(numEntries);
768 std::vector<ST> values(numEntries);
769 Teuchos::ArrayView<const LO> indices_av(indices);
770 Teuchos::ArrayView<const ST> values_av(values);
771 tCrsOp->getLocalRowView(i,indices_av,values_av);
774 for(LO j=0;j<numEntries;j++)
775 diag.sumIntoLocalValue(i,std::abs(values_av[j]));
777 diag.reciprocal(diag);
780 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),
"absRowSum( " + op->getObjectLabel() +
" ))");
793 ModifiableLinearOp getLumpedMatrix(
const LinearOp & op)
795 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
796 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
799 Thyra::assign(ones.ptr(),1.0);
803 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
805 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
816 ModifiableLinearOp getInvLumpedMatrix(
const LinearOp & op)
818 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
819 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
822 Thyra::assign(ones.ptr(),1.0);
825 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
826 Thyra::reciprocal(*diag,diag.ptr());
828 return rcp(
new Thyra::DefaultDiagonalLinearOp<ST>(diag));
842 const ModifiableLinearOp getDiagonalOp(
const LinearOp & op)
844 bool isTpetra =
false;
845 RCP<const Epetra_CrsMatrix> eCrsOp;
846 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
850 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
851 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
854 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
856 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
858 else if (!tOp.is_null()){
859 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
863 throw std::logic_error(
"Neither Epetra nor Tpetra");
865 catch (std::exception & e) {
866 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
868 *out <<
"Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
869 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
871 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
873 *out <<
"*** THROWN EXCEPTION ***\n";
874 *out << e.what() << std::endl;
875 *out <<
"************************\n";
882 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
883 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
886 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
890 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
891 tCrsOp->getLocalDiagCopy(*diag);
894 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
899 const MultiVector getDiagonal(
const LinearOp & op)
901 bool isTpetra =
false;
902 RCP<const Epetra_CrsMatrix> eCrsOp;
903 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
907 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
908 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
911 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
913 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
915 else if (!tOp.is_null()){
916 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
920 throw std::logic_error(
"Neither Epetra nor Tpetra");
922 catch (std::exception & e) {
923 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
925 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
926 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
930 *out <<
"Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
931 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
933 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
935 *out <<
"*** THROWN EXCEPTION ***\n";
936 *out << e.what() << std::endl;
937 *out <<
"************************\n";
944 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
945 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
947 return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
951 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
952 tCrsOp->getLocalDiagCopy(*diag);
955 return Thyra::createVector<ST,LO,GO,NT>(diag,Thyra::createVectorSpace<ST,LO,GO,NT>(tCrsOp->getRowMap()));
960 const MultiVector getDiagonal(
const Teko::LinearOp & A,
const DiagonalType & dt)
962 LinearOp diagOp = Teko::getDiagonalOp(A,dt);
964 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
965 Teuchos::rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
966 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
980 const ModifiableLinearOp getInvDiagonalOp(
const LinearOp & op)
982 bool isTpetra =
false;
983 RCP<const Epetra_CrsMatrix> eCrsOp;
984 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
988 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
989 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
992 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
994 eCrsOp = rcp_dynamic_cast<
const Epetra_CrsMatrix>(eOp->epetra_op(),
true);
996 else if (!tOp.is_null()){
997 tCrsOp = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),
true);
1001 throw std::logic_error(
"Neither Epetra nor Tpetra");
1003 catch (std::exception & e) {
1004 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1006 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<
const Thyra::EpetraLinearOp >(op);
1007 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
1011 *out <<
"Teko: getInvDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
1012 *out <<
" Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
1014 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
1016 *out <<
"*** THROWN EXCEPTION ***\n";
1017 *out << e.what() << std::endl;
1018 *out <<
"************************\n";
1025 const RCP<Epetra_Vector> diag = rcp(
new Epetra_Vector(eCrsOp->RowMap()));
1026 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
1027 diag->Reciprocal(*diag);
1030 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
1034 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
1035 tCrsOp->getLocalDiagCopy(*diag);
1036 diag->reciprocal(*diag);
1039 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),
"inv(diag( " + op->getObjectLabel() +
" ))");
1056 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr)
1058 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1059 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1060 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1062 if(isTpetral && isTpetram && isTpetrar){
1066 bool transpl =
false;
1067 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1069 bool transpm =
false;
1070 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1072 bool transpr =
false;
1073 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1076 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1077 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1080 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1081 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1082 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1083 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1084 explicitCrsOp->resumeFill();
1085 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1086 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1087 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1090 }
else if (isTpetral && !isTpetram && isTpetrar){
1094 bool transpl =
false;
1095 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1097 bool transpr =
false;
1098 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1101 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opm,
true);
1102 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1103 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1104 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()));
1107 tCrsOplm->rightScale(*diagPtr);
1110 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1111 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1114 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1115 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1116 explicitCrsOp->resumeFill();
1117 explicitCrsOp->scale(scalarl*scalarr);
1118 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1119 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1125 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1128 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1129 Thyra::epetraExtDiagScaledMatProdTransformer();
1132 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1133 prodTrans->transform(*implicitOp,explicitOp.ptr());
1134 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1135 " * " + opm->getObjectLabel() +
1136 " * " + opr->getObjectLabel() +
" )");
1157 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr,
1158 const ModifiableLinearOp & destOp)
1160 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1161 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1162 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1164 if(isTpetral && isTpetram && isTpetrar){
1168 bool transpl =
false;
1169 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1171 bool transpm =
false;
1172 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1174 bool transpr =
false;
1175 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1178 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1179 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1182 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1183 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1184 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1185 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1186 explicitCrsOp->resumeFill();
1187 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1188 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1189 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1192 }
else if (isTpetral && !isTpetram && isTpetrar){
1196 bool transpl =
false;
1197 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1199 bool transpr =
false;
1200 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1203 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opm,
true);
1204 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),
true);
1205 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1206 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()));
1209 tCrsOplm->rightScale(*diagPtr);
1212 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1213 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1216 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1217 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,
false,*tCrsOpr,transpr,*explicitCrsOp);
1218 explicitCrsOp->resumeFill();
1219 explicitCrsOp->scale(scalarl*scalarr);
1220 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1221 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1227 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1230 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1231 Thyra::epetraExtDiagScaledMatProdTransformer();
1234 ModifiableLinearOp explicitOp;
1237 if(destOp==Teuchos::null)
1238 explicitOp = prodTrans->createOutputOp();
1240 explicitOp = destOp;
1243 prodTrans->transform(*implicitOp,explicitOp.ptr());
1246 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1247 " * " + opm->getObjectLabel() +
1248 " * " + opr->getObjectLabel() +
" )");
1265 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr)
1267 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1268 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1270 if(isTpetral && isTpetrar){
1273 bool transpl =
false;
1274 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1276 bool transpr =
false;
1277 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1280 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1281 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1284 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1285 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1286 explicitCrsOp->resumeFill();
1287 explicitCrsOp->scale(scalarl*scalarr);
1288 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1289 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1292 }
else if (isTpetral && !isTpetrar){
1296 bool transpl =
false;
1297 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1300 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opr);
1301 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1302 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1303 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()));
1305 explicitCrsOp->rightScale(*diagPtr);
1306 explicitCrsOp->resumeFill();
1307 explicitCrsOp->scale(scalarl);
1308 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1310 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1312 }
else if (!isTpetral && isTpetrar){
1316 bool transpr =
false;
1317 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1320 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opl);
1321 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1322 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1323 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()));
1325 explicitCrsOp->leftScale(*diagPtr);
1326 explicitCrsOp->resumeFill();
1327 explicitCrsOp->scale(scalarr);
1328 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1330 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1335 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1338 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1339 = Thyra::epetraExtDiagScalingTransformer();
1343 if(not prodTrans->isCompatible(*implicitOp))
1344 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1347 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1348 prodTrans->transform(*implicitOp,explicitOp.ptr());
1349 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1350 " * " + opr->getObjectLabel() +
" )");
1369 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr,
1370 const ModifiableLinearOp & destOp)
1372 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1373 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1375 if(isTpetral && isTpetrar){
1379 bool transpl =
false;
1380 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1382 bool transpr =
false;
1383 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1386 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1387 if(destOp!=Teuchos::null)
1388 explicitOp = destOp;
1390 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1391 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1394 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1395 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1396 explicitCrsOp->resumeFill();
1397 explicitCrsOp->scale(scalarl*scalarr);
1398 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1399 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1402 }
else if (isTpetral && !isTpetrar){
1406 bool transpl =
false;
1407 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1410 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opr);
1411 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),
true);
1412 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1415 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()));
1416 explicitCrsOp->rightScale(*diagPtr);
1417 explicitCrsOp->resumeFill();
1418 explicitCrsOp->scale(scalarl);
1419 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1420 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1422 }
else if (!isTpetral && isTpetrar){
1426 bool transpr =
false;
1427 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1430 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<ST> >(opl,
true);
1431 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<
const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),
true);
1432 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<
const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),
true);
1435 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()));
1436 explicitCrsOp->leftScale(*diagPtr);
1437 explicitCrsOp->resumeFill();
1438 explicitCrsOp->scale(scalarr);
1439 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1440 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1445 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1449 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1450 = Thyra::epetraExtDiagScalingTransformer();
1454 if(not prodTrans->isCompatible(*implicitOp))
1455 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1458 ModifiableLinearOp explicitOp;
1461 if(destOp==Teuchos::null)
1462 explicitOp = prodTrans->createOutputOp();
1464 explicitOp = destOp;
1467 prodTrans->transform(*implicitOp,explicitOp.ptr());
1470 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1471 " * " + opr->getObjectLabel() +
" )");
1487 const LinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr)
1489 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1490 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1492 if(isTpetral && isTpetrar){
1496 bool transpl =
false;
1497 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1499 bool transpr =
false;
1500 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1503 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1504 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1507 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1508 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1513 const LinearOp implicitOp = Thyra::add(opl,opr);
1516 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1517 Thyra::epetraExtAddTransformer();
1520 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1521 prodTrans->transform(*implicitOp,explicitOp.ptr());
1522 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1523 " + " + opr->getObjectLabel() +
" )");
1541 const ModifiableLinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr,
1542 const ModifiableLinearOp & destOp)
1544 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1545 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1547 if(isTpetral && isTpetrar){
1551 bool transpl =
false;
1552 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1554 bool transpr =
false;
1555 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1558 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1559 if(destOp!=Teuchos::null)
1560 explicitOp = destOp;
1562 explicitOp = rcp(
new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1563 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1566 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1567 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1573 const LinearOp implicitOp = Thyra::add(opl,opr);
1576 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1577 Thyra::epetraExtAddTransformer();
1580 RCP<Thyra::LinearOpBase<double> > explicitOp;
1581 if(destOp!=Teuchos::null)
1582 explicitOp = destOp;
1584 explicitOp = prodTrans->createOutputOp();
1587 prodTrans->transform(*implicitOp,explicitOp.ptr());
1588 explicitOp->setObjectLabel(
"explicit( " + opl->getObjectLabel() +
1589 " + " + opr->getObjectLabel() +
" )");
1599 const ModifiableLinearOp explicitSum(
const LinearOp & op,
1600 const ModifiableLinearOp & destOp)
1603 const RCP<const Epetra_CrsMatrix> epetraOp =
1604 rcp_dynamic_cast<
const Epetra_CrsMatrix>(get_Epetra_Operator(*op),
true);
1606 if(destOp==Teuchos::null) {
1607 Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(
new Epetra_CrsMatrix(*epetraOp));
1609 return Thyra::nonconstEpetraLinearOp(epetraDest);
1612 const RCP<Epetra_CrsMatrix> epetraDest =
1613 rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp),
true);
1615 EpetraExt::MatrixMatrix::Add(*epetraOp,
false,1.0,*epetraDest,1.0);
1620 const LinearOp explicitTranspose(
const LinearOp & op)
1622 RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1623 TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
1624 "Teko::explicitTranspose Not an Epetra_Operator");
1625 RCP<const Epetra_RowMatrix> eRowMatrixOp
1626 = Teuchos::rcp_dynamic_cast<
const Epetra_RowMatrix>(eOp);
1627 TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
1628 "Teko::explicitTranspose Not an Epetra_RowMatrix");
1631 EpetraExt::RowMatrix_Transpose tranposeOp;
1632 Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
1636 Teuchos::RCP<Epetra_CrsMatrix> crsMat
1637 = Teuchos::rcp(
new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
1639 return Thyra::epetraLinearOp(crsMat);
1642 const LinearOp buildDiagonal(
const MultiVector & src,
const std::string & lbl)
1644 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
1645 Thyra::copy(*src->col(0),dst.ptr());
1647 return Thyra::diagonal<double>(dst,lbl);
1650 const LinearOp buildInvDiagonal(
const MultiVector & src,
const std::string & lbl)
1652 const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
1653 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
1654 Thyra::reciprocal<double>(*srcV,dst.ptr());
1656 return Thyra::diagonal<double>(dst,lbl);
1660 BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> & mvv)
1662 Teuchos::Array<MultiVector> mvA;
1663 Teuchos::Array<VectorSpace> vsA;
1666 std::vector<MultiVector>::const_iterator itr;
1667 for(itr=mvv.begin();itr!=mvv.end();++itr) {
1668 mvA.push_back(*itr);
1669 vsA.push_back((*itr)->range());
1673 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
1674 = Thyra::productVectorSpace<double>(vsA);
1676 return Thyra::defaultProductMultiVector<double>(vs,mvA);
1689 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
1690 const std::vector<int> & indices,
1691 const VectorSpace & vs,
1692 double onValue,
double offValue)
1698 RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
1699 Thyra::put_scalar<double>(offValue,v.ptr());
1702 for(std::size_t i=0;i<indices.size();i++)
1703 Thyra::set_ele<double>(indices[i],onValue,v.ptr());
1732 double computeSpectralRad(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
1733 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
1735 typedef Thyra::LinearOpBase<double> OP;
1736 typedef Thyra::MultiVectorBase<double> MV;
1738 int startVectors = 1;
1741 const RCP<MV> ivec = Thyra::createMember(A->domain());
1742 Thyra::randomize(-1.0,1.0,ivec.ptr());
1744 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1745 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1747 eigProb->setHermitian(isHermitian);
1750 if(not eigProb->setProblem()) {
1752 return Teuchos::ScalarTraits<double>::nan();
1756 std::string which(
"LM");
1760 Teuchos::ParameterList MyPL;
1761 MyPL.set(
"Verbosity", verbosity );
1762 MyPL.set(
"Which", which );
1763 MyPL.set(
"Block Size", startVectors );
1764 MyPL.set(
"Num Blocks", numBlocks );
1765 MyPL.set(
"Maximum Restarts", restart );
1766 MyPL.set(
"Convergence Tolerance", tol );
1774 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1777 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1779 if(solverreturn==Anasazi::Unconverged) {
1780 double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
1781 double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
1783 return -std::abs(std::sqrt(real*real+comp*comp));
1789 double real = eigProb->getSolution().Evals.begin()->realpart;
1790 double comp = eigProb->getSolution().Evals.begin()->imagpart;
1792 return std::abs(std::sqrt(real*real+comp*comp));
1822 double computeSmallestMagEig(
const RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
1823 bool isHermitian,
int numBlocks,
int restart,
int verbosity)
1825 typedef Thyra::LinearOpBase<double> OP;
1826 typedef Thyra::MultiVectorBase<double> MV;
1828 int startVectors = 1;
1831 const RCP<MV> ivec = Thyra::createMember(A->domain());
1832 Thyra::randomize(-1.0,1.0,ivec.ptr());
1834 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1835 = rcp(
new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1837 eigProb->setHermitian(isHermitian);
1840 if(not eigProb->setProblem()) {
1842 return Teuchos::ScalarTraits<double>::nan();
1846 std::string which(
"SM");
1849 Teuchos::ParameterList MyPL;
1850 MyPL.set(
"Verbosity", verbosity );
1851 MyPL.set(
"Which", which );
1852 MyPL.set(
"Block Size", startVectors );
1853 MyPL.set(
"Num Blocks", numBlocks );
1854 MyPL.set(
"Maximum Restarts", restart );
1855 MyPL.set(
"Convergence Tolerance", tol );
1863 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1866 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1868 if(solverreturn==Anasazi::Unconverged) {
1870 return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1874 return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1886 ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt)
1890 return getDiagonalOp(A);
1892 return getLumpedMatrix(A);
1894 return getAbsRowSumMatrix(A);
1897 TEUCHOS_TEST_FOR_EXCEPT(
true);
1900 return Teuchos::null;
1911 ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp & A,
const Teko::DiagonalType & dt)
1915 return getInvDiagonalOp(A);
1917 return getInvLumpedMatrix(A);
1919 return getAbsRowSumInvMatrix(A);
1922 TEUCHOS_TEST_FOR_EXCEPT(
true);
1925 return Teuchos::null;
1934 std::string getDiagonalName(
const DiagonalType & dt)
1962 if(name==
"Diagonal")
1966 if(name==
"AbsRowSum")
1974 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp & Op){
1975 #ifdef Teko_ENABLE_Isorropia 1976 Teuchos::ParameterList probeList;
1977 Prober prober(G,probeList,
true);
1978 Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(
new Epetra_CrsMatrix(Copy,*G));
1980 prober.probe(Mwrap,*Mat);
1981 return Thyra::epetraLinearOp(Mat);
1983 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
"Probe requires Isorropia");
1987 double norm_1(
const MultiVector & v,std::size_t col)
1989 Teuchos::Array<double> n(v->domain()->dim());
1990 Thyra::norms_1<double>(*v,n);
1995 double norm_2(
const MultiVector & v,std::size_t col)
1997 Teuchos::Array<double> n(v->domain()->dim());
1998 Thyra::norms_2<double>(*v,n);
2003 void putScalar(
const ModifiableLinearOp & op,
double scalar)
2007 RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2010 RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
2012 eCrsOp->PutScalar(scalar);
2014 catch (std::exception & e) {
2015 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2017 *out <<
"Teko: putScalar requires an Epetra_CrsMatrix\n";
2018 *out <<
" Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2020 *out <<
" Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2022 *out <<
"*** THROWN EXCEPTION ***\n";
2023 *out << e.what() << std::endl;
2024 *out <<
"************************\n";
2030 void clipLower(MultiVector & v,
double lowerBound)
2033 using Teuchos::rcp_dynamic_cast;
2039 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2040 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2041 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2043 Teuchos::ArrayRCP<double> values;
2045 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2046 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2047 if(values[j]<lowerBound)
2048 values[j] = lowerBound;
2053 void clipUpper(MultiVector & v,
double upperBound)
2056 using Teuchos::rcp_dynamic_cast;
2061 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2062 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2063 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2065 Teuchos::ArrayRCP<double> values;
2067 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2068 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2069 if(values[j]>upperBound)
2070 values[j] = upperBound;
2075 void replaceValue(MultiVector & v,
double currentValue,
double newValue)
2078 using Teuchos::rcp_dynamic_cast;
2083 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2084 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2085 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),
true);
2087 Teuchos::ArrayRCP<double> values;
2089 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2090 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2091 if(values[j]==currentValue)
2092 values[j] = newValue;
2097 void columnAverages(
const MultiVector & v,std::vector<double> & averages)
2099 averages.resize(v->domain()->dim());
2102 Thyra::sums<double>(*v,averages);
2105 Thyra::Ordinal rows = v->range()->dim();
2106 for(std::size_t i=0;i<averages.size();i++)
2107 averages[i] = averages[i]/rows;
2110 double average(
const MultiVector & v)
2112 Thyra::Ordinal rows = v->range()->dim();
2113 Thyra::Ordinal cols = v->domain()->dim();
2115 std::vector<double> averages;
2116 columnAverages(v,averages);
2119 for(std::size_t i=0;i<averages.size();i++)
2120 sum += averages[i] * rows;
2122 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 .