17#ifndef ANASAZI_SADDLE_CONTAINER_HPP 
   18#define ANASAZI_SADDLE_CONTAINER_HPP 
   21#include "Teuchos_VerboseObject.hpp" 
   23#ifdef HAVE_ANASAZI_BELOS 
   24#include "BelosConfigDefs.hpp" 
   25#include "BelosMultiVecTraits.hpp" 
   34template <
class ScalarType, 
class MV>
 
   37  template <
class Scalar_, 
class MV_, 
class OP_> 
friend class SaddleOperator;
 
   41  typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
 
   42  const ScalarType ONE, ZERO;
 
   44  RCP<SerialDenseMatrix> lowerRaw_;
 
   45  std::vector<int> indices_;
 
   47  RCP<SerialDenseMatrix> getLower() 
const;
 
   48  void setLower(
const RCP<SerialDenseMatrix> lower);
 
   52  SaddleContainer( ) : ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero()) { };
 
   53  SaddleContainer( 
const RCP<MV> X, 
bool eye=
false );
 
   57  SaddleContainer<ScalarType, MV> * Clone(
const int nvecs) 
const;
 
   59  SaddleContainer<ScalarType, MV> * CloneCopy() 
const;
 
   61  SaddleContainer<ScalarType, MV> * CloneCopy(
const std::vector<int> &index) 
const;
 
   62  SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const std::vector<int> &index) 
const;
 
   63  SaddleContainer<ScalarType, MV> * CloneViewNonConst(
const Teuchos::Range1D& index) 
const;
 
   64  const SaddleContainer<ScalarType, MV> * CloneView(
const std::vector<int> &index) 
const;
 
   65  const SaddleContainer<ScalarType, MV> * CloneView(
const Teuchos::Range1D& index) 
const;
 
   66  ptrdiff_t GetGlobalLength()
 const { 
return MVT::GetGlobalLength(*upper_) + lowerRaw_->numRows(); };
 
   69  void MvTimesMatAddMv(ScalarType alpha, 
const SaddleContainer<ScalarType,MV> &A,
 
   70                       const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
 
   73  void MvAddMv(ScalarType alpha, 
const SaddleContainer<ScalarType,MV>& A,
 
   74               ScalarType beta,  
const SaddleContainer<ScalarType,MV>& B);
 
   76  void MvScale( ScalarType alpha );
 
   78  void MvScale( 
const std::vector<ScalarType>& alpha );
 
   80  void MvTransMv (ScalarType alpha, 
const SaddleContainer<ScalarType, MV>& A,
 
   81                  Teuchos::SerialDenseMatrix< int, ScalarType >& B) 
const;
 
   83  void MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b) 
const;
 
   85  void MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec) 
const;
 
   89  void SetBlock (
const SaddleContainer<ScalarType, MV>& A, 
const std::vector<int> &index);
 
   91  void Assign (
const SaddleContainer<ScalarType, MV>&A);
 
   95  void  MvInit (ScalarType alpha);
 
   97  void MvPrint (std::ostream &os) 
const;
 
  103template <
class ScalarType, 
class MV>
 
  104RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > SaddleContainer<ScalarType, MV>::getLower()
 const 
  111  int nrows = lowerRaw_->numRows();
 
  112  int ncols = indices_.size();
 
  114  RCP<SerialDenseMatrix> lower = rcp(
new SerialDenseMatrix(nrows,ncols,
false));
 
  116  for(
int r=0; r<nrows; r++)
 
  118    for(
int c=0; c<ncols; c++)
 
  120      (*lower)(r,c) = (*lowerRaw_)(r,indices_[c]);
 
  130template <
class ScalarType, 
class MV>
 
  131void SaddleContainer<ScalarType, MV>::setLower(
const RCP<SerialDenseMatrix> lower)
 
  139  int nrows = lowerRaw_->numRows();
 
  140  int ncols = indices_.size();
 
  142  for(
int r=0; r<nrows; r++)
 
  144    for(
int c=0; c<ncols; c++)
 
  146      (*lowerRaw_)(r,indices_[c]) = (*lower)(r,c);
 
  154template <
