Anasazi  Version of the Day
AnasaziSaddleOperator.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
35 #ifndef ANASAZI_SADDLE_OPERATOR_HPP
36 #define ANASAZI_SADDLE_OPERATOR_HPP
37 
38 #include "AnasaziConfigDefs.hpp"
41 
42 #include "Teuchos_SerialDenseSolver.hpp"
43 
44 using Teuchos::RCP;
45 
46 enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
47 
48 namespace Anasazi {
49 namespace Experimental {
50 
51 template <class ScalarType, class MV, class OP>
52 class SaddleOperator : public TraceMinOp<ScalarType,SaddleContainer<ScalarType,MV>,OP>
53 {
55  typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
56 
57 public:
58  // Default constructor
59  SaddleOperator( ) { };
60  SaddleOperator( const Teuchos::RCP<OP> A, const Teuchos::RCP<const MV> B, PrecType pt=NO_PREC, const ScalarType alpha=1. );
61 
62  // Applies the saddle point operator to a "multivector"
63  void Apply(const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y) const;
64 
65  void removeIndices(const std::vector<int>& indicesToRemove) { A_->removeIndices(indicesToRemove); };
66 
67 private:
68  // A is the 1-1 block, and B the 1-2 block
69  Teuchos::RCP<OP> A_;
70  Teuchos::RCP<const MV> B_;
71  Teuchos::RCP<SerialDenseMatrix> Schur_;
72  PrecType pt_;
73  ScalarType alpha_;
74 };
75 
76 
77 
78 // Default constructor
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 )
81 {
82  // Get a pointer to A and B
83  A_ = A;
84  B_ = B;
85  pt_ = pt;
86  alpha_ = alpha;
87 
88  if(pt == BD_PREC)
89  {
90  // Form the Schur complement
91  int nvecs = MVT::GetNumberVecs(*B);
92  Teuchos::RCP<MV> AinvB = MVT::Clone(*B,nvecs);
93  Schur_ = rcp(new SerialDenseMatrix(nvecs,nvecs));
94 
95  A_->Apply(*B_,*AinvB);
96 
97  MVT::MvTransMv(1., *B_, *AinvB, *Schur_);
98  }
99 }
100 
101 // Applies the saddle point operator to a "multivector"
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
104 {
105  RCP<SerialDenseMatrix> Xlower = X.getLower();
106  RCP<SerialDenseMatrix> Ylower = Y.getLower();
107 
108  if(pt_ == NO_PREC)
109  {
110  // trans does literally nothing, because the operator is symmetric
111  // Y.bottom = B'X.top
112  MVT::MvTransMv(1., *B_, *(X.upper_), *Ylower);
113 
114  // Y.top = A*X.top+B*X.bottom
115  A_->Apply(*(X.upper_), *(Y.upper_));
116  MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
117  }
118  else if(pt_ == NONSYM)
119  {
120  // Y.bottom = -B'X.top
121  MVT::MvTransMv(-1., *B_, *(X.upper_), *Ylower);
122 
123  // Y.top = A*X.top+B*X.bottom
124  A_->Apply(*(X.upper_), *(Y.upper_));
125  MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
126  }
127  else if(pt_ == BD_PREC)
128  {
129  Teuchos::SerialDenseSolver<int,ScalarType> MySolver;
130 
131  // Solve A Y.X = X.X
132  A_->Apply(*(X.upper_),*(Y.upper_));
133 
134  // So, let me tell you a funny story about how the SerialDenseSolver destroys the original matrix...
135  Teuchos::RCP<SerialDenseMatrix> localSchur = Teuchos::rcp(new SerialDenseMatrix(*Schur_));
136 
137  // Solve the small system
138  MySolver.setMatrix(localSchur);
139  MySolver.setVectors(Ylower, Xlower);
140  MySolver.solve();
141  }
142  // Hermitian-Skew Hermitian splitting has some extra requirements
143  // We need B'B = I, which is true for standard eigenvalue problems, but not generalized
144  // We also need to use gmres, because our operator is no longer symmetric
145  else if(pt_ == HSS_PREC)
146  {
147 // std::cout << "applying preconditioner to";
148 // X.MvPrint(std::cout);
149 
150  // Solve (H + alpha I) Y1 = X
151  // 1. Apply preconditioner
152  A_->Apply(*(X.upper_),*(Y.upper_));
153  // 2. Scale by 1/alpha
154  *Ylower = *Xlower;
155  Ylower->scale(1./alpha_);
156 
157 // std::cout << "H preconditioning produced";
158 // Y.setLower(Ylower);
159 // Y.MvPrint(std::cout);
160 
161  // Solve (S + alpha I) Y = Y1
162  // 1. Y_lower = (B' Y1_upper + alpha Y1_lower) / (1 + alpha^2)
163  Teuchos::RCP<SerialDenseMatrix> Y1_lower = Teuchos::rcp(new SerialDenseMatrix(*Ylower));
164  MVT::MvTransMv(1,*B_,*(Y.upper_),*Ylower);
165 // std::cout << "Y'b1 " << *Ylower;
166  Y1_lower->scale(alpha_);
167 // std::cout << "alpha b2 " << *Y1_lower;
168  *Ylower += *Y1_lower;
169 // std::cout << "alpha b2 + Y'b1 " << *Ylower;
170  Ylower->scale(1/(1+alpha_*alpha_));
171  // 2. Y_upper = (Y1_upper - B Y_lower) / alpha
172  MVT::MvTimesMatAddMv(-1/alpha_,*B_,*Ylower,1/alpha_,*(Y.upper_));
173 
174 // std::cout << "preconditioning produced";
175 // Y.setLower(Ylower);
176 // Y.MvPrint(std::cout);
177  }
178  else
179  {
180  std::cout << "Not a valid preconditioner type\n";
181  }
182 
183  Y.setLower(Ylower);
184 
185 // std::cout << "result of applying operator";
186 // Y.MvPrint(std::cout);
187 }
188 
189 } // End namespace Experimental
190 
191 template<class ScalarType, class MV, class OP>
192 class OperatorTraits<ScalarType, Experimental::SaddleContainer<ScalarType,MV>, Experimental::SaddleOperator<ScalarType,MV,OP> >
193 {
194 public:
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); };
199 };
200 
201 } // end namespace Anasazi
202 
203 #ifdef HAVE_ANASAZI_BELOS
204 namespace Belos {
205 
206 template<class ScalarType, class MV, class OP>
207 class OperatorTraits<ScalarType, Anasazi::Experimental::SaddleContainer<ScalarType,MV>, Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP> >
208 {
209 public:
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); };
214 };
215 
216 } // end namespace Belos
217 #endif
218 
219 #endif
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...