55 #ifndef __Teko_Utilities_hpp__ 56 #define __Teko_Utilities_hpp__ 58 #include "Epetra_CrsMatrix.h" 59 #include "Tpetra_CrsMatrix.hpp" 62 #include "Teuchos_VerboseObject.hpp" 65 #include "Thyra_LinearOpBase.hpp" 66 #include "Thyra_PhysicallyBlockedLinearOpBase.hpp" 67 #include "Thyra_ProductVectorSpaceBase.hpp" 68 #include "Thyra_VectorSpaceBase.hpp" 69 #include "Thyra_ProductMultiVectorBase.hpp" 70 #include "Thyra_MultiVectorStdOps.hpp" 71 #include "Thyra_MultiVectorBase.hpp" 72 #include "Thyra_VectorBase.hpp" 73 #include "Thyra_VectorStdOps.hpp" 74 #include "Thyra_DefaultBlockedLinearOp.hpp" 75 #include "Thyra_DefaultMultipliedLinearOp.hpp" 76 #include "Thyra_DefaultScaledAdjointLinearOp.hpp" 77 #include "Thyra_DefaultAddedLinearOp.hpp" 78 #include "Thyra_DefaultIdentityLinearOp.hpp" 79 #include "Thyra_DefaultZeroLinearOp.hpp" 81 #include "Teko_ConfigDefs.hpp" 84 #ifndef _MSC_EXTENSIONS 85 #define _MSC_EXTENSIONS 86 #define TEKO_DEFINED_MSC_EXTENSIONS 92 #define Teko_DEBUG_INT 5 96 using Thyra::multiply;
99 using Thyra::identity;
101 using Thyra::block2x2;
102 using Thyra::block2x1;
103 using Thyra::block1x2;
123 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double * coords,
const Epetra_CrsMatrix & stencil);
124 Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(
int dim,ST * coords,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
148 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double * x,
double * y,
double * z,
int stride,
const Epetra_CrsMatrix & stencil);
149 Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
157 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
161 #ifndef Teko_DEBUG_OFF 163 #define Teko_DEBUG_EXPR(str) str 164 #define Teko_DEBUG_MSG(str,level) if(level<=Teko_DEBUG_INT) { \ 165 Teuchos::RCP<Teuchos::FancyOStream> out = Teko::getOutputStream(); \ 166 *out << "Teko: " << str << std::endl; } 167 #define Teko_DEBUG_MSG_BEGIN(level) if(level<=Teko_DEBUG_INT) { \ 168 Teko::getOutputStream()->pushTab(3); \ 169 *Teko::getOutputStream() << "Teko: Begin debug MSG\n"; \ 170 std::ostream & DEBUG_STREAM = *Teko::getOutputStream(); \ 171 Teko::getOutputStream()->pushTab(3); 172 #define Teko_DEBUG_MSG_END() Teko::getOutputStream()->popTab(); \ 173 *Teko::getOutputStream() << "Teko: End debug MSG\n"; \ 174 Teko::getOutputStream()->popTab(); } 175 #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3) 176 #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab() 177 #define Teko_DEBUG_SCOPE(str,level) 188 #define Teko_DEBUG_EXPR(str) 189 #define Teko_DEBUG_MSG(str,level) 190 #define Teko_DEBUG_MSG_BEGIN(level) if(false) { \ 191 std::ostream & DEBUG_STREAM = *Teko::getOutputStream(); 192 #define Teko_DEBUG_MSG_END() } 193 #define Teko_DEBUG_PUSHTAB() 194 #define Teko_DEBUG_POPTAB() 195 #define Teko_DEBUG_SCOPE(str,level) 199 typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
206 typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
207 typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
213 inline const MultiVector
toMultiVector(
const BlockedMultiVector & bmv) {
return bmv; }
217 {
return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv); }
221 {
return bmv->productSpace()->numBlocks(); }
224 inline MultiVector
getBlock(
int i,
const BlockedMultiVector & bmv)
225 {
return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i)); }
229 {
return v->clone_mv(); }
232 inline MultiVector
copyAndInit(
const MultiVector & v,
double scalar)
233 { MultiVector mv = v->clone_mv(); Thyra::assign(mv.ptr(),scalar);
return mv; }
236 inline BlockedMultiVector
deepcopy(
const BlockedMultiVector & v)
252 inline MultiVector
datacopy(
const MultiVector & src,MultiVector & dst)
254 if(dst==Teuchos::null)
257 bool rangeCompat = src->range()->isCompatible(*dst->range());
258 bool domainCompat = src->domain()->isCompatible(*dst->domain());
260 if(not (rangeCompat && domainCompat))
264 Thyra::assign<double>(dst.ptr(),*src);
281 inline BlockedMultiVector
datacopy(
const BlockedMultiVector & src,BlockedMultiVector & dst)
283 if(dst==Teuchos::null)
286 bool rangeCompat = src->range()->isCompatible(*dst->range());
287 bool domainCompat = src->domain()->isCompatible(*dst->domain());
289 if(not (rangeCompat && domainCompat))
293 Thyra::assign<double>(dst.ptr(),*src);
298 BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> & mvs);
310 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
311 const std::vector<int> & indices,
312 const VectorSpace & vs,
313 double onValue=1.0,
double offValue=0.0);
321 typedef Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > BlockedLinearOp;
322 typedef Teuchos::RCP<const Thyra::LinearOpBase<ST> > LinearOp;
323 typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > InverseLinearOp;
324 typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > ModifiableLinearOp;
327 inline LinearOp
zero(
const VectorSpace & vs)
328 {
return Thyra::zero<ST>(vs,vs); }
331 void putScalar(
const ModifiableLinearOp & op,
double scalar);
335 {
return lo->range(); }
339 {
return lo->domain(); }
344 Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
345 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
351 Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
352 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
356 inline LinearOp
toLinearOp(BlockedLinearOp & blo) {
return blo; }
359 inline const LinearOp
toLinearOp(
const BlockedLinearOp & blo) {
return blo; }
362 inline LinearOp
toLinearOp(ModifiableLinearOp & blo) {
return blo; }
365 inline const LinearOp
toLinearOp(
const ModifiableLinearOp & blo) {
return blo; }
369 {
return blo->productRange()->numBlocks(); }
373 {
return blo->productDomain()->numBlocks(); }
376 inline LinearOp
getBlock(
int i,
int j,
const BlockedLinearOp & blo)
377 {
return blo->getBlock(i,j); }
380 inline void setBlock(
int i,
int j,BlockedLinearOp & blo,
const LinearOp & lo)
381 {
return blo->setBlock(i,j,lo); }
385 {
return rcp(
new Thyra::DefaultBlockedLinearOp<double>()); }
399 { blo->beginBlockFill(rowCnt,colCnt); }
411 { blo->beginBlockFill(); }
415 { blo->endBlockFill(); }
418 BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill=
true);
421 BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill=
true);
442 BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp & blo);
445 bool isZeroOp(
const LinearOp op);
455 ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp & op);
465 ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp & op);
474 ModifiableLinearOp getLumpedMatrix(
const LinearOp & op);
484 ModifiableLinearOp getInvLumpedMatrix(
const LinearOp & op);
509 void applyOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha=1.0,
double beta=0.0);
529 inline void applyOp(
const LinearOp & A,
const BlockedMultiVector & x,BlockedMultiVector & y,
double alpha=1.0,
double beta=0.0)
531 applyOp(A,x_mv,y_mv,alpha,beta); }
542 void update(
double alpha,
const MultiVector & x,
double beta,MultiVector & y);
545 inline void update(
double alpha,
const BlockedMultiVector & x,
double beta,BlockedMultiVector & y)
547 update(alpha,x_mv,beta,y_mv); }
550 inline void scale(
const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
553 inline void scale(
const double alpha,BlockedMultiVector & x)
557 inline LinearOp
scale(
const double alpha,ModifiableLinearOp & a)
558 {
return Thyra::nonconstScale(alpha,a); }
561 inline LinearOp
adjoint(ModifiableLinearOp & a)
562 {
return Thyra::nonconstAdjoint(a); }
580 const ModifiableLinearOp getDiagonalOp(
const LinearOp & op);
587 const MultiVector getDiagonal(
const LinearOp & op);
600 const ModifiableLinearOp getInvDiagonalOp(
const LinearOp & op);
614 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr);
630 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr,
631 const ModifiableLinearOp & destOp);
643 const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr);
658 const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr,
659 const ModifiableLinearOp & destOp);
671 const LinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr);
685 const ModifiableLinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr,
686 const ModifiableLinearOp & destOp);
690 const ModifiableLinearOp explicitSum(
const LinearOp & opl,
691 const ModifiableLinearOp & destOp);
696 const LinearOp explicitTranspose(
const LinearOp & op);
701 const LinearOp buildDiagonal(
const MultiVector & v,
const std::string & lbl=
"ANYM");
706 const LinearOp buildInvDiagonal(
const MultiVector & v,
const std::string & lbl=
"ANYM");
733 double computeSpectralRad(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
734 bool isHermitian=
false,
int numBlocks=5,
int restart=0,
int verbosity=0);
759 double computeSmallestMagEig(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
760 bool isHermitian=
false,
int numBlocks=5,
int restart=0,
int verbosity=0);
778 ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt);
788 ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt);
795 const MultiVector getDiagonal(
const LinearOp & op,
const DiagonalType & dt);
803 std::string getDiagonalName(
const DiagonalType & dt);
813 DiagonalType getDiagonalType(std::string name);
815 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp & Op);
819 double norm_1(
const MultiVector & v,std::size_t col);
823 double norm_2(
const MultiVector & v,std::size_t col);
828 void clipLower(MultiVector & v,
double lowerBound);
833 void clipUpper(MultiVector & v,
double upperBound);
838 void replaceValue(MultiVector & v,
double currentValue,
double newValue);
842 void columnAverages(
const MultiVector & v,std::vector<double> & averages);
846 double average(
const MultiVector & v);
851 #ifdef TEKO_DEFINED_MSC_EXTENSIONS 852 #undef _MSC_EXTENSIONS LinearOp adjoint(ModifiableLinearOp &a)
Construct an implicit adjoint of the linear operators.
MultiVector toMultiVector(BlockedMultiVector &bmv)
Convert to a MultiVector from a BlockedMultiVector.
int blockCount(const BlockedMultiVector &bmv)
Get the column count in a block linear operator.
Specifies that a block diagonal approximation is to be used.
Specifies that just the diagonal is used.
MultiVector deepcopy(const MultiVector &v)
Perform a deep copy of the vector.
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.
BlockedLinearOp toBlockedLinearOp(LinearOp &clo)
Converts a LinearOp to a BlockedLinearOp.
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
MultiVector copyAndInit(const MultiVector &v, double scalar)
Perform a deep copy of the vector.
VectorSpace domainSpace(const LinearOp &lo)
Get the domain space of a linear operator.
void beginBlockFill(BlockedLinearOp &blo, int rowCnt, int colCnt)
Let the blocked operator know that you are going to set the sub blocks.
const BlockedMultiVector toBlockedMultiVector(const MultiVector &bmv)
Convert to a BlockedMultiVector from a MultiVector.
VectorSpace rangeSpace(const LinearOp &lo)
Get the range space of a linear operator.
LinearOp toLinearOp(BlockedLinearOp &blo)
Convert to a LinearOp from a BlockedLinearOp.
LinearOp zero(const VectorSpace &vs)
Build a square zero operator from a single vector space.
MultiVector datacopy(const MultiVector &src, MultiVector &dst)
Copy the contents of a multivector to a destination vector.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
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.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
Specifies that the diagonal entry is .