Teko_Utilities.hpp

Go to the documentation of this file.

/*
// @HEADER
// 
// ***********************************************************************
// 
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 Sandia Corporation 
//  
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//  
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//  
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//  
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//  
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission. 
//  
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//  
// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
// 
// ***********************************************************************
// 
// @HEADER
 
*/
 
#ifndef __Teko_Utilities_hpp__
#define __Teko_Utilities_hpp__
 
#include "Epetra_CrsMatrix.h"
#include "Tpetra_CrsMatrix.hpp"
 
// Teuchos includes
#include "Teuchos_VerboseObject.hpp"
 
// Thyra includes
#include "Thyra_LinearOpBase.hpp"
#include "Thyra_PhysicallyBlockedLinearOpBase.hpp"
#include "Thyra_ProductVectorSpaceBase.hpp"
#include "Thyra_VectorSpaceBase.hpp"
#include "Thyra_ProductMultiVectorBase.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
#include "Thyra_MultiVectorBase.hpp"
#include "Thyra_VectorBase.hpp"
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_DefaultBlockedLinearOp.hpp"
#include "Thyra_DefaultMultipliedLinearOp.hpp"
#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
#include "Thyra_DefaultAddedLinearOp.hpp"
#include "Thyra_DefaultIdentityLinearOp.hpp"
#include "Thyra_DefaultZeroLinearOp.hpp"
 
#include "Teko_ConfigDefs.hpp"
 
#ifdef _MSC_VER
#ifndef _MSC_EXTENSIONS
#define _MSC_EXTENSIONS
#define TEKO_DEFINED_MSC_EXTENSIONS
#endif
#include <iso646.h> // For C alternative tokens
#endif
 
// #define Teko_DEBUG_OFF
#define Teko_DEBUG_INT 5
 
namespace Teko {
 
using Thyra::multiply;
using Thyra::scale;
using Thyra::add;
using Thyra::identity;
using Thyra::zero; // make it to take one argument (square matrix)
using Thyra::block2x2;
using Thyra::block2x1;
using Thyra::block1x2;
 
Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil);
Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(int dim,ST * coords,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
 
Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil);
Teuchos::RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil);
 
const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
// inline const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream();
// { return Teuchos::VerboseObjectBase::getDefaultOStream(); }
 
#ifndef Teko_DEBUG_OFF
//#if 0
   #define Teko_DEBUG_EXPR(str) str
   #define Teko_DEBUG_MSG(str,level) if(level<=Teko_DEBUG_INT) { \
      Teuchos::RCP<Teuchos::FancyOStream> out = Teko::getOutputStream(); \
      *out << "Teko: " << str << std::endl; }
   #define Teko_DEBUG_MSG_BEGIN(level) if(level<=Teko_DEBUG_INT) { \
      Teko::getOutputStream()->pushTab(3); \
      *Teko::getOutputStream() << "Teko: Begin debug MSG\n"; \
      std::ostream & DEBUG_STREAM = *Teko::getOutputStream(); \
      Teko::getOutputStream()->pushTab(3);
   #define Teko_DEBUG_MSG_END() Teko::getOutputStream()->popTab(); \
                             *Teko::getOutputStream() << "Teko: End debug MSG\n"; \
                              Teko::getOutputStream()->popTab(); }
   #define Teko_DEBUG_PUSHTAB() Teko::getOutputStream()->pushTab(3)
   #define Teko_DEBUG_POPTAB() Teko::getOutputStream()->popTab()
   #define Teko_DEBUG_SCOPE(str,level)
 
//   struct __DebugScope__ {
//      __DebugScope__(const std::string & str,int level)
//         : str_(str), level_(level)
//      { Teko_DEBUG_MSG("BEGIN "+str_,level_); Teko_DEBUG_PUSHTAB(); }      
//      ~__DebugScope__()
//      { Teko_DEBUG_POPTAB(); Teko_DEBUG_MSG("END "+str_,level_); } 
//      std::string str_; int level_; };
//   #define Teko_DEBUG_SCOPE(str,level) __DebugScope__ __dbgScope__(str,level);
#else 
   #define Teko_DEBUG_EXPR(str)
   #define Teko_DEBUG_MSG(str,level)
   #define Teko_DEBUG_MSG_BEGIN(level) if(false) { \
      std::ostream & DEBUG_STREAM = *Teko::getOutputStream();
   #define Teko_DEBUG_MSG_END() }
   #define Teko_DEBUG_PUSHTAB() 
   #define Teko_DEBUG_POPTAB() 
   #define Teko_DEBUG_SCOPE(str,level)
