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