class ScalarType, 
class MV>
 
  155SaddleContainer<ScalarType, MV>::SaddleContainer( 
const RCP<MV> X, 
bool eye )
 
  156: ONE(Teuchos::ScalarTraits<ScalarType>::one()), ZERO(Teuchos::ScalarTraits<ScalarType>::zero())
 
  158  int nvecs = MVT::GetNumberVecs(*X);
 
  163    upper_ = MVT::Clone(*X, nvecs);
 
  164    MVT::MvInit(*upper_);
 
  167    lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
 
  168    for(
int i=0; i < nvecs; i++)
 
  169      (*lowerRaw_)(i,i) = ONE;
 
  177    lowerRaw_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
 
  184template <
class ScalarType, 
class MV>
 
  185SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::Clone(
const int nvecs)
 const 
  187  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  189  newSC->upper_ = MVT::Clone(*upper_,nvecs);
 
  190  newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(lowerRaw_->numRows(),nvecs));
 
  198template <
class ScalarType, 
class MV>
 
  199SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy()
 const 
  201  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  203  newSC->upper_ = MVT::CloneCopy(*upper_);
 
  204  newSC->lowerRaw_ = getLower();
 
  212template <
class ScalarType, 
class MV>
 
  213SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneCopy(
const std::vector< int > &index)
 const 
  215  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  217  newSC->upper_ = MVT::CloneCopy(*upper_,index);
 
  219  int ncols = index.size();
 
  220  int nrows = lowerRaw_->numRows();
 
  221  RCP<SerialDenseMatrix> lower = getLower();
 
  222  newSC->lowerRaw_ = rcp(
new SerialDenseMatrix(nrows,ncols));
 
  223  for(
int c=0; c < ncols; c++)
 
  225    for(
int r=0; r < nrows; r++)
 
  226      (*newSC->lowerRaw_)(r,c) = (*lower)(r,index[c]);
 
  235template <
class ScalarType, 
class MV>
 
  236SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const std::vector<int> &index)
 const 
  238  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  240  newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
 
  242  newSC->lowerRaw_ = lowerRaw_;
 
  244  if(!indices_.empty())
 
  246    newSC->indices_.resize(index.size());
 
  247    for(
unsigned int i=0; i<index.size(); i++)
 
  249      newSC->indices_[i] = indices_[index[i]];
 
  254    newSC->indices_ = index;
 
  261template <
class ScalarType, 
class MV>
 
  262SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneViewNonConst(
const Teuchos::Range1D& index)
 const 
  264  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  266  newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
 
  268  newSC->lowerRaw_ = lowerRaw_;
 
  270  newSC->indices_.resize(index.size());
 
  271  for(
unsigned int i=0; i<index.size(); i++)
 
  273    newSC->indices_[i] = indices_[index.lbound()+i];
 
  282template <
class ScalarType, 
class MV>
 
  283const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const std::vector<int> &index)
 const 
  285  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  287  newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
 
  289  newSC->lowerRaw_ = lowerRaw_;
 
  291  if(!indices_.empty())
 
  293    newSC->indices_.resize(index.size());
 
  294    for(
unsigned int i=0; i<index.size(); i++)
 
  296      newSC->indices_[i] = indices_[index[i]];
 
  301    newSC->indices_ = index;
 
  308template <
class ScalarType, 
class MV>
 
  309const SaddleContainer<ScalarType, MV> * SaddleContainer<ScalarType, MV>::CloneView(
const Teuchos::Range1D& index)
 const 
  311  SaddleContainer<ScalarType, MV> * newSC = 
new SaddleContainer<ScalarType, MV>();
 
  313  newSC->upper_ = MVT::CloneViewNonConst(*upper_,index);
 
  315  newSC->lowerRaw_ = lowerRaw_;
 
  317  newSC->indices_.resize(index.size());
 
  318  for(
unsigned int i=0; i<index.size(); i++)
 
  320    newSC->indices_[i] = indices_[index.lbound()+i];
 
  330template <
