Teko Version of the Day
Teko_Utilities.hpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00055 #ifndef __Teko_Utilities_hpp__
00056 #define __Teko_Utilities_hpp__
00057 
00058 #include "Epetra_CrsMatrix.h"
00059 
00060 // Teuchos includes
00061 #include "Teuchos_VerboseObject.hpp"
00062 
00063 // Thyra includes
00064 #include "Thyra_LinearOpBase.hpp"
00065 #include "Thyra_PhysicallyBlockedLinearOpBase.hpp"
00066 #include "Thyra_ProductVectorSpaceBase.hpp"
00067 #include "Thyra_VectorSpaceBase.hpp"
00068 #include "Thyra_ProductMultiVectorBase.hpp"
00069 #include "Thyra_MultiVectorStdOps.hpp"
00070 #include "Thyra_MultiVectorBase.hpp"
00071 #include "Thyra_VectorBase.hpp"
00072 #include "Thyra_VectorStdOps.hpp"
00073 #include "Thyra_DefaultBlockedLinearOp.hpp"
00074 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00075 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00076 #include "Thyra_DefaultAddedLinearOp.hpp"
00077 #include "Thyra_DefaultIdentityLinearOp.hpp"
00078 #include "Thyra_DefaultZeroLinearOp.hpp"
00079 
00080 // #define Teko_DEBUG_OFF
00081 #define Teko_DEBUG_INT 5
00082 
00083 namespace Teko {
00084 
00085 using Thyra::multiply;
00086 using Thyra::scale;
00087 using Thyra::add;
00088 using Thyra::identity;
00089 using Thyra::zero; // make it to take one argument (square matrix)
00090 using Thyra::block2x2;
00091 using Thyra::block2x1;
00092 using Thyra::block1x2;
00093 
00112 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil);
00113 
00136 Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil);
00137 
00144 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
00145 // inline const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
00146 // { return Teuchos::VerboseObjectBase::getDefaultOStream(); }
00147 
00148 #ifndef Teko_DEBUG_OFF
00149    #define Teko_DEBUG_EXPR(str) str
00150    #define Teko_DEBUG_MSG(str,level) if(level<=Teko_DEBUG_INT) { \
00151       Teuchos::RCP<Teuchos::FancyOStream> out = Teko::getOutputStream(); \
00152       *out << "Teko: " << str << std::endl; }
00153    #define Teko_DEBUG_MSG_BEGIN(level) if(level<=Teko_DEBUG_INT) { \
00154       Teko::getOutputStream()->pushTab(3); \
00155       *Teko::getOutputStream() << "Teko: Begin debug MSG\n"; \
00156       std::ostream & DEBUG_STREAM = *Teko::getOutputStream(); \
00157       Teko::getOutputStream()->pushTab(3);
00158    #define Teko_DEBUG_MSG_END() Teko::getOutputStream()->popTab(); \
00159                              *Teko::getOutputStream() << "Teko: End debug MSG\n"; \
00160                               Teko::getOutputStream()->popTab(); }
00161    #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3)
00162    #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab()
00163 
00164    struct __DebugScope__ {
00165       __DebugScope__(const std::string & str,int level)
00166          : str_(str), level_(level)
00167       { Teko_DEBUG_MSG("BEGIN "+str_,level_); Teko_DEBUG_PUSHTAB(); }      
00168       ~__DebugScope__()
00169       { Teko_DEBUG_POPTAB(); Teko_DEBUG_MSG("END "+str_,level_); } 
00170       std::string str_; int level_; };
00171    #define Teko_DEBUG_SCOPE(str,level) __DebugScope__ __dbgScope__(str,level);
00172 #else 
00173    #define Teko_DEBUG_EXPR(str)
00174    #define Teko_DEBUG_MSG(str,level)
00175    #define Teko_DEBUG_MSG_BEGIN(level) if(false) { \
00176       std::ostream & DEBUG_STREAM = *Teko::getOutputStream();
00177    #define Teko_DEBUG_MSG_END() }
00178    #define Teko_DEBUG_PUSHTAB() 
00179    #define Teko_DEBUG_POPTAB() 
00180    #define Teko_DEBUG_SCOPE(str,level)
00181 #endif
00182 
00183 // typedefs for increased simplicity
00184 typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
00185 
00186 // ----------------------------------------------------------------------------
00187 
00189 
00190 
00191 typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
00192 typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
00193 
00195 inline MultiVector toMultiVector(BlockedMultiVector & bmv) { return bmv; }
00196 
00198 inline const MultiVector toMultiVector(const BlockedMultiVector & bmv) { return bmv; }
00199 
00201 inline const BlockedMultiVector toBlockedMultiVector(const MultiVector & bmv) 
00202 { return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv); }
00203 
00205 inline int blockCount(const BlockedMultiVector & bmv)
00206 { return bmv->productSpace()->numBlocks(); }
00207 
00209 inline MultiVector getBlock(int i,const BlockedMultiVector & bmv)
00210 { return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i)); }
00211 
00213 inline MultiVector deepcopy(const MultiVector & v)
00214 { return v->clone_mv(); }
00215 
00217 inline MultiVector copyAndInit(const MultiVector & v,double scalar)
00218 { MultiVector mv = v->clone_mv(); Thyra::assign(mv.ptr(),scalar); return mv; }
00219 
00221 inline BlockedMultiVector deepcopy(const BlockedMultiVector & v)
00222 { return toBlockedMultiVector(v->clone_mv()); }
00223 
00237 inline MultiVector datacopy(const MultiVector & src,MultiVector & dst)
00238 { 
00239    if(dst==Teuchos::null)
00240       return deepcopy(src);
00241 
00242    bool rangeCompat = src->range()->isCompatible(*dst->range());
00243    bool domainCompat = src->domain()->isCompatible(*dst->domain());
00244  
00245    if(not (rangeCompat && domainCompat))
00246       return deepcopy(src);
00247 
00248    // perform data copy
00249    Thyra::assign<double>(dst.ptr(),*src);
00250    return dst;
00251 }
00252 
00266 inline BlockedMultiVector datacopy(const BlockedMultiVector & src,BlockedMultiVector & dst)
00267 { 
00268    if(dst==Teuchos::null)
00269       return deepcopy(src);
00270 
00271    bool rangeCompat = src->range()->isCompatible(*dst->range());
00272    bool domainCompat = src->domain()->isCompatible(*dst->domain());
00273 
00274    if(not (rangeCompat && domainCompat))
00275       return deepcopy(src);
00276 
00277    // perform data copy
00278    Thyra::assign<double>(dst.ptr(),*src);
00279    return dst;
00280 }
00281 
00283 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvs);
00284 
00295 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
00296                                  const std::vector<int> & indices,
00297                                  const VectorSpace & vs,
00298                                  double onValue=1.0, double offValue=0.0);
00299 
00301 
00302 // ----------------------------------------------------------------------------
00303 
00305 
00306 typedef Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > BlockedLinearOp;
00307 typedef Teuchos::RCP<const Thyra::LinearOpBase<double> > LinearOp;
00308 typedef Teuchos::RCP<Thyra::LinearOpBase<double> > InverseLinearOp;
00309 typedef Teuchos::RCP<Thyra::LinearOpBase<double> > ModifiableLinearOp;
00310 
00312 inline LinearOp zero(const VectorSpace & vs)
00313 { return Thyra::zero<double>(vs,vs); }
00314 
00316 inline VectorSpace rangeSpace(const LinearOp & lo)
00317 { return lo->range(); }
00318 
00320 inline VectorSpace domainSpace(const LinearOp & lo)
00321 { return lo->domain(); }
00322 
00324 inline BlockedLinearOp toBlockedLinearOp(LinearOp & clo)
00325 {
00326    Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
00327    return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
00328 }
00329 
00331 inline const BlockedLinearOp toBlockedLinearOp(const LinearOp & clo)
00332 {
00333    Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
00334    return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
00335 }
00336 
00338 inline LinearOp toLinearOp(BlockedLinearOp & blo) { return blo; }
00339 
00341 inline const LinearOp toLinearOp(const BlockedLinearOp & blo) { return blo; }
00342 
00344 inline LinearOp toLinearOp(ModifiableLinearOp & blo) { return blo; }
00345 
00347 inline const LinearOp toLinearOp(const ModifiableLinearOp & blo) { return blo; }
00348 
00350 inline int blockRowCount(const BlockedLinearOp & blo)
00351 { return blo->productRange()->numBlocks(); }
00352 
00354 inline int blockColCount(const BlockedLinearOp & blo)
00355 { return blo->productDomain()->numBlocks(); }
00356 
00358 inline LinearOp getBlock(int i,int j,const BlockedLinearOp & blo)
00359 { return blo->getBlock(i,j); }
00360 
00362 inline void setBlock(int i,int j,BlockedLinearOp & blo, const LinearOp & lo)
00363 { return blo->setBlock(i,j,lo); }
00364 
00366 inline BlockedLinearOp createBlockedOp()
00367 { return rcp(new Thyra::DefaultBlockedLinearOp<double>()); }
00368 
00380 inline void beginBlockFill(BlockedLinearOp & blo,int rowCnt,int colCnt)
00381 { blo->beginBlockFill(rowCnt,colCnt); }
00382 
00392 inline void beginBlockFill(BlockedLinearOp & blo)
00393 { blo->beginBlockFill(); }
00394 
00396 inline void endBlockFill(BlockedLinearOp & blo)
00397 { blo->endBlockFill(); }
00398 
00400 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
00401 
00403 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
00404 
00424 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo);
00425 
00427 bool isZeroOp(const LinearOp op);
00428 
00437 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op);
00438 
00447 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op);
00448 
00456 ModifiableLinearOp getLumpedMatrix(const LinearOp & op);
00457 
00466 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op);
00467 
00469 
00471 
00472 
00491 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
00492 
00511 inline void applyOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
00512 { const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
00513   applyOp(A,x_mv,y_mv,alpha,beta); }
00514 
00524 void update(double alpha,const MultiVector & x,double beta,MultiVector & y);
00525 
00527 inline void update(double alpha,const BlockedMultiVector & x,double beta,BlockedMultiVector & y)
00528 { MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
00529   update(alpha,x_mv,beta,y_mv); }
00530 
00532 inline void scale(const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
00533 
00535 inline void scale(const double alpha,BlockedMultiVector & x) 
00536 {  MultiVector x_mv = toMultiVector(x); scale(alpha,x_mv); }
00537 
00539 inline LinearOp scale(const double alpha,ModifiableLinearOp & a) 
00540 {  return Thyra::nonconstScale(alpha,a); }
00541 
00543 inline LinearOp adjoint(ModifiableLinearOp & a) 
00544 {  return Thyra::nonconstAdjoint(a); }
00545 
00547 
00549 
00550 
00562 const ModifiableLinearOp getDiagonalOp(const LinearOp & op);
00563 
00569 const MultiVector getDiagonal(const LinearOp & op);
00570 
00582 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op);
00583 
00596 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr);
00597 
00612 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
00613                                           const ModifiableLinearOp & destOp);
00614 
00625 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr);
00626 
00640 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
00641                                           const ModifiableLinearOp & destOp);
00642 
00653 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr);
00654 
00667 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
00668                                      const ModifiableLinearOp & destOp);
00669 
00673 const LinearOp explicitTranspose(const LinearOp & op);
00674 
00678 const LinearOp buildDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
00679 
00683 const LinearOp buildInvDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
00684 
00686 
00710 double computeSpectralRad(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,double tol,
00711                           bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
00712 
00736 double computeSmallestMagEig(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A, double tol,
00737                           bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
00738 
00740 typedef enum {  Diagonal     
00741               , Lumped       
00742               , AbsRowSum    
00743         , BlkDiag      
00744               , NotDiag      
00745               } DiagonalType;
00746 
00755 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
00756 
00765 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
00766 
00772 const MultiVector getDiagonal(const LinearOp & op,const DiagonalType & dt);
00773 
00780 std::string getDiagonalName(const DiagonalType & dt);
00781 
00790 DiagonalType getDiagonalType(std::string name);
00791 
00792 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp & Op);
00793 
00796 double norm_1(const MultiVector & v,std::size_t col);
00797 
00800 double norm_2(const MultiVector & v,std::size_t col);
00801 
00802 } // end namespace Teko
00803 
00804 #endif
 All Classes Files Functions Variables