48#ifndef ANASAZI_SADDLE_OPERATOR_HPP
49#define ANASAZI_SADDLE_OPERATOR_HPP
55#include "Teuchos_SerialDenseSolver.hpp"
59enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
64template <
class ScalarType,
class MV,
class OP>
65class SaddleOperator :
public TraceMinOp<ScalarType,SaddleContainer<ScalarType,MV>,OP>
68 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
72 SaddleOperator( ) { };
73 SaddleOperator(
const Teuchos::RCP<OP> A,
const Teuchos::RCP<const MV> B, PrecType pt=NO_PREC,
const ScalarType alpha=1. );
76 void Apply(
const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y)
const;
78 void removeIndices(
const std::vector<int>& indicesToRemove) { A_->removeIndices(indicesToRemove); };
83 Teuchos::RCP<const MV> B_;
84 Teuchos::RCP<SerialDenseMatrix> Schur_;
92template <
class ScalarType,
class MV,
class OP>
93SaddleOperator<ScalarType, MV, OP>::SaddleOperator(
const Teuchos::RCP<OP> A,
const Teuchos::RCP<const MV> B, PrecType pt,
const ScalarType alpha )
104 int nvecs = MVT::GetNumberVecs(*B);
105 Teuchos::RCP<MV> AinvB = MVT::Clone(*B,nvecs);
106 Schur_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
108 A_->Apply(*B_,*AinvB);
110 MVT::MvTransMv(1., *B_, *AinvB, *Schur_);
115template <
class ScalarType,
class MV,
class OP>
116void SaddleOperator<ScalarType, MV, OP>::Apply(
const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y)
const
118 RCP<SerialDenseMatrix> Xlower = X.getLower();
119 RCP<SerialDenseMatrix> Ylower = Y.getLower();
125 MVT::MvTransMv(1., *B_, *(X.upper_), *Ylower);
128 A_->Apply(*(X.upper_), *(Y.upper_));
129 MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
131 else if(pt_ == NONSYM)
134 MVT::MvTransMv(-1., *B_, *(X.upper_), *Ylower);
137 A_->Apply(*(X.upper_), *(Y.upper_));
138 MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
140 else if(pt_ == BD_PREC)
142 Teuchos::SerialDenseSolver<int,ScalarType> MySolver;
145 A_->Apply(*(X.upper_),*(Y.upper_));
148 Teuchos::RCP<SerialDenseMatrix> localSchur = Teuchos::rcp(
new SerialDenseMatrix(*Schur_));
151 MySolver.setMatrix(localSchur);
152 MySolver.setVectors(Ylower, Xlower);
158 else if(pt_ == HSS_PREC)
165 A_->Apply(*(X.upper_),*(Y.upper_));
168 Ylower->scale(1./alpha_);
176 Teuchos::RCP<SerialDenseMatrix> Y1_lower = Teuchos::rcp(
new SerialDenseMatrix(*Ylower));
177 MVT::MvTransMv(1,*B_,*(Y.upper_),*Ylower);
179 Y1_lower->scale(alpha_);
181 *Ylower += *Y1_lower;
183 Ylower->scale(1/(1+alpha_*alpha_));
185 MVT::MvTimesMatAddMv(-1/alpha_,*B_,*Ylower,1/alpha_,*(Y.upper_));
193 std::cout <<
"Not a valid preconditioner type\n";
204template<
class ScalarType,
class MV,
class OP>
205class OperatorTraits<ScalarType,
Experimental::SaddleContainer<ScalarType,MV>, Experimental::SaddleOperator<ScalarType,MV,OP> >
208 static void Apply(
const Experimental::SaddleOperator<ScalarType,MV,OP>& Op,
209 const Experimental::SaddleContainer<ScalarType,MV>& x,
210 Experimental::SaddleContainer<ScalarType,MV>& y)
211 { Op.Apply( x, y); };
216#ifdef HAVE_ANASAZI_BELOS
219template<
class ScalarType,
class MV,
class OP>
220class OperatorTraits<ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV>, Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP> >
223 static void Apply(
const Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP>& Op,
224 const Anasazi::Experimental::SaddleContainer<ScalarType,MV>& x,
225 Anasazi::Experimental::SaddleContainer<ScalarType,MV>& y)
226 { Op.Apply( x, y); };
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Stores a set of vectors of the form [V; L] where V is a distributed multivector and L is a serialdens...
Traits class which defines basic operations on multivectors.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.