class ScalarType, 
class MV>
 
  331void SaddleContainer<ScalarType, MV>::MvTimesMatAddMv(ScalarType alpha, 
const SaddleContainer<ScalarType,MV> &A,
 
  332                                                      const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
 
  335  MVT::MvTimesMatAddMv(alpha,*(A.upper_),B,beta,*upper_);
 
  336  RCP<SerialDenseMatrix> lower = getLower();
 
  337  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  338  lower->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,alpha,*Alower,B,beta);
 
  345template <
class ScalarType, 
class MV>
 
  346void SaddleContainer<ScalarType, MV>::MvAddMv(ScalarType alpha, 
const SaddleContainer<ScalarType,MV>& A,
 
  347                                              ScalarType beta,  
const SaddleContainer<ScalarType,MV>& B)
 
  349  MVT::MvAddMv(alpha, *(A.upper_), beta, *(B.upper_), *upper_);
 
  351  RCP<SerialDenseMatrix> lower = getLower();
 
  352  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  353  RCP<SerialDenseMatrix> Blower = B.getLower();
 
  360  lower->assign(*Alower);
 
  366  else if(beta == -ONE)
 
  368  else if(beta != ZERO)
 
  370    SerialDenseMatrix scaledB(*Blower);
 
  381template <
class ScalarType, 
class MV>
 
  382void SaddleContainer<ScalarType, MV>::MvScale( ScalarType alpha )
 
  384  MVT::MvScale(*upper_, alpha);
 
  387  RCP<SerialDenseMatrix> lower = getLower();
 
  395template <
class ScalarType, 
class MV>
 
  396void SaddleContainer<ScalarType, MV>::MvScale( 
const std::vector<ScalarType>& alpha )
 
  398  MVT::MvScale(*upper_, alpha);
 
  400  RCP<SerialDenseMatrix> lower = getLower();
 
  402  int nrows = lower->numRows();
 
  403  int ncols = lower->numCols();
 
  405  for(
int c=0; c<ncols; c++)
 
  407    for(
int r=0; r<nrows; r++)
 
  408      (*lower)(r,c) *= alpha[c];
 
  418template <
class ScalarType, 
class MV>
 
  419void SaddleContainer<ScalarType, MV>::MvTransMv (ScalarType alpha, 
const SaddleContainer<ScalarType, MV>& A,
 
  420                                                 Teuchos::SerialDenseMatrix< int, ScalarType >& B)
 const 
  422  MVT::MvTransMv(alpha,*(A.upper_),*upper_,B);
 
  423  RCP<SerialDenseMatrix> lower = getLower();
 
  424  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  425  B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,alpha,*(Alower),*lower,ONE);
 
  431template <
class ScalarType, 
class MV>
 
  432void SaddleContainer<ScalarType, MV>::MvDot (
const SaddleContainer<ScalarType, MV>& A, std::vector<ScalarType> &b)
 const 
  434  MVT::MvDot(*upper_, *(A.upper_), b);
 
  436  RCP<SerialDenseMatrix> lower = getLower();
 
  437  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  439  int nrows = lower->numRows();
 
  440  int ncols = lower->numCols();
 
  442  for(
int c=0; c < ncols; c++)
 
  444    for(
int r=0; r < nrows; r++)
 
  446      b[c] += ((*Alower)(r,c) * (*lower)(r,c));
 
  455template <
class ScalarType, 
class MV>
 
  456void SaddleContainer<ScalarType, MV>::MvNorm ( std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec)
 const 
  459  MvDot(*
this,normvec);
 
  460  for(
unsigned int i=0; i<normvec.size(); i++)
 
  461    normvec[i] = sqrt(normvec[i]);
 
  469template <
class ScalarType, 
class MV>
 
  470void SaddleContainer<ScalarType, MV>::SetBlock (
const SaddleContainer<ScalarType, MV>& A, 
const std::vector<int> &index)
 
  472  MVT::SetBlock(*(A.upper_), index, *upper_);
 
  474  RCP<SerialDenseMatrix> lower = getLower();
 
  475  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  477  int nrows = lower->numRows();
 
  479  int nvecs = index.size();
 
  480  for(
int c=0; c<nvecs; c++)
 
  482    for(
int r=0; r<nrows; r++)
 
  483      (*lower)(r,index[c]) = (*Alower)(r,c);
 
  492template <
