55#ifndef __Teko_Utilities_hpp__
56#define __Teko_Utilities_hpp__
58#include "Teko_ConfigDefs.hpp"
60#ifdef TEKO_HAVE_EPETRA
61#include "Epetra_CrsMatrix.h"
64#include "Tpetra_CrsMatrix.hpp"
67#include "Teuchos_VerboseObject.hpp"
70#include "Thyra_LinearOpBase.hpp"
71#include "Thyra_PhysicallyBlockedLinearOpBase.hpp"
72#include "Thyra_ProductVectorSpaceBase.hpp"
73#include "Thyra_VectorSpaceBase.hpp"
74#include "Thyra_ProductMultiVectorBase.hpp"
75#include "Thyra_MultiVectorStdOps.hpp"
76#include "Thyra_MultiVectorBase.hpp"
77#include "Thyra_VectorBase.hpp"
78#include "Thyra_VectorStdOps.hpp"
79#include "Thyra_DefaultBlockedLinearOp.hpp"
80#include "Thyra_DefaultMultipliedLinearOp.hpp"
81#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
82#include "Thyra_DefaultAddedLinearOp.hpp"
83#include "Thyra_DefaultIdentityLinearOp.hpp"
84#include "Thyra_DefaultZeroLinearOp.hpp"
87#ifndef _MSC_EXTENSIONS
88#define _MSC_EXTENSIONS
89#define TEKO_DEFINED_MSC_EXTENSIONS
95#define Teko_DEBUG_INT 5
102using Thyra::identity;
104using Thyra::block2x2;
105using Thyra::block2x1;
106using Thyra::block1x2;
126#ifdef TEKO_HAVE_EPETRA
127Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double * coords,
const Epetra_CrsMatrix & stencil);
130Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(
int dim,ST * coords,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
154#ifdef TEKO_HAVE_EPETRA
155Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double * x,
double * y,
double * z,
int stride,
const Epetra_CrsMatrix & stencil);
158Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,
const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
166const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
170#ifndef Teko_DEBUG_OFF
172 #define Teko_DEBUG_EXPR(str) str
173 #define Teko_DEBUG_MSG(str,level) if(level<=Teko_DEBUG_INT) { \
174 Teuchos::RCP<Teuchos::FancyOStream> out = Teko::getOutputStream(); \
175 *out << "Teko: " << str << std::endl; }
176 #define Teko_DEBUG_MSG_BEGIN(level) if(level<=Teko_DEBUG_INT) { \
177 Teko::getOutputStream()->pushTab(3); \
178 *Teko::getOutputStream() << "Teko: Begin debug MSG\n"; \
179 std::ostream & DEBUG_STREAM = *Teko::getOutputStream(); \
180 Teko::getOutputStream()->pushTab(3);
181 #define Teko_DEBUG_MSG_END() Teko::getOutputStream()->popTab(); \
182 *Teko::getOutputStream() << "Teko: End debug MSG\n"; \
183 Teko::getOutputStream()->popTab(); }
184 #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3)
185 #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab()
186 #define Teko_DEBUG_SCOPE(str,level)
197 #define Teko_DEBUG_EXPR(str)
198 #define Teko_DEBUG_MSG(str,level)
199 #define Teko_DEBUG_MSG_BEGIN(level) if(false) { \
200 std::ostream & DEBUG_STREAM = *Teko::getOutputStream();
201 #define Teko_DEBUG_MSG_END() }
202 #define Teko_DEBUG_PUSHTAB()
203 #define Teko_DEBUG_POPTAB()
204 #define Teko_DEBUG_SCOPE(str,level)
208typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
215typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
216typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
222inline const MultiVector
toMultiVector(
const BlockedMultiVector & bmv) {
return bmv; }
226{
return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv); }
230{
return bmv->productSpace()->numBlocks(); }
233inline MultiVector
getBlock(
int i,
const BlockedMultiVector & bmv)
234{
return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i)); }
238{
return v->clone_mv(); }
241inline MultiVector
copyAndInit(
const MultiVector & v,
double scalar)
242{ MultiVector mv = v->clone_mv(); Thyra::assign(mv.ptr(),scalar);
return mv; }
245inline BlockedMultiVector
deepcopy(
const BlockedMultiVector & v)
261inline MultiVector
datacopy(
const MultiVector & src,MultiVector & dst)
263 if(dst==Teuchos::null)
266 bool rangeCompat = src->range()->isCompatible(*dst->range());
267 bool domainCompat = src->domain()->isCompatible(*dst->domain());
269 if(not (rangeCompat && domainCompat))
273 Thyra::assign<double>(dst.ptr(),*src);
290inline BlockedMultiVector
datacopy(
const BlockedMultiVector & src,BlockedMultiVector & dst)
292 if(dst==Teuchos::null)
295 bool rangeCompat = src->range()->isCompatible(*dst->range());
296 bool domainCompat = src->domain()->isCompatible(*dst->domain());
298 if(not (rangeCompat && domainCompat))
302 Thyra::assign<double>(dst.ptr(),*src);
307BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> & mvs);
319Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
320 const std::vector<int> & indices,
321 const VectorSpace & vs,
322 double onValue=1.0,
double offValue=0.0);
330typedef Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > BlockedLinearOp;
331typedef Teuchos::RCP<const Thyra::LinearOpBase<ST> > LinearOp;
332typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > InverseLinearOp;
333typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > ModifiableLinearOp;
336inline LinearOp
zero(
const VectorSpace & vs)
337{
return Thyra::zero<ST>(vs,vs); }
340#ifdef TEKO_HAVE_EPETRA
341void putScalar(
const ModifiableLinearOp & op,
double scalar);
346{
return lo->range(); }
350{
return lo->domain(); }
355 Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
356 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
362 Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
363 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
367inline LinearOp
toLinearOp(BlockedLinearOp & blo) {
return blo; }
370inline const LinearOp
toLinearOp(
const BlockedLinearOp & blo) {
return blo; }
373inline LinearOp
toLinearOp(ModifiableLinearOp & blo) {
return blo; }
376inline const LinearOp
toLinearOp(
const ModifiableLinearOp & blo) {
return blo; }
380{
return blo->productRange()->numBlocks(); }
384{
return blo->productDomain()->numBlocks(); }
387inline LinearOp
getBlock(
int i,
int j,
const BlockedLinearOp & blo)
388{
return blo->getBlock(i,j); }
391inline void setBlock(
int i,
int j,BlockedLinearOp & blo,
const LinearOp & lo)
392{
return blo->setBlock(i,j,lo); }
396{
return rcp(
new Thyra::DefaultBlockedLinearOp<double>()); }
410{ blo->beginBlockFill(rowCnt,colCnt); }
422{ blo->beginBlockFill(); }
426{ blo->endBlockFill(); }
429BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill=
true);
432BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp & blo,
bool callEndBlockFill=
true);
453BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp & blo);
456bool isZeroOp(
const LinearOp op);
466ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp & op);
476ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp & op);
485ModifiableLinearOp getLumpedMatrix(
const LinearOp & op);
495ModifiableLinearOp getInvLumpedMatrix(
const LinearOp & op);
520void applyOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha=1.0,
double beta=0.0);
541void applyTransposeOp(
const LinearOp & A,
const MultiVector & x,MultiVector & y,
double alpha=1.0,
double beta=0.0);
561inline void applyOp(
const LinearOp & A,
const BlockedMultiVector & x,BlockedMultiVector & y,
double alpha=1.0,
double beta=0.0)
563 applyOp(A,x_mv,y_mv,alpha,beta); }
583inline void applyTransposeOp(
const LinearOp & A,
const BlockedMultiVector & x,BlockedMultiVector & y,
double alpha=1.0,
double beta=0.0)
585 applyTransposeOp(A,x_mv,y_mv,alpha,beta); }
596void update(
double alpha,
const MultiVector & x,
double beta,MultiVector & y);
599inline void update(
double alpha,
const BlockedMultiVector & x,
double beta,BlockedMultiVector & y)
601 update(alpha,x_mv,beta,y_mv); }
604inline void scale(
const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
607inline void scale(
const double alpha,BlockedMultiVector & x)
611inline LinearOp
scale(
const double alpha,ModifiableLinearOp & a)
612{
return Thyra::nonconstScale(alpha,a); }
615inline LinearOp
adjoint(ModifiableLinearOp & a)
616{
return Thyra::nonconstAdjoint(a); }
634const ModifiableLinearOp getDiagonalOp(
const LinearOp & op);
641const MultiVector getDiagonal(
const LinearOp & op);
654const ModifiableLinearOp getInvDiagonalOp(
const LinearOp & op);
668const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr);
684const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opm,
const LinearOp & opr,
685 const ModifiableLinearOp & destOp);
697const LinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr);
712const ModifiableLinearOp explicitMultiply(
const LinearOp & opl,
const LinearOp & opr,
713 const ModifiableLinearOp & destOp);
725const LinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr);
739const ModifiableLinearOp explicitAdd(
const LinearOp & opl,
const LinearOp & opr,
740 const ModifiableLinearOp & destOp);
744const ModifiableLinearOp explicitSum(
const LinearOp & opl,
745 const ModifiableLinearOp & destOp);
750const LinearOp explicitTranspose(
const LinearOp & op);
754double frobeniusNorm(
const LinearOp & op);
755double oneNorm(
const LinearOp & op);
756double infNorm(
const LinearOp & op);
761const LinearOp buildDiagonal(
const MultiVector & v,
const std::string & lbl=
"ANYM");
766const LinearOp buildInvDiagonal(
const MultiVector & v,
const std::string & lbl=
"ANYM");
793double computeSpectralRad(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
794 bool isHermitian=
false,
int numBlocks=5,
int restart=0,
int verbosity=0);
819double computeSmallestMagEig(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & A,
double tol,
820 bool isHermitian=
false,
int numBlocks=5,
int restart=0,
int verbosity=0);
838ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt);
848ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp & A,
const DiagonalType & dt);
855const MultiVector getDiagonal(
const LinearOp & op,
const DiagonalType & dt);
863std::string getDiagonalName(
const DiagonalType & dt);
873DiagonalType getDiagonalType(std::string name);
875#ifdef TEKO_HAVE_EPETRA
876LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp & Op);
881double norm_1(
const MultiVector & v,std::size_t col);
885double norm_2(
const MultiVector & v,std::size_t col);
890void clipLower(MultiVector & v,
double lowerBound);
895void clipUpper(MultiVector & v,
double upperBound);
900void replaceValue(MultiVector & v,
double currentValue,
double newValue);
904void columnAverages(
const MultiVector & v,std::vector<double> & averages);
908double average(
const MultiVector & v);
912bool isPhysicallyBlockedLinearOp(
const LinearOp & op);
916Teuchos::RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(
const LinearOp & op, ST *scalar,
bool *transp);
921#ifdef TEKO_DEFINED_MSC_EXTENSIONS
922#undef _MSC_EXTENSIONS
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
MultiVector datacopy(const MultiVector &src, MultiVector &dst)
Copy the contents of a multivector to a destination vector.
void beginBlockFill(BlockedLinearOp &blo, int rowCnt, int colCnt)
Let the blocked operator know that you are going to set the sub blocks.
@ 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.
const BlockedMultiVector toBlockedMultiVector(const MultiVector &bmv)
Convert to a BlockedMultiVector from a MultiVector.
LinearOp toLinearOp(BlockedLinearOp &blo)
Convert to a LinearOp from a BlockedLinearOp.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
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.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
VectorSpace rangeSpace(const LinearOp &lo)
Replace nonzeros with a scalar value, used to zero out an operator.
MultiVector copyAndInit(const MultiVector &v, double scalar)
Perform a deep copy of the vector.
LinearOp zero(const VectorSpace &vs)
Build a square zero operator from a single vector space.
MultiVector deepcopy(const MultiVector &v)
Perform a deep copy of the vector.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
LinearOp adjoint(ModifiableLinearOp &a)
Construct an implicit adjoint of the linear operators.
BlockedLinearOp toBlockedLinearOp(LinearOp &clo)
Converts a LinearOp to a BlockedLinearOp.
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
VectorSpace domainSpace(const LinearOp &lo)
Get the domain space of a linear operator.