35 #ifndef ANASAZI_SADDLE_OPERATOR_HPP 36 #define ANASAZI_SADDLE_OPERATOR_HPP 42 #include "Teuchos_SerialDenseSolver.hpp" 46 enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
51 template <
class ScalarType,
class MV,
class OP>
52 class SaddleOperator :
public TraceMinOp<ScalarType,SaddleContainer<ScalarType,MV>,OP>
55 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
59 SaddleOperator( ) { };
60 SaddleOperator(
const Teuchos::RCP<OP> A,
const Teuchos::RCP<const MV> B, PrecType pt=NO_PREC,
const ScalarType alpha=1. );
63 void Apply(
const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y)
const;
65 void removeIndices(
const std::vector<int>& indicesToRemove) { A_->removeIndices(indicesToRemove); };
70 Teuchos::RCP<const MV> B_;
71 Teuchos::RCP<SerialDenseMatrix> Schur_;
79 template <
class ScalarType,
class MV,
class OP>
80 SaddleOperator<ScalarType, MV, OP>::SaddleOperator(
const Teuchos::RCP<OP> A,
const Teuchos::RCP<const MV> B, PrecType pt,
const ScalarType alpha )
91 int nvecs = MVT::GetNumberVecs(*B);
92 Teuchos::RCP<MV> AinvB = MVT::Clone(*B,nvecs);
93 Schur_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
95 A_->Apply(*B_,*AinvB);
97 MVT::MvTransMv(1., *B_, *AinvB, *Schur_);
102 template <
class ScalarType,
class MV,
class OP>
103 void SaddleOperator<ScalarType, MV, OP>::Apply(
const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y)
const 105 RCP<SerialDenseMatrix> Xlower = X.getLower();
106 RCP<SerialDenseMatrix> Ylower = Y.getLower();
112 MVT::MvTransMv(1., *B_, *(X.upper_), *Ylower);
115 A_->Apply(*(X.upper_), *(Y.upper_));
116 MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
118 else if(pt_ == NONSYM)
121 MVT::MvTransMv(-1., *B_, *(X.upper_), *Ylower);
124 A_->Apply(*(X.upper_), *(Y.upper_));
125 MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
127 else if(pt_ == BD_PREC)
129 Teuchos::SerialDenseSolver<int,ScalarType> MySolver;
132 A_->Apply(*(X.upper_),*(Y.upper_));
135 Teuchos::RCP<SerialDenseMatrix> localSchur = Teuchos::rcp(
new SerialDenseMatrix(*Schur_));
138 MySolver.setMatrix(localSchur);
139 MySolver.setVectors(Ylower, Xlower);
145 else if(pt_ == HSS_PREC)
152 A_->Apply(*(X.upper_),*(Y.upper_));
155 Ylower->scale(1./alpha_);
163 Teuchos::RCP<SerialDenseMatrix> Y1_lower = Teuchos::rcp(
new SerialDenseMatrix(*Ylower));
164 MVT::MvTransMv(1,*B_,*(Y.upper_),*Ylower);
166 Y1_lower->scale(alpha_);
168 *Ylower += *Y1_lower;
170 Ylower->scale(1/(1+alpha_*alpha_));
172 MVT::MvTimesMatAddMv(-1/alpha_,*B_,*Ylower,1/alpha_,*(Y.upper_));
180 std::cout <<
"Not a valid preconditioner type\n";
191 template<
class ScalarType,
class MV,
class OP>
192 class OperatorTraits<ScalarType,
Experimental::SaddleContainer<ScalarType,MV>, Experimental::SaddleOperator<ScalarType,MV,OP> >
195 static void Apply(
const Experimental::SaddleOperator<ScalarType,MV,OP>& Op,
196 const Experimental::SaddleContainer<ScalarType,MV>& x,
197 Experimental::SaddleContainer<ScalarType,MV>& y)
198 { Op.Apply( x, y); };
203 #ifdef HAVE_ANASAZI_BELOS 206 template<
class ScalarType,
class MV,
class OP>
207 class OperatorTraits<ScalarType,
Anasazi::Experimental::SaddleContainer<ScalarType,MV>, Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP> >
210 static void Apply(
const Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP>& Op,
211 const Anasazi::Experimental::SaddleContainer<ScalarType,MV>& x,
212 Anasazi::Experimental::SaddleContainer<ScalarType,MV>& y)
213 { Op.Apply( x, y); };
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...
Traits class which defines basic operations on multivectors.
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...