Teko  Version of the Day
 All Classes Files Functions Variables Pages
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 
62 namespace Teko {
63 
64 using Teuchos::RCP;
65 using Teuchos::rcp;
66 
67 static RCP<Teuchos::StringToIntegralParameterEntryValidator<Teko::DiagonalType> > scalingTypeVdtor;
68 
69 template <typename ScalarT>
70 NeumannSeriesPreconditionerFactory<ScalarT>::NeumannSeriesPreconditionerFactory()
71  : numberOfTerms_(1), scalingType_(Teko::NotDiag)
72 {
73 }
74 
76 template <typename ScalarT>
77 bool NeumannSeriesPreconditionerFactory<ScalarT>::isCompatible(const Thyra::LinearOpSourceBase<ScalarT> &/* fwdOpSrc */) const
78 {
79  return true;
80 }
81 
83 template <typename ScalarT>
84 RCP<Thyra::PreconditionerBase<ScalarT> > NeumannSeriesPreconditionerFactory<ScalarT>::createPrec() const
85 {
86  return rcp(new Thyra::DefaultPreconditioner<ScalarT>());
87 }
88 
97 template <typename ScalarT>
98 void 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 
142 template <typename ScalarT>
143 void 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 
156 template <typename ScalarT>
157 void 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 
177 template <typename ScalarT>
178 RCP<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 
204 template <typename ScalarT>
205 RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::unsetParameterList()
206 {
207  Teuchos::RCP<Teuchos::ParameterList> oldList = paramList_;
208  paramList_ = Teuchos::null;
209  return oldList;
210 }
211 
213 template <typename ScalarT>
214 RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getParameterList() const
215 {
216  return paramList_;
217 }
218 
220 template <typename ScalarT>
221 RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getNonconstParameterList()
222 {
223  return paramList_;
224 }
225 
226 template <typename ScalarT>
227 std::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