Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_Utilities.hpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 //
4 // ***********************************************************************
5 //
6 // Teko: A package for block and physics based preconditioning
7 // Copyright 2010 Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40 //
41 // ***********************************************************************
42 //
43 // @HEADER
44 
45 */
46 
55 #ifndef __Teko_Utilities_hpp__
56 #define __Teko_Utilities_hpp__
57 
58 #include "Teko_ConfigDefs.hpp"
59 
60 #ifdef TEKO_HAVE_EPETRA
61 #include "Epetra_CrsMatrix.h"
62 #endif
63 
64 #include "Tpetra_CrsMatrix.hpp"
65 
66 // Teuchos includes
67 #include "Teuchos_VerboseObject.hpp"
68 
69 // Thyra includes
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"
85 
86 #ifdef _MSC_VER
87 #ifndef _MSC_EXTENSIONS
88 #define _MSC_EXTENSIONS
89 #define TEKO_DEFINED_MSC_EXTENSIONS
90 #endif
91 #include <iso646.h> // For C alternative tokens
92 #endif
93 
94 // #define Teko_DEBUG_OFF
95 #define Teko_DEBUG_INT 5
96 
97 namespace Teko {
98 
99 using Thyra::add;
100 using Thyra::block1x2;
101 using Thyra::block2x1;
102 using Thyra::block2x2;
103 using Thyra::identity;
104 using Thyra::multiply;
105 using Thyra::scale;
106 using Thyra::zero; // make it to take one argument (square matrix)
107 
126 #ifdef TEKO_HAVE_EPETRA
127 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim, double *coords,
128  const Epetra_CrsMatrix &stencil);
129 #endif
130 
131 Teuchos::RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildGraphLaplacian(
132  int dim, ST *coords, const Tpetra::CrsMatrix<ST, LO, GO, NT> &stencil);
133 
156 #ifdef TEKO_HAVE_EPETRA
157 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(double *x, double *y, double *z, int stride,
158  const Epetra_CrsMatrix &stencil);
159 #endif
160 
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);
163 
170 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
171 // inline const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
172 // { return Teuchos::VerboseObjectBase::getDefaultOStream(); }
173 
174 #ifndef Teko_DEBUG_OFF
175 //#if 0
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; \
181  }
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(); \
192  }
193 #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3)
194 #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab()
195 #define Teko_DEBUG_SCOPE(str, level)
196 
197 // struct __DebugScope__ {
198 // __DebugScope__(const std::string & str,int level)
199 // : str_(str), level_(level)
200 // { Teko_DEBUG_MSG("BEGIN "+str_,level_); Teko_DEBUG_PUSHTAB(); }
201 // ~__DebugScope__()
202 // { Teko_DEBUG_POPTAB(); Teko_DEBUG_MSG("END "+str_,level_); }
203 // std::string str_; int level_; };
204 // #define Teko_DEBUG_SCOPE(str,level) __DebugScope__ __dbgScope__(str,level);
205 #else
206 #define Teko_DEBUG_EXPR(str)
207 #define Teko_DEBUG_MSG(str, level)
208 #define Teko_DEBUG_MSG_BEGIN(level) \
209  if (false) { \
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)
215 #endif
216 
217 // typedefs for increased simplicity
218 typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
219 
220 // ----------------------------------------------------------------------------
221 
223 
224 
225 typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
226 typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
227 
229 inline MultiVector toMultiVector(BlockedMultiVector &bmv) { return bmv; }
230 
232 inline const MultiVector toMultiVector(const BlockedMultiVector &bmv) { return bmv; }
233 
235 inline const BlockedMultiVector toBlockedMultiVector(const MultiVector &bmv) {
236  return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv);
237 }
238 
240 inline int blockCount(const BlockedMultiVector &bmv) { return bmv->productSpace()->numBlocks(); }
241 
243 inline MultiVector getBlock(int i, const BlockedMultiVector &bmv) {
244  return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i));
245 }
246 
248 inline MultiVector deepcopy(const MultiVector &v) { return v->clone_mv(); }
249 
251 inline MultiVector copyAndInit(const MultiVector &v, double scalar) {
252  MultiVector mv = v->clone_mv();
253  Thyra::assign(mv.ptr(), scalar);
254  return mv;
255 }
256 
258 inline BlockedMultiVector deepcopy(const BlockedMultiVector &v) {
259  return toBlockedMultiVector(v->clone_mv());
260 }
261 
275 inline MultiVector datacopy(const MultiVector &src, MultiVector &dst) {
276  if (dst == Teuchos::null) return deepcopy(src);
277 
278  bool rangeCompat = src->range()->isCompatible(*dst->range());
279  bool domainCompat = src->domain()->isCompatible(*dst->domain());
280 
281  if (not(rangeCompat && domainCompat)) return deepcopy(src);
282 
283  // perform data copy
284  Thyra::assign<double>(dst.ptr(), *src);
285  return dst;
286 }
287 
301 inline BlockedMultiVector datacopy(const BlockedMultiVector &src, BlockedMultiVector &dst) {
302  if (dst == Teuchos::null) return deepcopy(src);
303 
304  bool rangeCompat = src->range()->isCompatible(*dst->range());
305  bool domainCompat = src->domain()->isCompatible(*dst->domain());
306 
307  if (not(rangeCompat && domainCompat)) return deepcopy(src);
308 
309  // perform data copy
310  Thyra::assign<double>(dst.ptr(), *src);
311  return dst;
312 }
313 
315 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> &mvs);
316 
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);
331 
333 
334 // ----------------------------------------------------------------------------
335 
337 
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;
342 
344 inline LinearOp zero(const VectorSpace &vs) { return Thyra::zero<ST>(vs, vs); }
345 
347 #ifdef TEKO_HAVE_EPETRA
348 void putScalar(const ModifiableLinearOp &op, double scalar);
349 #endif
350 
352 inline VectorSpace rangeSpace(const LinearOp &lo) { return lo->range(); }
353 
355 inline VectorSpace domainSpace(const LinearOp &lo) { return lo->domain(); }
356 
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);
362 }
363 
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);
369 }
370 
372 inline LinearOp toLinearOp(BlockedLinearOp &blo) { return blo; }
373 
375 inline const LinearOp toLinearOp(const BlockedLinearOp &blo) { return blo; }
376 
378 inline LinearOp toLinearOp(ModifiableLinearOp &blo) { return blo; }
379 
381 inline const LinearOp toLinearOp(const ModifiableLinearOp &blo) { return blo; }
382 
384 inline int blockRowCount(const BlockedLinearOp &blo) { return blo->productRange()->numBlocks(); }
385 
387 inline int blockColCount(const BlockedLinearOp &blo) { return blo->productDomain()->numBlocks(); }
388 
390 inline LinearOp getBlock(int i, int j, const BlockedLinearOp &blo) { return blo->getBlock(i, j); }
391 
393 inline void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo) {
394  return blo->setBlock(i, j, lo);
395 }
396 
398 inline BlockedLinearOp createBlockedOp() {
399  return rcp(new Thyra::DefaultBlockedLinearOp<double>());
400 }
401 
413 inline void beginBlockFill(BlockedLinearOp &blo, int rowCnt, int colCnt) {
414  blo->beginBlockFill(rowCnt, colCnt);
415 }
416 
426 inline void beginBlockFill(BlockedLinearOp &blo) { blo->beginBlockFill(); }
427 
429 inline void endBlockFill(BlockedLinearOp &blo) { blo->endBlockFill(); }
430 
432 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp &blo, bool callEndBlockFill = true);
433 
435 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp &blo, bool callEndBlockFill = true);
436 
456 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp &blo);
457 
459 bool isZeroOp(const LinearOp op);
460 
469 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp &op);
470 
479 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp &op);
480 
488 ModifiableLinearOp getLumpedMatrix(const LinearOp &op);
489 
498 ModifiableLinearOp getInvLumpedMatrix(const LinearOp &op);
499 
501 
503 
504 
523 void applyOp(const LinearOp &A, const MultiVector &x, MultiVector &y, double alpha = 1.0,
524  double beta = 0.0);
525 
544 void applyTransposeOp(const LinearOp &A, const MultiVector &x, MultiVector &y, double alpha = 1.0,
545  double beta = 0.0);
546 
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);
570 }
571 
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);
595 }
596 
606 void update(double alpha, const MultiVector &x, double beta, MultiVector &y);
607 
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);
613 }
614 
616 inline void scale(const double alpha, MultiVector &x) { Thyra::scale<double>(alpha, x.ptr()); }
617 
619 inline void scale(const double alpha, BlockedMultiVector &x) {
620  MultiVector x_mv = toMultiVector(x);
621  scale(alpha, x_mv);
622 }
623 
625 inline LinearOp scale(const double alpha, ModifiableLinearOp &a) {
626  return Thyra::nonconstScale(alpha, a);
627 }
628 
630 inline LinearOp adjoint(ModifiableLinearOp &a) { return Thyra::nonconstAdjoint(a); }
631 
633 
635 
636 
648 const ModifiableLinearOp getDiagonalOp(const LinearOp &op);
649 
655 const MultiVector getDiagonal(const LinearOp &op);
656 
668 const ModifiableLinearOp getInvDiagonalOp(const LinearOp &op);
669 
682 const LinearOp explicitMultiply(const LinearOp &opl, const LinearOp &opm, const LinearOp &opr);
683 
698 const ModifiableLinearOp explicitMultiply(const LinearOp &opl, const LinearOp &opm,
699  const LinearOp &opr, const ModifiableLinearOp &destOp);
700 
711 const LinearOp explicitMultiply(const LinearOp &opl, const LinearOp &opr);
712 
726 const ModifiableLinearOp explicitMultiply(const LinearOp &opl, const LinearOp &opr,
727  const ModifiableLinearOp &destOp);
728 
739 const LinearOp explicitAdd(const LinearOp &opl, const LinearOp &opr);
740 
753 const ModifiableLinearOp explicitAdd(const LinearOp &opl, const LinearOp &opr,
754  const ModifiableLinearOp &destOp);
755 
758 const ModifiableLinearOp explicitSum(const LinearOp &opl, const ModifiableLinearOp &destOp);
759 
763 const LinearOp explicitTranspose(const LinearOp &op);
764 
767 const LinearOp explicitScale(double scalar, const LinearOp &op);
768 
771 double frobeniusNorm(const LinearOp &op);
772 double oneNorm(const LinearOp &op);
773 double infNorm(const LinearOp &op);
774 
778 const LinearOp buildDiagonal(const MultiVector &v, const std::string &lbl = "ANYM");
779 
783 const LinearOp buildInvDiagonal(const MultiVector &v, const std::string &lbl = "ANYM");
784 
786 
810 double computeSpectralRad(const Teuchos::RCP<const Thyra::LinearOpBase<double> > &A, double tol,
811  bool isHermitian = false, int numBlocks = 5, int restart = 0,
812  int verbosity = 0);
813 
837 double computeSmallestMagEig(const Teuchos::RCP<const Thyra::LinearOpBase<double> > &A, double tol,
838  bool isHermitian = false, int numBlocks = 5, int restart = 0,
839  int verbosity = 0);
840 
842 typedef enum {
843  Diagonal
844  ,
845  Lumped
846  ,
847  AbsRowSum
848  ,
849  BlkDiag
850  ,
851  NotDiag
852 } DiagonalType;
853 
862 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp &A, const DiagonalType &dt);
863 
872 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp &A, const DiagonalType &dt);
873 
879 const MultiVector getDiagonal(const LinearOp &op, const DiagonalType &dt);
880 
887 std::string getDiagonalName(const DiagonalType &dt);
888 
897 DiagonalType getDiagonalType(std::string name);
898 
899 #ifdef TEKO_HAVE_EPETRA
900 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp &Op);
901 #endif
902 
905 double norm_1(const MultiVector &v, std::size_t col);
906 
909 double norm_2(const MultiVector &v, std::size_t col);
910 
914 void clipLower(MultiVector &v, double lowerBound);
915 
919 void clipUpper(MultiVector &v, double upperBound);
920 
924 void replaceValue(MultiVector &v, double currentValue, double newValue);
925 
928 void columnAverages(const MultiVector &v, std::vector<double> &averages);
929 
932 double average(const MultiVector &v);
933 
936 bool isPhysicallyBlockedLinearOp(const LinearOp &op);
937 
940 Teuchos::RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(
941  const LinearOp &op, ST *scalar, bool *transp);
942 
944 std::string formatBlockName(const std::string &prefix, int i, int j, int nrow);
945 
946 } // end namespace Teko
947 
948 #ifdef _MSC_VER
949 #ifdef TEKO_DEFINED_MSC_EXTENSIONS
950 #undef _MSC_EXTENSIONS
951 #endif
952 #endif
953 
954 #endif