class ScalarType, 
class MV>
 
  493void SaddleContainer<ScalarType, MV>::Assign (
const SaddleContainer<ScalarType, MV>&A)
 
  495  MVT::Assign(*(A.upper_),*(upper_));
 
  497  RCP<SerialDenseMatrix> lower = getLower();
 
  498  RCP<SerialDenseMatrix> Alower = A.getLower();
 
  509template <
class ScalarType, 
class MV>
 
  510void  SaddleContainer<ScalarType, MV>::MvRandom ()
 
  512  MVT::MvRandom(*upper_);
 
  514  RCP<SerialDenseMatrix> lower = getLower();
 
  522template <
class ScalarType, 
class MV>
 
  523void SaddleContainer<ScalarType, MV>::MvInit (ScalarType alpha)
 
  525  MVT::MvInit(*upper_,alpha);
 
  527  RCP<SerialDenseMatrix> lower = getLower();
 
  528  lower->putScalar(alpha);
 
  535template <
class ScalarType, 
class MV>
 
  536void SaddleContainer<ScalarType, MV>::MvPrint (std::ostream &os)
 const 
  538  RCP<SerialDenseMatrix> lower = getLower();
 
  542  os << 
"Object SaddleContainer" << std::endl;
 
  544  upper_->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
 
  547  os << 
"Y\n" << *lower << std::endl;
 
  552template<
class ScalarType, 
class MV >
 
  553class MultiVecTraits<ScalarType,
