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
100 using Thyra::block1x2;
101 using Thyra::block2x1;
102 using Thyra::block2x2;
103 using Thyra::identity;
104 using Thyra::multiply;
126 #ifdef TEKO_HAVE_EPETRA
127 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
int dim,
double *coords,
128 const Epetra_CrsMatrix &stencil);
131 Teuchos::RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildGraphLaplacian(
132 int dim, ST *coords,
const Tpetra::CrsMatrix<ST, LO, GO, NT> &stencil);
156 #ifdef TEKO_HAVE_EPETRA
157 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(
double *x,
double *y,
double *z,
int stride,
158 const Epetra_CrsMatrix &stencil);
161 Teuchos::RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildGraphLaplacian(
162 ST *x, ST *y, ST *z, GO stride,
const Tpetra::CrsMatrix<ST, LO, GO, NT> &stencil);
170 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
174 #ifndef Teko_DEBUG_OFF
176 #define Teko_DEBUG_EXPR(str) str
177 #define Teko_DEBUG_MSG(str, level) \
178 if (level <= Teko_DEBUG_INT) { \
179 Teuchos::RCP<Teuchos::FancyOStream> out = Teko::getOutputStream(); \
180 *out << "Teko: " << str << std::endl; \
182 #define Teko_DEBUG_MSG_BEGIN(level) \
183 if (level <= Teko_DEBUG_INT) { \
184 Teko::getOutputStream()->pushTab(3); \
185 *Teko::getOutputStream() << "Teko: Begin debug MSG\n"; \
186 std::ostream &DEBUG_STREAM = *Teko::getOutputStream(); \
187 Teko::getOutputStream()->pushTab(3);
188 #define Teko_DEBUG_MSG_END() \
189 Teko::getOutputStream()->popTab(); \
190 *Teko::getOutputStream() << "Teko: End debug MSG\n"; \
191 Teko::getOutputStream()->popTab(); \
193 #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3)
194 #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab()
195 #define Teko_DEBUG_SCOPE(str, level)
206 #define Teko_DEBUG_EXPR(str)
207 #define Teko_DEBUG_MSG(str, level)
208 #define Teko_DEBUG_MSG_BEGIN(level) \
210 std::ostream &DEBUG_STREAM = *Teko::getOutputStream();
211 #define Teko_DEBUG_MSG_END() }
212 #define Teko_DEBUG_PUSHTAB()
213 #define Teko_DEBUG_POPTAB()
214 #define Teko_DEBUG_SCOPE(str, level)
218 typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
225 typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
226 typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
229 inline MultiVector toMultiVector(BlockedMultiVector &bmv) {
return bmv; }
232 inline const MultiVector toMultiVector(
const BlockedMultiVector &bmv) {
return bmv; }
235 inline const BlockedMultiVector toBlockedMultiVector(
const MultiVector &bmv) {
236 return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv);
240 inline int blockCount(
const BlockedMultiVector &bmv) {
return bmv->productSpace()->numBlocks(); }
243 inline MultiVector getBlock(
int i,
const BlockedMultiVector &bmv) {
244 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i));
248 inline MultiVector deepcopy(
const MultiVector &v) {
return v->clone_mv(); }
251 inline MultiVector copyAndInit(
const MultiVector &v,
double scalar) {
252 MultiVector mv = v->clone_mv();
253 Thyra::assign(mv.ptr(), scalar);
258 inline BlockedMultiVector deepcopy(
const BlockedMultiVector &v) {
259 return toBlockedMultiVector(v->clone_mv());
275 inline MultiVector datacopy(
const MultiVector &src, MultiVector &dst) {
276 if (dst == Teuchos::null)
return deepcopy(src);
278 bool rangeCompat = src->range()->isCompatible(*dst->range());
279 bool domainCompat = src->domain()->isCompatible(*dst->domain());
281 if (not(rangeCompat && domainCompat))
return deepcopy(src);
284 Thyra::assign<double>(dst.ptr(), *src);
301 inline BlockedMultiVector datacopy(
const BlockedMultiVector &src, BlockedMultiVector &dst) {
302 if (dst == Teuchos::null)
return deepcopy(src);
304 bool rangeCompat = src->range()->isCompatible(*dst->range());
305 bool domainCompat = src->domain()->isCompatible(*dst->domain());
307 if (not(rangeCompat && domainCompat))
return deepcopy(src);
310 Thyra::assign<double>(dst.ptr(), *src);
315 BlockedMultiVector buildBlockedMultiVector(
const std::vector<MultiVector> &mvs);
327 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
const std::vector<int> &indices,
328 const VectorSpace &vs,
329 double onValue = 1.0,
330 double offValue = 0.0);
338 typedef Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > BlockedLinearOp;
339 typedef Teuchos::RCP<const Thyra::LinearOpBase<ST> > LinearOp;
340 typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > InverseLinearOp;
341 typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > ModifiableLinearOp;
344 inline LinearOp zero(
const VectorSpace &vs) {
return Thyra::zero<ST>(vs, vs); }
347 #ifdef TEKO_HAVE_EPETRA
348 void putScalar(
const ModifiableLinearOp &op,
double scalar);
352 inline VectorSpace rangeSpace(
const LinearOp &lo) {
return lo->range(); }
355 inline VectorSpace domainSpace(
const LinearOp &lo) {
return lo->domain(); }
358 inline BlockedLinearOp toBlockedLinearOp(LinearOp &clo) {
359 Teuchos::RCP<Thyra::LinearOpBase<double> > lo =
360 Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
361 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(lo);
365 inline const BlockedLinearOp toBlockedLinearOp(
const LinearOp &clo) {
366 Teuchos::RCP<Thyra::LinearOpBase<double> > lo =
367 Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
368 return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(lo);
372 inline LinearOp toLinearOp(BlockedLinearOp &blo) {
return blo; }
375 inline const LinearOp toLinearOp(
const BlockedLinearOp &blo) {
return blo; }
378 inline LinearOp toLinearOp(ModifiableLinearOp &blo) {
return blo; }
381 inline const LinearOp toLinearOp(
const ModifiableLinearOp &blo) {
return blo; }
384 inline int blockRowCount(
const BlockedLinearOp &blo) {
return blo->productRange()->numBlocks(); }
387 inline int blockColCount(
const BlockedLinearOp &blo) {
return blo->productDomain()->numBlocks(); }
390 inline LinearOp getBlock(
int i,
int j,
const BlockedLinearOp &blo) {
return blo->getBlock(i, j); }
393 inline void setBlock(
int i,
int j, BlockedLinearOp &blo,
const LinearOp &lo) {
394 return blo->setBlock(i, j, lo);
398 inline BlockedLinearOp createBlockedOp() {
399 return rcp(
new Thyra::DefaultBlockedLinearOp<double>());
413 inline void beginBlockFill(BlockedLinearOp &blo,
int rowCnt,
int colCnt) {
414 blo->beginBlockFill(rowCnt, colCnt);
426 inline void beginBlockFill(BlockedLinearOp &blo) { blo->beginBlockFill(); }
429 inline void endBlockFill(BlockedLinearOp &blo) { blo->endBlockFill(); }
432 BlockedLinearOp getUpperTriBlocks(
const BlockedLinearOp &blo,
bool callEndBlockFill =
true);
435 BlockedLinearOp getLowerTriBlocks(
const BlockedLinearOp &blo,
bool callEndBlockFill =
true);
456 BlockedLinearOp zeroBlockedOp(
const BlockedLinearOp &blo);
459 bool isZeroOp(
const LinearOp op);
469 ModifiableLinearOp getAbsRowSumMatrix(
const LinearOp &op);
479 ModifiableLinearOp getAbsRowSumInvMatrix(
const LinearOp &op);
488 ModifiableLinearOp getLumpedMatrix(
const LinearOp &op);
498 ModifiableLinearOp getInvLumpedMatrix(
const LinearOp &op);
523 void applyOp(
const LinearOp &A,
const MultiVector &x, MultiVector &y,
double alpha = 1.0,
544 void applyTransposeOp(
const LinearOp &A,
const MultiVector &x, MultiVector &y,
double alpha = 1.0,
565 inline void applyOp(
const LinearOp &A,
const BlockedMultiVector &x, BlockedMultiVector &y,
566 double alpha = 1.0,
double beta = 0.0) {
567 const MultiVector x_mv = toMultiVector(x);
568 MultiVector y_mv = toMultiVector(y);
569 applyOp(A, x_mv, y_mv, alpha, beta);
590 inline void applyTransposeOp(
const LinearOp &A,
const BlockedMultiVector &x, BlockedMultiVector &y,
591 double alpha = 1.0,
double beta = 0.0) {
592 const MultiVector x_mv = toMultiVector(x);
593 MultiVector y_mv = toMultiVector(y);
594 applyTransposeOp(A, x_mv, y_mv, alpha, beta);
606 void update(
double alpha,
const MultiVector &x,
double beta, MultiVector &y);
609 inline void update(
double alpha,
const BlockedMultiVector &x,
double beta, BlockedMultiVector &y) {
610 MultiVector x_mv = toMultiVector(x);
611 MultiVector y_mv = toMultiVector(y);
612 update(alpha, x_mv, beta, y_mv);
616 inline void scale(
const double alpha, MultiVector &x) { Thyra::scale<double>(alpha, x.ptr()); }
619 inline void scale(
const double alpha, BlockedMultiVector &x) {
620 MultiVector x_mv = toMultiVector(x);
625 inline LinearOp scale(
const double alpha, ModifiableLinearOp &a) {
626 return Thyra::nonconstScale(alpha, a);
630 inline LinearOp adjoint(ModifiableLinearOp &a) {
return Thyra::nonconstAdjoint(a); }
648 const ModifiableLinearOp getDiagonalOp(
const LinearOp &op);
655 const MultiVector getDiagonal(
const LinearOp &op);
668 const ModifiableLinearOp getInvDiagonalOp(
const LinearOp &op);
682 const LinearOp explicitMultiply(
const LinearOp &opl,
const LinearOp &opm,
const LinearOp &opr);
698 const ModifiableLinearOp explicitMultiply(
const LinearOp &opl,
const LinearOp &opm,
699 const LinearOp &opr,
const ModifiableLinearOp &destOp);
711 const LinearOp explicitMultiply(
const LinearOp &opl,
const LinearOp &opr);
726 const ModifiableLinearOp explicitMultiply(
const LinearOp &opl,
const LinearOp &opr,
727 const ModifiableLinearOp &destOp);
739 const LinearOp explicitAdd(
const LinearOp &opl,
const LinearOp &opr);
753 const ModifiableLinearOp explicitAdd(
const LinearOp &opl,
const LinearOp &opr,
754 const ModifiableLinearOp &destOp);
758 const ModifiableLinearOp explicitSum(
const LinearOp &opl,
const ModifiableLinearOp &destOp);
763 const LinearOp explicitTranspose(
const LinearOp &op);
767 const LinearOp explicitScale(
double scalar,
const LinearOp &op);
771 double frobeniusNorm(
const LinearOp &op);
772 double oneNorm(
const LinearOp &op);
773 double infNorm(
const LinearOp &op);
778 const LinearOp buildDiagonal(
const MultiVector &v,
const std::string &lbl =
"ANYM");
783 const LinearOp buildInvDiagonal(
const MultiVector &v,
const std::string &lbl =
"ANYM");
810 double computeSpectralRad(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > &A,
double tol,
811 bool isHermitian =
false,
int numBlocks = 5,
int restart = 0,
837 double computeSmallestMagEig(
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > &A,
double tol,
838 bool isHermitian =
false,
int numBlocks = 5,
int restart = 0,
862 ModifiableLinearOp getDiagonalOp(
const Teko::LinearOp &A,
const DiagonalType &dt);
872 ModifiableLinearOp getInvDiagonalOp(
const Teko::LinearOp &A,
const DiagonalType &dt);
879 const MultiVector getDiagonal(
const LinearOp &op,
const DiagonalType &dt);
887 std::string getDiagonalName(
const DiagonalType &dt);
897 DiagonalType getDiagonalType(std::string name);
899 #ifdef TEKO_HAVE_EPETRA
900 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,
const LinearOp &Op);
905 double norm_1(
const MultiVector &v, std::size_t col);
909 double norm_2(
const MultiVector &v, std::size_t col);
914 void clipLower(MultiVector &v,
double lowerBound);
919 void clipUpper(MultiVector &v,
double upperBound);
924 void replaceValue(MultiVector &v,
double currentValue,
double newValue);
928 void columnAverages(
const MultiVector &v, std::vector<double> &averages);
932 double average(
const MultiVector &v);
936 bool isPhysicallyBlockedLinearOp(
const LinearOp &op);
940 Teuchos::RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(
941 const LinearOp &op, ST *scalar,
bool *transp);
944 std::string formatBlockName(
const std::string &prefix,
int i,
int j,
int nrow);
949 #ifdef TEKO_DEFINED_MSC_EXTENSIONS
950 #undef _MSC_EXTENSIONS