Teko Version of the Day
Loading...
Searching...
No Matches
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
97namespace Teko {
98
99using Thyra::multiply;
100using Thyra::scale;
101using Thyra::add;
102using Thyra::identity;
103using Thyra::zero; // make it to take one argument (square matrix)
104using Thyra::block2x2;
105using Thyra::block2x1;
106using Thyra::block1x2;
107
126#ifdef TEKO_HAVE_EPETRA
127Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil);
128#endif
129
130Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(int dim,ST * coords,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
131
154#ifdef TEKO_HAVE_EPETRA
155Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil);
156#endif
157
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);
159
166const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
167// inline const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
168// { return Teuchos::VerboseObjectBase::getDefaultOStream(); }
169
170#ifndef Teko_DEBUG_OFF
171//#if 0
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)
187
188// struct __DebugScope__ {
189// __DebugScope__(const std::string & str,int level)
190// : str_(str), level_(level)
191// { Teko_DEBUG_MSG("BEGIN "+str_,level_); Teko_DEBUG_PUSHTAB(); }
192// ~__DebugScope__()
193// { Teko_DEBUG_POPTAB(); Teko_DEBUG_MSG("END "+str_,level_); }
194// std::string str_; int level_; };
195// #define Teko_DEBUG_SCOPE(str,level) __DebugScope__ __dbgScope__(str,level);
196#else
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)
205#endif
206
207// typedefs for increased simplicity
208typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
209
210// ----------------------------------------------------------------------------
211
213
214
215typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
216typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
217
219inline MultiVector toMultiVector(BlockedMultiVector & bmv) { return bmv; }
220
222inline const MultiVector toMultiVector(const BlockedMultiVector & bmv) { return bmv; }
223
225inline const BlockedMultiVector toBlockedMultiVector(const MultiVector & bmv)
226{ return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv); }
227
229inline int blockCount(const BlockedMultiVector & bmv)
230{ return bmv->productSpace()->numBlocks(); }
231
233inline MultiVector getBlock(int i,const BlockedMultiVector & bmv)
234{ return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i)); }
235
237inline MultiVector deepcopy(const MultiVector & v)
238{ return v->clone_mv(); }
239
241inline MultiVector copyAndInit(const MultiVector & v,double scalar)
242{ MultiVector mv = v->clone_mv(); Thyra::assign(mv.ptr(),scalar); return mv; }
243
245inline BlockedMultiVector deepcopy(const BlockedMultiVector & v)
246{ return toBlockedMultiVector(v->clone_mv()); }
247
261inline MultiVector datacopy(const MultiVector & src,MultiVector & dst)
262{
263 if(dst==Teuchos::null)
264 return deepcopy(src);
265
266 bool rangeCompat = src->range()->isCompatible(*dst->range());
267 bool domainCompat = src->domain()->isCompatible(*dst->domain());
268
269 if(not (rangeCompat && domainCompat))
270 return deepcopy(src);
271
272 // perform data copy
273 Thyra::assign<double>(dst.ptr(),*src);
274 return dst;
275}
276
290inline BlockedMultiVector datacopy(const BlockedMultiVector & src,BlockedMultiVector & dst)
291{
292 if(dst==Teuchos::null)
293 return deepcopy(src);
294
295 bool rangeCompat = src->range()->isCompatible(*dst->range());
296 bool domainCompat = src->domain()->isCompatible(*dst->domain());
297
298 if(not (rangeCompat && domainCompat))
299 return deepcopy(src);
300
301 // perform data copy
302 Thyra::assign<double>(dst.ptr(),*src);
303 return dst;
304}
305
307BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvs);
308
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);
323
325
326// ----------------------------------------------------------------------------
327
329
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;
334
336inline LinearOp zero(const VectorSpace & vs)
337{ return Thyra::zero<ST>(vs,vs); }
338
340#ifdef TEKO_HAVE_EPETRA
341void putScalar(const ModifiableLinearOp & op,double scalar);
342#endif
343
345inline VectorSpace rangeSpace(const LinearOp & lo)
346{ return lo->range(); }
347
349inline VectorSpace domainSpace(const LinearOp & lo)
350{ return lo->domain(); }
351
353inline BlockedLinearOp toBlockedLinearOp(LinearOp & clo)
354{
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);
357}
358
360inline const BlockedLinearOp toBlockedLinearOp(const LinearOp & clo)
361{
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);
364}
365
367inline LinearOp toLinearOp(BlockedLinearOp & blo) { return blo; }
368
370inline const LinearOp toLinearOp(const BlockedLinearOp & blo) { return blo; }
371
373inline LinearOp toLinearOp(ModifiableLinearOp & blo) { return blo; }
374
376inline const LinearOp toLinearOp(const ModifiableLinearOp & blo) { return blo; }
377
379inline int blockRowCount(const BlockedLinearOp & blo)
380{ return blo->productRange()->numBlocks(); }
381
383inline int blockColCount(const BlockedLinearOp & blo)
384{ return blo->productDomain()->numBlocks(); }
385
387inline LinearOp getBlock(int i,int j,const BlockedLinearOp & blo)
388{ return blo->getBlock(i,j); }
389
391inline void setBlock(int i,int j,BlockedLinearOp & blo, const LinearOp & lo)
392{ return blo->setBlock(i,j,lo); }
393
395inline BlockedLinearOp createBlockedOp()
396{ return rcp(new Thyra::DefaultBlockedLinearOp<double>()); }
397
409inline void beginBlockFill(BlockedLinearOp & blo,int rowCnt,int colCnt)
410{ blo->beginBlockFill(rowCnt,colCnt); }
411
421inline void beginBlockFill(BlockedLinearOp & blo)
422{ blo->beginBlockFill(); }
423
425inline void endBlockFill(BlockedLinearOp & blo)
426{ blo->endBlockFill(); }
427
429BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
430
432BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
433
453BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo);
454
456bool isZeroOp(const LinearOp op);
457
466ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op);
467
476ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op);
477
485ModifiableLinearOp getLumpedMatrix(const LinearOp & op);
486
495ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op);
496
498
500
501
520void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
521
522
541void applyTransposeOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
542
561inline void applyOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
562{ const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
563 applyOp(A,x_mv,y_mv,alpha,beta); }
564
583inline void applyTransposeOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
584{ const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
585 applyTransposeOp(A,x_mv,y_mv,alpha,beta); }
586
596void update(double alpha,const MultiVector & x,double beta,MultiVector & y);
597
599inline void update(double alpha,const BlockedMultiVector & x,double beta,BlockedMultiVector & y)
600{ MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
601 update(alpha,x_mv,beta,y_mv); }
602
604inline void scale(const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
605
607inline void scale(const double alpha,BlockedMultiVector & x)
608{ MultiVector x_mv = toMultiVector(x); scale(alpha,x_mv); }
609
611inline LinearOp scale(const double alpha,ModifiableLinearOp & a)
612{ return Thyra::nonconstScale(alpha,a); }
613
615inline LinearOp adjoint(ModifiableLinearOp & a)
616{ return Thyra::nonconstAdjoint(a); }
617
619
621
622
634const ModifiableLinearOp getDiagonalOp(const LinearOp & op);
635
641const MultiVector getDiagonal(const LinearOp & op);
642
654const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op);
655
668const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr);
669
684const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
685 const ModifiableLinearOp & destOp);
686
697const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr);
698
712const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
713 const ModifiableLinearOp & destOp);
714
725const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr);
726
739const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
740 const ModifiableLinearOp & destOp);
741
744const ModifiableLinearOp explicitSum(const LinearOp & opl,
745 const ModifiableLinearOp & destOp);
746
750const LinearOp explicitTranspose(const LinearOp & op);
751
754double frobeniusNorm(const LinearOp & op);
755double oneNorm(const LinearOp & op);
756double infNorm(const LinearOp & op);
757
761const LinearOp buildDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
762
766const LinearOp buildInvDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
767
769
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);
795
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);
821
823typedef enum { Diagonal
828 } DiagonalType;
829
838ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
839
848ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
849
855const MultiVector getDiagonal(const LinearOp & op,const DiagonalType & dt);
856
863std::string getDiagonalName(const DiagonalType & dt);
864
873DiagonalType getDiagonalType(std::string name);
874
875#ifdef TEKO_HAVE_EPETRA
876LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp & Op);
877#endif
878
881double norm_1(const MultiVector & v,std::size_t col);
882
885double norm_2(const MultiVector & v,std::size_t col);
886
890void clipLower(MultiVector & v,double lowerBound);
891
895void clipUpper(MultiVector & v,double upperBound);
896
900void replaceValue(MultiVector & v,double currentValue,double newValue);
901
904void columnAverages(const MultiVector & v,std::vector<double> & averages);
905
908double average(const MultiVector & v);
909
912bool isPhysicallyBlockedLinearOp(const LinearOp & op);
913
916Teuchos::RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(const LinearOp & op, ST *scalar, bool *transp);
917
918} // end namespace Teko
919
920#ifdef _MSC_VER
921#ifdef TEKO_DEFINED_MSC_EXTENSIONS
922#undef _MSC_EXTENSIONS
923#endif
924#endif
925
926#endif
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.