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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
48 #ifndef ANASAZI_SADDLE_OPERATOR_HPP
49 #define ANASAZI_SADDLE_OPERATOR_HPP
50 
51 #include "AnasaziConfigDefs.hpp"
54 
55 #include "Teuchos_SerialDenseSolver.hpp"
56 
57 using Teuchos::RCP;
58 
59 enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
60 
61 namespace Anasazi {
62 namespace Experimental {
63 
64 template <class ScalarType, class MV, class OP>
65 class SaddleOperator : public TraceMinOp<ScalarType,SaddleContainer<ScalarType,MV>,OP>
66 {
68  typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
69 
70 public:
71  // Default constructor
72  SaddleOperator( ) { };
73  SaddleOperator( const Teuchos::RCP<OP> A, const Teuchos::RCP<const MV> B, PrecType pt=NO_PREC, const ScalarType alpha=1. );
74 
75  // Applies the saddle point operator to a "multivector"
76  void Apply(const SaddleContainer<ScalarType,MV>& X, SaddleContainer<ScalarType,MV>& Y) const;
77 
78  void removeIndices(const std::vector<int>& indicesToRemove) { A_->removeIndices(indicesToRemove); };
79 
80 private:
81  // A is the 1-1 block, and B the 1-2 block
82  Teuchos::RCP<OP> A_;
83  Teuchos::RCP<const MV> B_;
84  Teuchos::RCP<SerialDenseMatrix> Schur_;
85  PrecType pt_;
86  ScalarType alpha_;
87 };
88 
89 
90 
91 // Default constructor
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 )
94 {
95  // Get a pointer to A and B
96  A_ = A;
97  B_ = B;
98  pt_ = pt;
99  alpha_ = alpha;
100 
101  if(pt == BD_PREC)
102  {
103  // Form the Schur complement
104  int nvecs = MVT::GetNumberVecs(*B);
105  Teuchos::RCP<MV> AinvB = MVT::Clone(*B,nvecs);
106  Schur_ = rcp(new SerialDenseMatrix(nvecs,nvecs));
107 
108  A_->Apply(*B_,*AinvB);
109 
110  MVT::MvTransMv(1., *B_, *AinvB, *Schur_);
111  }
112 }
113 
114 // Applies the saddle point operator to a "multivector"
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
117 {
118  RCP<SerialDenseMatrix> Xlower = X.getLower();
119  RCP<SerialDenseMatrix> Ylower = Y.getLower();
120 
121  if(pt_ == NO_PREC)
122  {
123  // trans does literally nothing, because the operator is symmetric
124  // Y.bottom = B'X.top
125  MVT::MvTransMv(1., *B_, *(X.upper_), *Ylower);
126 
127  // Y.top = A*X.top+B*X.bottom
128  A_->Apply(*(X.upper_), *(Y.upper_));
129  MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
130  }
131  else if(pt_ == NONSYM)
132  {
133  // Y.bottom = -B'X.top
134  MVT::MvTransMv(-1., *B_, *(X.upper_), *Ylower);
135 
136  // Y.top = A*X.top+B*X.bottom
137  A_->Apply(*(X.upper_), *(Y.upper_));
138  MVT::MvTimesMatAddMv(1., *B_, *Xlower, 1., *(Y.upper_));
139  }
140  else if(pt_ == BD_PREC)
141  {
142  Teuchos::SerialDenseSolver<int,ScalarType> MySolver;
143 
144  // Solve A Y.X = X.X
145  A_->Apply(*(X.upper_),*(Y.upper_));
146 
147  // So, let me tell you a funny story about how the SerialDenseSolver destroys the original matrix...
148  Teuchos::RCP<SerialDenseMatrix> localSchur = Teuchos::rcp(new SerialDenseMatrix(*Schur_));
149 
150  // Solve the small system
151  MySolver.setMatrix(localSchur);
152  MySolver.setVectors(Ylower, Xlower);
153  MySolver.solve();
154  }
155  // Hermitian-Skew Hermitian splitting has some extra requirements
156  // We need B'B = I, which is true for standard eigenvalue problems, but not generalized
157  // We also need to use gmres, because our operator is no longer symmetric
158  else if(pt_ == HSS_PREC)
159  {
160 // std::cout << "applying preconditioner to";
161 // X.MvPrint(std::cout);
162 
163  // Solve (H + alpha I) Y1 = X
164  // 1. Apply preconditioner
165  A_->Apply(*(X.upper_),*(Y.upper_));
166  // 2. Scale by 1/alpha
167  *Ylower = *Xlower;
168  Ylower->scale(1./alpha_);
169 
170 // std::cout << "H preconditioning produced";
171 // Y.setLower(Ylower);
172 // Y.MvPrint(std::cout);
173 
174  // Solve (S + alpha I) Y = Y1
175  // 1. Y_lower = (B' Y1_upper + alpha Y1_lower) / (1 + alpha^2)
176  Teuchos::RCP<SerialDenseMatrix> Y1_lower = Teuchos::rcp(new SerialDenseMatrix(*Ylower));
177  MVT::MvTransMv(1,*B_,*(Y.upper_),*Ylower);
178 // std::cout << "Y'b1 " << *Ylower;
179  Y1_lower->scale(alpha_);
180 // std::cout << "alpha b2 " << *Y1_lower;
181  *Ylower += *Y1_lower;
182 // std::cout << "alpha b2 + Y'b1 " << *Ylower;
183  Ylower->scale(1/(1+alpha_*alpha_));
184  // 2. Y_upper = (Y1_upper - B Y_lower) / alpha
185  MVT::MvTimesMatAddMv(-1/alpha_,*B_,*Ylower,1/alpha_,*(Y.upper_));
186 
187 // std::cout << "preconditioning produced";
188 // Y.setLower(Ylower);
189 // Y.MvPrint(std::cout);
190  }
191  else
192  {
193  std::cout << "Not a valid preconditioner type\n";
194  }
195 
196  Y.setLower(Ylower);
197 
198 // std::cout << "result of applying operator";
199 // Y.MvPrint(std::cout);
200 }
201 
202 } // End namespace Experimental
203 
204 template<class ScalarType, class MV, class OP>
205 class OperatorTraits<ScalarType, Experimental::SaddleContainer<ScalarType,MV>, Experimental::SaddleOperator<ScalarType,MV,OP> >
206 {
207 public:
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); };
212 };
213 
214 } // end namespace Anasazi
215 
216 #ifdef HAVE_ANASAZI_BELOS
217 namespace Belos {
218 
219 template<class ScalarType, class MV, class OP>
220 class OperatorTraits<ScalarType, Anasazi::Experimental::SaddleContainer<ScalarType,MV>, Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP> >
221 {
222 public:
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); };
227 };
228 
229 } // end namespace Belos
230 #endif
231 
232 #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...