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 )
93 Schur_ = rcp(
new SerialDenseMatrix(nvecs,nvecs));
95 A_->Apply(*B_,*AinvB);
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();
115 A_->Apply(*(X.upper_), *(Y.upper_));
118 else if(pt_ == NONSYM)
124 A_->Apply(*(X.upper_), *(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));
166 Y1_lower->scale(alpha_);
168 *Ylower += *Y1_lower;
170 Ylower->scale(1/(1+alpha_*alpha_));
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 MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Stores a set of vectors of the form [V; L] where V is a distributed multivector and L is a serialdens...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.