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 "Teko_ConfigDefs.hpp"

#ifdef TEKO_HAVE_EPETRA
#include "Epetra_CrsMatrix.h"
#endif

#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"

#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::add;
using Thyra::block1x2;
using Thyra::block2x1;
using Thyra::block2x2;
using Thyra::identity;
using Thyra::multiply;
using Thyra::scale;
using Thyra::zero;  // make it to take one argument (square matrix)

#ifdef TEKO_HAVE_EPETRA
Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim, double *coords,
                                                   const Epetra_CrsMatrix &stencil);
#endif

Teuchos::RCP<Tpetra::CrsMatrix<ST, LO, GO, NT> > buildGraphLaplacian(
    int dim, ST *coords, const Tpetra::CrsMatrix<ST, LO, GO, NT> &stencil);

#ifdef TEKO_HAVE_EPETRA
Teuchos::RCP<Epetra_CrsMatrix> buildGraphLaplacian(double *x, double *y, double *z, int stride,
                                                   const Epetra_CrsMatrix &stencil);
#endif

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); }

#ifdef TEKO_HAVE_EPETRA
void putScalar(const ModifiableLinearOp &op, double scalar);
#endif

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);

const LinearOp explicitScale(double scalar, 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);

#ifdef TEKO_HAVE_EPETRA
LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G, const LinearOp &Op);
#endif

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);

std::string formatBlockName(const std::string &prefix, int i, int j, int nrow);

}  // end namespace Teko

#ifdef _MSC_VER
#ifdef TEKO_DEFINED_MSC_EXTENSIONS
#undef _MSC_EXTENSIONS
#endif
#endif

#endif