Experimental::SaddleContainer<ScalarType, MV> >
 
  555typedef Experimental::SaddleContainer<ScalarType,MV>  SC;
 
  558  static RCP<SC > Clone( 
const SC& mv, 
const int numvecs )
 
  559    { 
return rcp( 
const_cast<SC&
>(mv).Clone(numvecs) ); }
 
  561  static RCP<SC > CloneCopy( 
const SC& mv )
 
  562    { 
return rcp( 
const_cast<SC&
>(mv).CloneCopy() ); }
 
  564  static RCP<SC > CloneCopy( 
const SC& mv, 
const std::vector<int>& index )
 
  565    { 
return rcp( 
const_cast<SC&
>(mv).CloneCopy(index) ); }
 
  567  static ptrdiff_t GetGlobalLength( 
const SC& mv )
 
  568    { 
return mv.GetGlobalLength(); }
 
  570  static int GetNumberVecs( 
const SC& mv )
 
  571    { 
return mv.GetNumberVecs(); }
 
  573  static void MvTimesMatAddMv( ScalarType alpha, 
const SC& A,
 
  574                               const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
 
  575                               ScalarType beta, SC& mv )
 
  576    { mv.MvTimesMatAddMv(alpha, A, B, beta); }
 
  578  static void MvAddMv( ScalarType alpha, 
const SC& A, ScalarType beta, 
const SC& B, SC& mv )
 
  579    { mv.MvAddMv(alpha, A, beta, B); }
 
  581  static void MvTransMv( ScalarType alpha, 
const SC& A, 
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
 
  582    { mv.MvTransMv(alpha, A, B); }
 
  584  static void MvDot( 
const SC& mv, 
const SC& A, std::vector<ScalarType> & b)
 
  587  static void MvScale ( SC& mv, ScalarType alpha )
 
  588    { mv.MvScale( alpha ); }
 
  590  static void MvScale ( SC& mv, 
const std::vector<ScalarType>& alpha )
 
  591    { mv.MvScale( alpha ); }
 
  593  static void MvNorm( 
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec)
 
  594    { mv.MvNorm(normvec); }
 
  596  static void SetBlock( 
const SC& A, 
const std::vector<int>& index, SC& mv )
 
  597    { mv.SetBlock(A, index); }
 
  599  static void Assign( 
const SC& A, SC& mv )
 
  602  static void MvRandom( SC& mv )
 
  605  static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
 
  606    { mv.MvInit(alpha); }
 
  608  static void MvPrint( 
const SC& mv, std::ostream& os )
 
  614#ifdef HAVE_ANASAZI_BELOS 
  618template<
class ScalarType, 
class MV >
 
  619class MultiVecTraits< ScalarType, 
Anasazi::Experimental::SaddleContainer<ScalarType,MV> >
 
  621typedef Anasazi::Experimental::SaddleContainer<ScalarType,MV>  SC;
 
  623  static RCP<SC > Clone( 
const SC& mv, 
const int numvecs )
 
  624    { 
return rcp( 
const_cast<SC&
>(mv).Clone(numvecs) ); }
 
  626  static RCP<SC > CloneCopy( 
const SC& mv )
 
  627    { 
return rcp( 
const_cast<SC&
>(mv).CloneCopy() ); }
 
  629  static RCP<SC > CloneCopy( 
const SC& mv, 
const std::vector<int>& index )
 
  630    { 
return rcp( 
const_cast<SC&
>(mv).CloneCopy(index) ); }
 
  632  static RCP<SC> CloneViewNonConst( SC& mv, 
const std::vector<int>& index )
 
  633    { 
return rcp( mv.CloneViewNonConst(index) ); }
 
  635  static RCP<SC> CloneViewNonConst( SC& mv, 
const Teuchos::Range1D& index )
 
  636    { 
return rcp( mv.CloneViewNonConst(index) ); }
 
  638  static RCP<const SC> CloneView( 
const SC& mv, 
const std::vector<int> & index )
 
  639    { 
return rcp( mv.CloneView(index) ); }
 
  641  static RCP<const SC> CloneView( 
const SC& mv, 
const Teuchos::Range1D& index )
 
  642    { 
return rcp( mv.CloneView(index) ); }
 
  644  static ptrdiff_t GetGlobalLength( 
const SC& mv )
 
  645    { 
return mv.GetGlobalLength(); }
 
  647  static int GetNumberVecs( 
const SC& mv )
 
  648    { 
return mv.GetNumberVecs(); }
 
  650  static void MvTimesMatAddMv( ScalarType alpha, 
const SC& A,
 
  651                               const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
 
  652                               ScalarType beta, SC& mv )
 
  653    { mv.MvTimesMatAddMv(alpha, A, B, beta); }
 
  655  static void MvAddMv( ScalarType alpha, 
const SC& A, ScalarType beta, 
const SC& B, SC& mv )
 
  656    { mv.MvAddMv(alpha, A, beta, B); }
 
  658  static void MvTransMv( ScalarType alpha, 
const SC& A, 
const SC& mv, Teuchos::SerialDenseMatrix<int,ScalarType>& B)
 
  659    { mv.MvTransMv(alpha, A, B); }
 
  661  static void MvDot( 
const SC& mv, 
const SC& A, std::vector<ScalarType> & b)
 
  664  static void MvScale ( SC& mv, ScalarType alpha )
 
  665    { mv.MvScale( alpha ); }
 
  667  static void MvScale ( SC& mv, 
const std::vector<ScalarType>& alpha )
 
  668    { mv.MvScale( alpha ); }
 
  671  static void MvNorm( 
const SC& mv, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec, NormType type=TwoNorm)
 
  672    { mv.MvNorm(normvec); }
 
  674  static void SetBlock( 
const SC& A, 
const std::vector<int>& index, SC& mv )
 
  675    { mv.SetBlock(A, index); }
 
  677  static void Assign( 
const SC& A, SC& mv )
 
  680  static void MvRandom( SC& mv )
 
  683  static void MvInit( SC& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
 
  684    { mv.MvInit(alpha); }
 
  686  static void MvPrint( 
const SC& mv, std::ostream& os )
 
  689#ifdef HAVE_BELOS_TSQR 
  702  typedef Belos::details::StubTsqrAdapter<SC> tsqr_adaptor_type;
 
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
 
Traits class which defines basic operations on multivectors.
 
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
 
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
 
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.