48 #ifndef ANASAZI_SADDLE_OPERATOR_HPP 49 #define ANASAZI_SADDLE_OPERATOR_HPP 55 #include "Teuchos_SerialDenseSolver.hpp" 59 enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
64 template <
class ScalarType,
class MV,
class OP>
65 class 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_;
92 template <
class ScalarType,
class MV,
class OP>
93 SaddleOperator<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_);
115 template <
class ScalarType,
class MV,
class OP>
116 void 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";
204 template<
class ScalarType,
class MV,
class OP>
205 class 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 219 template<
class ScalarType,
class MV,
class OP>
220 class 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); };
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...