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 void putScalar(const ModifiableLinearOp & op,double scalar);
00317 
00319 inline VectorSpace rangeSpace(const LinearOp & lo)
00320 { return lo->range(); }
00321 
00323 inline VectorSpace domainSpace(const LinearOp & lo)
00324 { return lo->domain(); }
00325 
00327 inline BlockedLinearOp toBlockedLinearOp(LinearOp & clo)
00328 {
00329    Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
00330    return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
00331 }
00332 
00334 inline const BlockedLinearOp toBlockedLinearOp(const LinearOp & clo)
00335 {
00336    Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
00337    return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
00338 }
00339 
00341 inline LinearOp toLinearOp(BlockedLinearOp & blo) { return blo; }
00342 
00344 inline const LinearOp toLinearOp(const BlockedLinearOp & blo) { return blo; }
00345 
00347 inline LinearOp toLinearOp(ModifiableLinearOp & blo) { return blo; }
00348 
00350 inline const LinearOp toLinearOp(const ModifiableLinearOp & blo) { return blo; }
00351 
00353 inline int blockRowCount(const BlockedLinearOp & blo)
00354 { return blo->productRange()->numBlocks(); }
00355 
00357 inline int blockColCount(const BlockedLinearOp & blo)
00358 { return blo->productDomain()->numBlocks(); }
00359 
00361 inline LinearOp getBlock(int i,int j,const BlockedLinearOp & blo)
00362 { return blo->getBlock(i,j); }
00363 
00365 inline void setBlock(int i,int j,BlockedLinearOp & blo, const LinearOp & lo)
00366 { return blo->setBlock(i,j,lo); }
00367 
00369 inline BlockedLinearOp createBlockedOp()
00370 { return rcp(new Thyra::DefaultBlockedLinearOp<double>()); }
00371 
00383 inline void beginBlockFill(BlockedLinearOp & blo,int rowCnt,int colCnt)
00384 { blo->beginBlockFill(rowCnt,colCnt); }
00385 
00395 inline void beginBlockFill(BlockedLinearOp & blo)
00396 { blo->beginBlockFill(); }
00397 
00399 inline void endBlockFill(BlockedLinearOp & blo)
00400 { blo->endBlockFill(); }
00401 
00403 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
00404 
00406 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
00407 
00427 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo);
00428 
00430 bool isZeroOp(const LinearOp op);
00431 
00440 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op);
00441 
00450 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op);
00451 
00459 ModifiableLinearOp getLumpedMatrix(const LinearOp & op);
00460 
00469 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op);
00470 
00472 
00474 
00475 
00494 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
00495 
00514 inline void applyOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
00515 { const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
00516   applyOp(A,x_mv,y_mv,alpha,beta); }
00517 
00527 void update(double alpha,const MultiVector & x,double beta,MultiVector & y);
00528 
00530 inline void update(double alpha,const BlockedMultiVector & x,double beta,BlockedMultiVector & y)
00531 { MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
00532   update(alpha,x_mv,beta,y_mv); }
00533 
00535 inline void scale(const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
00536 
00538 inline void scale(const double alpha,BlockedMultiVector & x) 
00539 {  MultiVector x_mv = toMultiVector(x); scale(alpha,x_mv); }
00540 
00542 inline LinearOp scale(const double alpha,ModifiableLinearOp & a) 
00543 {  return Thyra::nonconstScale(alpha,a); }
00544 
00546 inline LinearOp adjoint(ModifiableLinearOp & a) 
00547 {  return Thyra::nonconstAdjoint(a); }
00548 
00550 
00552 
00553 
00565 const ModifiableLinearOp getDiagonalOp(const LinearOp & op);
00566 
00572 const MultiVector getDiagonal(const LinearOp & op);
00573 
00585 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op);
00586 
00599 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr);
00600 
00615 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
00616                                           const ModifiableLinearOp & destOp);
00617 
00628 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr);
00629 
00643 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
00644                                           const ModifiableLinearOp & destOp);
00645 
00656 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr);
00657 
00670 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
00671                                      const ModifiableLinearOp & destOp);
00672 
00675 const ModifiableLinearOp explicitSum(const LinearOp & opl,
00676                                      const ModifiableLinearOp & destOp);
00677 
00681 const LinearOp explicitTranspose(const LinearOp & op);
00682 
00686 const LinearOp buildDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
00687 
00691 const LinearOp buildInvDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
00692 
00694 
00718 double computeSpectralRad(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,double tol,
00719                           bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
00720 
00744 double computeSmallestMagEig(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A, double tol,
00745                           bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
00746 
00748 typedef enum {  Diagonal     
00749               , Lumped       
00750               , AbsRowSum    
00751         , BlkDiag      
00752               , NotDiag      
00753               } DiagonalType;
00754 
00763 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
00764 
00773 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
00774 
00780 const MultiVector getDiagonal(const LinearOp & op,const DiagonalType & dt);
00781 
00788 std::string getDiagonalName(const DiagonalType & dt);
00789 
00798 DiagonalType getDiagonalType(std::string name);
00799 
00800 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp & Op);
00801 
00804 double norm_1(const MultiVector & v,std::size_t col);
00805 
00808 double norm_2(const MultiVector & v,std::size_t col);
00809 
00813 void clipLower(MultiVector & v,double lowerBound);
00814 
00818 void clipUpper(MultiVector & v,double upperBound);
00819 
00823 void replaceValue(MultiVector & v,double currentValue,double newValue);
00824 
00827 void columnAverages(const MultiVector & v,std::vector<double> & averages);
00828 
00831 double average(const MultiVector & v);
00832 
00833 } // end namespace Teko
00834 
00835 #endif
 All Classes Files Functions Variables