Thyra_LinearOpWithSolveBaseDecl.hpp

// @HEADER
// ***********************************************************************
// 
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//                 Copyright (2004) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef THYRA_LINEAR_OP_WITH_SOLVE_BASE_DECL_HPP
#define THYRA_LINEAR_OP_WITH_SOLVE_BASE_DECL_HPP

#include "Thyra_LinearOpBaseDecl.hpp"
#include "Thyra_SolveSupportTypes.hpp"
#include "Teuchos_VerboseObject.hpp"
#include "Teuchos_PromotionTraits.hpp"

namespace Thyra {

template <class RangeScalar, class DomainScalar = RangeScalar>
class LinearOpWithSolveBase
  : virtual public LinearOpBase<RangeScalar,DomainScalar>
  , virtual public Teuchos::VerboseObject<LinearOpWithSolveBase<RangeScalar,DomainScalar> >
{
public:

  typedef typename Teuchos::PromotionTraits<RangeScalar,DomainScalar>::promote PromotedScalar;
  

  virtual void solve(
    const EConj                                   conj
    ,const MultiVectorBase<RangeScalar>           &B
    ,MultiVectorBase<DomainScalar>                *X
    ,const int                                    numBlocks             = 0
    ,const BlockSolveCriteria<PromotedScalar>     blockSolveCriteria[]  = NULL
    ,SolveStatus<PromotedScalar>                  blockSolveStatus[]    = NULL
    ) const = 0;



  virtual bool solveSupportsConj(EConj conj) const;

  virtual bool solveTransposeSupportsConj(EConj conj) const;

  virtual bool solveSupportsSolveMeasureType(EConj conj, const SolveMeasureType& solveMeasureType) const;

  virtual bool solveTransposeSupportsSolveMeasureType(EConj conj, const SolveMeasureType& solveMeasureType) const;

  virtual void solveTranspose(
    const EConj                                   conj
    ,const MultiVectorBase<DomainScalar>          &B
    ,MultiVectorBase<RangeScalar>                 *X
    ,const int                                    numBlocks             = 0
    ,const BlockSolveCriteria<PromotedScalar>     blockSolveCriteria[]  = NULL
    ,SolveStatus<PromotedScalar>                  blockSolveStatus[]    = NULL
    ) const;


};


template<class Scalar>
bool solveSupports(
  const LinearOpWithSolveBase<Scalar>   &A
  ,const EOpTransp                        A_trans
  )
{
  if( real_trans(A_trans) == NOTRANS ) {
    return A.solveSupportsConj(
      A_trans == NOTRANS ? NONCONJ_ELE : CONJ_ELE
      );
  }
  return A.solveTransposeSupportsConj(
    A_trans == TRANS ? NONCONJ_ELE : CONJ_ELE
    );
}


template<class Scalar>
bool solveSupportsSolveMeasureType(
  const LinearOpWithSolveBase<Scalar> &A,
  const EOpTransp A_trans,
  const SolveMeasureType &solveMeasureType
  )
{
  if( real_trans(A_trans) == NOTRANS ) {
    return A.solveSupportsSolveMeasureType(
      A_trans == NOTRANS ? NONCONJ_ELE : CONJ_ELE,
      solveMeasureType
      );
  }
  return A.solveTransposeSupportsSolveMeasureType(
    A_trans == TRANS ? NONCONJ_ELE : CONJ_ELE,
      solveMeasureType
    );
}


template<class Scalar>
void solve(
  const LinearOpWithSolveBase<Scalar> &M,
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &B,
  MultiVectorBase<Scalar> *X,
  const int numBlocks,
  const BlockSolveCriteria<Scalar> blockSolveCriteria[],
  SolveStatus<Scalar> blockSolveStatus[]
  )
{
  if (real_trans(M_trans) == NOTRANS)
    M.solve(
      transToConj(M_trans),B,X,numBlocks,blockSolveCriteria,blockSolveStatus);
  else
    M.solveTranspose(
      transToConj(M_trans),B,X,numBlocks,blockSolveCriteria,blockSolveStatus);
}


template<class Scalar>
SolveStatus<Scalar>
solve(
  const LinearOpWithSolveBase<Scalar>   &A
  ,const EOpTransp                        A_trans
  ,const MultiVectorBase<Scalar>        &B
  ,MultiVectorBase<Scalar>              *X
  ,const SolveCriteria<Scalar>          *solveCriteria = NULL
  )
{
  if(real_trans(A_trans)==NOTRANS)
    return solve(A,transToConj(A_trans),B,X,solveCriteria);
  return solveTranspose(A,transToConj(A_trans),B,X,solveCriteria);
}

template<class RangeScalar, class DomainScalar>
SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
solve(
  const LinearOpWithSolveBase<RangeScalar,DomainScalar>   &A
  ,const EConj                                            conj
  ,const MultiVectorBase<RangeScalar>                     &B
  ,MultiVectorBase<DomainScalar>                          *X
  ,const SolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          *solveCriteria  = NULL
  )
{
  typedef SolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>       SC;
  typedef BlockSolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>  BSC;
  typedef SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>         BSS;
  SC  defaultSolveCriteria;
  BSC blockSolveCriteria[1];
  BSS blockSolveStatus[1];
  blockSolveCriteria[0] = BSC(solveCriteria?*solveCriteria:defaultSolveCriteria,B.domain()->dim());
  A.solve(
    conj,B,X,1
    ,blockSolveCriteria
    ,blockSolveStatus
    );
  return blockSolveStatus[0];
}


template <class RangeScalar, class DomainScalar>
SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
solveTranspose(
  const LinearOpWithSolveBase<RangeScalar,DomainScalar>   &A
  ,const EConj                                            conj
  ,const MultiVectorBase<DomainScalar>                    &B
  ,MultiVectorBase<RangeScalar>                           *X
  ,const SolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          *solveCriteria  = NULL
  )
{
  typedef SolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>       SC;
  typedef BlockSolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>  BSC;
  typedef SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>         BSS;
  SC  defaultSolveCriteria;
  BSC blockSolveCriteria[1];
  BSS blockSolveStatus[1];
  blockSolveCriteria[0] = BSC(solveCriteria?*solveCriteria:defaultSolveCriteria,B.domain()->dim());
  A.solveTranspose(
    conj,B,X,1
    ,blockSolveCriteria
    ,blockSolveStatus
    );
  return blockSolveStatus[0];
}


template<class RangeScalar, class DomainScalar>
void solve(
  const LinearOpWithSolveBase<RangeScalar,DomainScalar>   &A
  ,const EConj                                            conj
  ,const MultiVectorBase<RangeScalar>                     &B
  ,MultiVectorBase<DomainScalar>                          *X
  ,const int                                              numBlocks
  ,const BlockSolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          blockSolveCriteria[] = NULL
  ,SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          blockSolveStatus[] = NULL
  )
{
  A.solve(conj,B,X,numBlocks,blockSolveCriteria,blockSolveStatus);
}


template <class RangeScalar, class DomainScalar>
void solveTranspose(
  const LinearOpWithSolveBase<RangeScalar,DomainScalar>   &A
  ,const EConj                                            conj
  ,const MultiVectorBase<DomainScalar>                    &B
  ,MultiVectorBase<RangeScalar>                           *X
  ,const int                                              numBlocks
  ,const BlockSolveCriteria<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          blockSolveCriteria[] = NULL
  ,SolveStatus<typename LinearOpWithSolveBase<RangeScalar,DomainScalar>::PromotedScalar>
                                                          blockSolveStatus[] = NULL
  )
{
  A.solveTranspose(conj,B,X,numBlocks,blockSolveCriteria,blockSolveStatus);
}


} // namespace Thyra

#endif // THYRA_LINEAR_OP_WITH_SOLVE_BASE_DECL_HPP

---