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 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.