#endif
 
// typedefs for increased simplicity
typedef Teuchos::RCP<const Thyra::VectorSpaceBase<double> > VectorSpace;
 
// ----------------------------------------------------------------------------
 
 
 
typedef Teuchos::RCP<Thyra::ProductMultiVectorBase<double> > BlockedMultiVector;
typedef Teuchos::RCP<Thyra::MultiVectorBase<double> > MultiVector;
 
inline MultiVector toMultiVector(BlockedMultiVector & bmv) { return bmv; }
 
inline const MultiVector toMultiVector(const BlockedMultiVector & bmv) { return bmv; }
 
inline const BlockedMultiVector toBlockedMultiVector(const MultiVector & bmv) 
{ return Teuchos::rcp_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(bmv); }
 
inline int blockCount(const BlockedMultiVector & bmv)
{ return bmv->productSpace()->numBlocks(); }
 
inline MultiVector getBlock(int i,const BlockedMultiVector & bmv)
{ return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(bmv->getMultiVectorBlock(i)); }
 
inline MultiVector deepcopy(const MultiVector & v)
{ return v->clone_mv(); }
 
inline MultiVector copyAndInit(const MultiVector & v,double scalar)
{ MultiVector mv = v->clone_mv(); Thyra::assign(mv.ptr(),scalar); return mv; }
 
inline BlockedMultiVector deepcopy(const BlockedMultiVector & v)
{ return toBlockedMultiVector(v->clone_mv()); }
 
inline MultiVector datacopy(const MultiVector & src,MultiVector & dst)
{ 
   if(dst==Teuchos::null)
      return deepcopy(src);
 
   bool rangeCompat = src->range()->isCompatible(*dst->range());
   bool domainCompat = src->domain()->isCompatible(*dst->domain());
 
   if(not (rangeCompat && domainCompat))
      return deepcopy(src);
 
   // perform data copy
   Thyra::assign<double>(dst.ptr(),*src);
   return dst;
}
 
inline BlockedMultiVector datacopy(const BlockedMultiVector & src,BlockedMultiVector & dst)
{ 
   if(dst==Teuchos::null)
      return deepcopy(src);
 
   bool rangeCompat = src->range()->isCompatible(*dst->range());
   bool domainCompat = src->domain()->isCompatible(*dst->domain());
 
   if(not (rangeCompat && domainCompat))
      return deepcopy(src);
 
   // perform data copy
   Thyra::assign<double>(dst.ptr(),*src);
   return dst;
}
 
BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvs);
 
Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
                                 const std::vector<int> & indices,
                                 const VectorSpace & vs,
                                 double onValue=1.0, double offValue=0.0);
 
 
// ----------------------------------------------------------------------------
 
 
typedef Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > BlockedLinearOp;
typedef Teuchos::RCP<const Thyra::LinearOpBase<ST> > LinearOp;
typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > InverseLinearOp;
typedef Teuchos::RCP<Thyra::LinearOpBase<ST> > ModifiableLinearOp;
 
inline LinearOp zero(const VectorSpace & vs)
{ return Thyra::zero<ST>(vs,vs); }
 
void putScalar(const ModifiableLinearOp & op,double scalar);
 
inline VectorSpace rangeSpace(const LinearOp & lo)
{ return lo->range(); }
 
inline VectorSpace domainSpace(const LinearOp & lo)
{ return lo->domain(); }
 
inline BlockedLinearOp toBlockedLinearOp(LinearOp & clo)
{
   Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
   return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
}
 
inline const BlockedLinearOp toBlockedLinearOp(const LinearOp & clo)
{
   Teuchos::RCP<Thyra::LinearOpBase<double> > lo = Teuchos::rcp_const_cast<Thyra::LinearOpBase<double> >(clo);
   return Teuchos::rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> > (lo);
}
 
inline LinearOp toLinearOp(BlockedLinearOp & blo) { return blo; }
 
inline const LinearOp toLinearOp(const BlockedLinearOp & blo) { return blo; }
 
inline LinearOp toLinearOp(ModifiableLinearOp & blo) { return blo; }
 
inline const LinearOp toLinearOp(const ModifiableLinearOp & blo) { return blo; }
 
inline int blockRowCount(const BlockedLinearOp & blo)
{ return blo->productRange()->numBlocks(); }
 
inline int blockColCount(const BlockedLinearOp & blo)
{ return blo->productDomain()->numBlocks(); }
 
inline LinearOp getBlock(int i,int j,const BlockedLinearOp & blo)
{ return blo->getBlock(i,j); }
 
