Teko Version of the Day
Loading...
Searching...
No Matches
Teko_NeumannSeriesPreconditionerFactory.hpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#ifndef __Teko_NeumannSeriesPreconditionerFactory_hpp__
48#define __Teko_NeumannSeriesPreconditionerFactory_hpp__
49
50#include "Teko_NeumannSeriesPreconditionerFactoryDecl.hpp"
51
52#include "Thyra_DefaultPreconditioner.hpp"
53#include "Thyra_DefaultPreconditioner.hpp"
54#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
55#include "Thyra_DefaultAddedLinearOp.hpp"
56#include "Thyra_DefaultMultipliedLinearOp.hpp"
57#include "Thyra_DefaultIdentityLinearOp.hpp"
58
59#include "Teuchos_Array.hpp"
60#include "Teuchos_StandardParameterEntryValidators.hpp"
61
62namespace Teko {
63
64using Teuchos::RCP;
65using Teuchos::rcp;
66
67static RCP<Teuchos::StringToIntegralParameterEntryValidator<Teko::DiagonalType> > scalingTypeVdtor;
68
69template <typename ScalarT>
70NeumannSeriesPreconditionerFactory<ScalarT>::NeumannSeriesPreconditionerFactory()
71 : numberOfTerms_(1), scalingType_(Teko::NotDiag)
72{
73}
74
76template <typename ScalarT>
77bool NeumannSeriesPreconditionerFactory<ScalarT>::isCompatible(const Thyra::LinearOpSourceBase<ScalarT> &/* fwdOpSrc */) const
78{
79 return true;
80}
81
83template <typename ScalarT>
84RCP<Thyra::PreconditionerBase<ScalarT> > NeumannSeriesPreconditionerFactory<ScalarT>::createPrec() const
85{
86 return rcp(new Thyra::DefaultPreconditioner<ScalarT>());
87}
88
97template <typename ScalarT>
98void NeumannSeriesPreconditionerFactory<ScalarT>::initializePrec(const RCP<const Thyra::LinearOpSourceBase<ScalarT> > & fwdOpSrc,
99 Thyra::PreconditionerBase<ScalarT> * prec,
100 const Thyra::ESupportSolveUse /* supportSolveUse */) const
101{
102 using Thyra::scale;
103 using Thyra::add;
104 using Thyra::multiply;
105
106 RCP<const Thyra::LinearOpBase<ScalarT> > M; // left-preconditioner
107 RCP<const Thyra::LinearOpBase<ScalarT> > A = fwdOpSrc->getOp();
108 if(scalingType_!=Teko::NotDiag) {
109 M = Teko::getInvDiagonalOp(A,scalingType_);
110 A = Thyra::multiply(M,A);
111 }
112
113 RCP<const Thyra::LinearOpBase<ScalarT> > id = Thyra::identity<ScalarT>(A->range()); // I
114 RCP<const Thyra::LinearOpBase<ScalarT> > idMA = add(id,scale(-1.0,A)); // I - A
115
116
117 RCP<const Thyra::LinearOpBase<ScalarT> > precOp;
118 if(numberOfTerms_==1) {
119 // no terms requested, just return identity matrix
120 precOp = id;
121 }
122 else {
123 int iters = numberOfTerms_-1;
124 // use Horner's method to computed higher order polynomials
125 precOp = add(scale(2.0,id),scale(-1.0,A)); // I + (I - A)
126 for(int i=0;i<iters;i++)
127 precOp = add(id,multiply(idMA,precOp)); // I + (I - A) * p
128 }
129
130 // multiply by the preconditioner if it exists
131 if(M!=Teuchos::null)
132 precOp = Thyra::multiply(precOp,M);
133
134 // must first cast that to be initialized
135 Thyra::DefaultPreconditioner<ScalarT> & dPrec = Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
136
137 // this const-cast is unfortunate...needs to be fixed (larger than it seems!) ECC FIXME!
138 dPrec.initializeUnspecified(Teuchos::rcp_const_cast<Thyra::LinearOpBase<ScalarT> >(precOp));
139}
140
142template <typename ScalarT>
143void NeumannSeriesPreconditionerFactory<ScalarT>::uninitializePrec(Thyra::PreconditionerBase<ScalarT> * prec,
144 RCP<const Thyra::LinearOpSourceBase<ScalarT> > * /* fwdOpSrc */,
145 Thyra::ESupportSolveUse * /* supportSolveUse */) const
146{
147 Thyra::DefaultPreconditioner<ScalarT> & dPrec = Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
148
149 // wipe out any old preconditioner
150 dPrec.uninitialize();
151}
152
153// for ParameterListAcceptor
154
156template <typename ScalarT>
157void NeumannSeriesPreconditionerFactory<ScalarT>::setParameterList(const RCP<Teuchos::ParameterList> & paramList)
158{
159 TEUCHOS_TEST_FOR_EXCEPT(paramList==Teuchos::null);
160
161 // check the parameters are correct
162 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
163
164 // store the parameter list
165 paramList_ = paramList;
166
167 numberOfTerms_ = paramList_->get<int>("Number of Terms");
168
169 // get the scaling type
170 scalingType_ = Teko::NotDiag;
171 const Teuchos::ParameterEntry * entry = paramList_->getEntryPtr("Scaling Type");
172 if(entry!=NULL)
173 scalingType_ = scalingTypeVdtor->getIntegralValue(*entry);
174}
175
177template <typename ScalarT>
178RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getValidParameters() const
179{
180 static RCP<Teuchos::ParameterList> validPL;
181
182 // only fill valid parameter list once
183 if(validPL==Teuchos::null) {
184 RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
185
186 // build the validator for scaling type
187 scalingTypeVdtor = Teuchos::stringToIntegralParameterEntryValidator<DiagonalType>(
188 Teuchos::tuple<std::string>("Diagonal","Lumped","AbsRowSum","None"),
189 Teuchos::tuple<Teko::DiagonalType>(Teko::Diagonal,Teko::Lumped,Teko::AbsRowSum,Teko::NotDiag),
190 "Scaling Type");
191
192 pl->set<int>("Number of Terms",1,
193 "The number of terms to use in the Neumann series expansion.");
194 pl->set("Scaling Type","None","The number of terms to use in the Neumann series expansion.",
195 scalingTypeVdtor);
196
197 validPL = pl;
198 }
199
200 return validPL;
201}
202
204template <typename ScalarT>
205RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::unsetParameterList()
206{
207 Teuchos::RCP<Teuchos::ParameterList> oldList = paramList_;
208 paramList_ = Teuchos::null;
209 return oldList;
210}
211
213template <typename ScalarT>
214RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getParameterList() const
215{
216 return paramList_;
217}
218
220template <typename ScalarT>
221RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getNonconstParameterList()
222{
223 return paramList_;
224}
225
226template <typename ScalarT>
227std::string NeumannSeriesPreconditionerFactory<ScalarT>::description() const
228{
229 std::ostringstream oss;
230 oss << "Teko::NeumannSeriesPreconditionerFactory";
231 return oss.str();
232}
233
234} // end namespace Teko
235
236#endif
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.