Anasazi Version of the Day
Loading...
Searching...
No Matches
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
57using Teuchos::RCP;
58
59enum PrecType {NO_PREC, NONSYM, BD_PREC, HSS_PREC};
60
61namespace Anasazi {
62namespace Experimental {
63
64template <class ScalarType, class MV, class OP>
65class SaddleOperator : public TraceMinOp<ScalarType,SaddleContainer<ScalarType,MV>,OP>
66{
68 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SerialDenseMatrix;
69
70public:
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
80private:
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
92template <class ScalarType, class MV, class OP>
93SaddleOperator<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"
115template <class ScalarType, class MV, class OP>
116void 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
204template<class ScalarType, class MV, class OP>
205class OperatorTraits<ScalarType, Experimental::SaddleContainer<ScalarType,MV>, Experimental::SaddleOperator<ScalarType,MV,OP> >
206{
207public:
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
217namespace Belos {
218
219template<class ScalarType, class MV, class OP>
220class OperatorTraits<ScalarType, Anasazi::Experimental::SaddleContainer<ScalarType,MV>, Anasazi::Experimental::SaddleOperator<ScalarType,MV,OP> >
221{
222public:
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
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...
Traits class which defines basic operations on multivectors.
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.