inline void setBlock(int i,int j,BlockedLinearOp & blo, const LinearOp & lo)
{ return blo->setBlock(i,j,lo); }
 
inline BlockedLinearOp createBlockedOp()
{ return rcp(new Thyra::DefaultBlockedLinearOp<double>()); }
 
inline void beginBlockFill(BlockedLinearOp & blo,int rowCnt,int colCnt)
{ blo->beginBlockFill(rowCnt,colCnt); }
 
inline void beginBlockFill(BlockedLinearOp & blo)
{ blo->beginBlockFill(); }
 
inline void endBlockFill(BlockedLinearOp & blo)
{ blo->endBlockFill(); }
 
BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
 
BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill=true);
 
BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo);
 
bool isZeroOp(const LinearOp op);
 
ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op);
 
ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op);
 
ModifiableLinearOp getLumpedMatrix(const LinearOp & op);
 
ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op);
 
 
 
 
void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
 
 
void applyTransposeOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha=1.0,double beta=0.0);
 
inline void applyOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
{ const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
  applyOp(A,x_mv,y_mv,alpha,beta); }
 
inline void applyTransposeOp(const LinearOp & A,const BlockedMultiVector & x,BlockedMultiVector & y,double alpha=1.0,double beta=0.0)
{ const MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
  applyTransposeOp(A,x_mv,y_mv,alpha,beta); }
 
void update(double alpha,const MultiVector & x,double beta,MultiVector & y);
 
inline void update(double alpha,const BlockedMultiVector & x,double beta,BlockedMultiVector & y)
{ MultiVector x_mv = toMultiVector(x); MultiVector y_mv = toMultiVector(y);
  update(alpha,x_mv,beta,y_mv); }
 
inline void scale(const double alpha,MultiVector & x) { Thyra::scale<double>(alpha,x.ptr()); }
 
inline void scale(const double alpha,BlockedMultiVector & x) 
{  MultiVector x_mv = toMultiVector(x); scale(alpha,x_mv); }
 
inline LinearOp scale(const double alpha,ModifiableLinearOp & a) 
{  return Thyra::nonconstScale(alpha,a); }
 
inline LinearOp adjoint(ModifiableLinearOp & a) 
{  return Thyra::nonconstAdjoint(a); }
 
 
 
 
const ModifiableLinearOp getDiagonalOp(const LinearOp & op);
 
const MultiVector getDiagonal(const LinearOp & op);
 
const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op);
 
const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr);
 
const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
                                          const ModifiableLinearOp & destOp);
 
const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr);
 
const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
                                          const ModifiableLinearOp & destOp);
 
const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr);
 
const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
                                     const ModifiableLinearOp & destOp);
 
const ModifiableLinearOp explicitSum(const LinearOp & opl,
                                     const ModifiableLinearOp & destOp);
 
const LinearOp explicitTranspose(const LinearOp & op);
 
double frobeniusNorm(const LinearOp & op);
double oneNorm(const LinearOp & op);
double infNorm(const LinearOp & op);
 
const LinearOp buildDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
 
const LinearOp buildInvDiagonal(const MultiVector & v,const std::string & lbl="ANYM");
 
 
double computeSpectralRad(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,double tol,
                          bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
 
double computeSmallestMagEig(const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A, double tol,
                          bool isHermitian=false,int numBlocks=5,int restart=0,int verbosity=0);
 
typedef enum {  Diagonal     
              , Lumped       
              , AbsRowSum    
        , BlkDiag      
              , NotDiag      
              } DiagonalType;
 
ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
 
ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt);
 
const MultiVector getDiagonal(const LinearOp & op,const DiagonalType & dt);
 
std::string getDiagonalName(const DiagonalType & dt);
 
DiagonalType getDiagonalType(std::string name);
 
LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp & Op);
 
double norm_1(const MultiVector & v,std::size_t col);
 
double norm_2(const MultiVector & v,std::size_t col);
 
void clipLower(MultiVector & v,double lowerBound);
 
void clipUpper(MultiVector & v,double upperBound);
 
void replaceValue(MultiVector & v,double currentValue,double newValue);
 
void columnAverages(const MultiVector & v,std::vector<double> & averages);
 
double average(const MultiVector & v);
 
bool isPhysicallyBlockedLinearOp(const LinearOp & op);
 
Teuchos::RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(const LinearOp & op, ST *scalar, bool *transp);
 
} // end namespace Teko
 
#ifdef _MSC_VER
#ifdef TEKO_DEFINED_MSC_EXTENSIONS
#undef _MSC_EXTENSIONS
#endif
#endif
 
#endif