Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_SmootherPreconditionerFactory.cpp
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 // Teko includes
48 #include "Teko_SmootherPreconditionerFactory.hpp"
49 
50 #include "Teko_PreconditionerInverseFactory.hpp"
51 #include "Thyra_MultiVectorStdOps.hpp"
52 
53 namespace Teko {
54 
55 // A small utility class to help distinguish smoothers
56 class SmootherRequestMesg : public RequestMesg {
57  public:
58  SmootherRequestMesg(unsigned int block)
59  : RequestMesg("__smoother_request_message__"), block_(block) {}
60 
61  unsigned int getBlock() const { return block_; }
62 
63  private:
64  unsigned int block_;
65 };
66 
70 
71 SmootherLinearOp::SmootherLinearOp(const LinearOp& A, const LinearOp& invM,
72  unsigned int applications, bool useDestAsInitialGuess)
73  : A_(A),
74  invM_(invM),
75  applications_(applications),
76  initialGuessType_(Unspecified),
77  requestMesg_(Teuchos::null) {
78  // set initial guess type
79  initialGuessType_ = useDestAsInitialGuess ? DestAsInitialGuess : NoInitialGuess;
80 }
81 
82 SmootherLinearOp::SmootherLinearOp(const LinearOp& A, const LinearOp& invM,
83  unsigned int applications, unsigned int block)
84  : A_(A),
85  invM_(invM),
86  applications_(applications),
87  initialGuessType_(RequestInitialGuess),
88  requestMesg_(Teuchos::null) {
89  requestMesg_ = Teuchos::rcp(new SmootherRequestMesg(block));
90 }
91 
92 void SmootherLinearOp::implicitApply(const MultiVector& b, MultiVector& x, const double alpha,
93  const double beta) const {
94  using Teuchos::RCP;
95 
96  MultiVector residual = deepcopy(b); // residual = b
97  MultiVector scrap = deepcopy(b); // scrap = b
98  MultiVector error; // location for initial guess
99 
100  // construct initial guess: required to assign starting point for destination
101  // vector appropriately
102  switch (initialGuessType_) {
103  case RequestInitialGuess:
104  // get initial guess from request handler
105  error = deepcopy(getRequestHandler()->request<MultiVector>(*requestMesg_));
106  Thyra::assign<double>(x.ptr(), *error); // x = error (initial guess)
107  break;
108  case DestAsInitialGuess:
109  error = deepcopy(x); // error = x
110  break;
111  case NoInitialGuess:
112  Thyra::assign<double>(x.ptr(), 0.0); // x = 0
113  error = deepcopy(x); // error = x
114  break;
115  case Unspecified:
116  default: TEUCHOS_ASSERT(false);
117  }
118 
119  for (unsigned int current = 0; current < applications_; ++current) {
120  // compute current residual
121  Teko::applyOp(A_, error, scrap);
122  Teko::update(-1.0, scrap, 1.0, residual); // residual = residual-A*error
123 
124  // compute appoximate correction using residual
125  Thyra::assign(error.ptr(), 0.0); // set error to zero
126  Teko::applyOp(invM_, residual, error);
127 
128  // update solution with error
129  Teko::update(1.0, error, 1.0, x); // x = x+error
130  }
131 }
132 
134 void SmootherLinearOp::setRequestHandler(const Teuchos::RCP<RequestHandler>& rh) {
135  Teko_DEBUG_SCOPE("SmootherLinearOp::setRequestHandler", 10);
136  requestHandler_ = rh;
137 }
138 
140 Teuchos::RCP<RequestHandler> SmootherLinearOp::getRequestHandler() const {
141  Teko_DEBUG_SCOPE("SmootherLinearOp::getRequestHandler", 10);
142  return requestHandler_;
143 }
144 
145 LinearOp buildSmootherLinearOp(const LinearOp& A, const LinearOp& invM, unsigned int applications,
146  bool useDestAsInitialGuess) {
147  return Teuchos::rcp(new SmootherLinearOp(A, invM, applications, useDestAsInitialGuess));
148 }
149 
150 LinearOp buildSmootherLinearOp(const LinearOp& A, const LinearOp& invM, unsigned int applications,
151  unsigned int initialGuessBlock) {
152  return Teuchos::rcp(new SmootherLinearOp(A, invM, applications, initialGuessBlock));
153 }
154 
158 
160 SmootherPreconditionerFactory::SmootherPreconditionerFactory()
161  : sweepCount_(0),
162  initialGuessType_(Unspecified),
163  initialGuessBlock_(0),
164  precFactory_(Teuchos::null) {}
165 
169 LinearOp SmootherPreconditionerFactory::buildPreconditionerOperator(
170  LinearOp& lo, PreconditionerState& state) const {
171  TEUCHOS_TEST_FOR_EXCEPTION(
172  precFactory_ == Teuchos::null, std::runtime_error,
173  "ERROR: Teko::SmootherPreconditionerFactory::buildPreconditionerOperator requires that a "
174  << "preconditioner factory has been set. Currently it is null!");
175 
176  // preconditions
177  TEUCHOS_ASSERT(sweepCount_ > 0);
178  TEUCHOS_ASSERT(initialGuessType_ != Unspecified);
179  TEUCHOS_ASSERT(precFactory_ != Teuchos::null);
180 
181  // build user specified preconditioner
182  ModifiableLinearOp& invM = state.getModifiableOp("prec");
183  if (invM == Teuchos::null)
184  invM = Teko::buildInverse(*precFactory_, lo);
185  else
186  Teko::rebuildInverse(*precFactory_, lo, invM);
187 
188  // conditional on initial guess type, build the smoother
189  switch (initialGuessType_) {
190  case RequestInitialGuess:
191  return buildSmootherLinearOp(lo, invM, sweepCount_, initialGuessBlock_);
192  case DestAsInitialGuess:
193  return buildSmootherLinearOp(lo, invM, sweepCount_, true); // use an initial guess
194  case NoInitialGuess:
195  return buildSmootherLinearOp(lo, invM, sweepCount_, false); // no initial guess
196  case Unspecified:
197  default: TEUCHOS_ASSERT(false);
198  }
199 
200  // should never get here
201  TEUCHOS_ASSERT(false);
202  return Teuchos::null;
203 }
204 
208 void SmootherPreconditionerFactory::initializeFromParameterList(
209  const Teuchos::ParameterList& settings) {
210  // declare all strings used by this initialization routine
212 
213  const std::string str_sweepCount = "Sweep Count";
214  const std::string str_initialGuessBlock = "Initial Guess Block";
215  const std::string str_destAsInitialGuess = "Destination As Initial Guess";
216  const std::string str_precType = "Preconditioner Type";
217 
218  // default parameters
220 
221  initialGuessType_ = Unspecified;
222  initialGuessBlock_ = 0;
223  sweepCount_ = 0;
224  precFactory_ = Teuchos::null;
225 
226  // get sweep counts
228 
229  if (settings.isParameter(str_sweepCount)) sweepCount_ = settings.get<int>(str_sweepCount);
230 
231  // get initial guess
233 
234  if (settings.isParameter(str_initialGuessBlock)) {
235  initialGuessBlock_ = settings.get<int>(str_initialGuessBlock);
236  initialGuessType_ = RequestInitialGuess;
237  }
238 
239  if (settings.isParameter(str_destAsInitialGuess)) {
240  bool useDest = settings.get<bool>(str_destAsInitialGuess);
241  if (useDest) {
242  TEUCHOS_TEST_FOR_EXCEPTION(initialGuessType_ != Unspecified, std::runtime_error,
243  "Cannot set both \"" << str_initialGuessBlock << "\" and \""
244  << str_destAsInitialGuess << "\"");
245 
246  initialGuessType_ = DestAsInitialGuess;
247  }
248  }
249 
250  // safe to assume if the other values are not turned on there is no initial guess
251  if (initialGuessType_ == Unspecified) initialGuessType_ = NoInitialGuess;
252 
253  // get preconditioner factory
255 
256  TEUCHOS_TEST_FOR_EXCEPTION(
257  not settings.isParameter(str_precType), std::runtime_error,
258  "Parameter \"" << str_precType << "\" is required by a Teko::SmootherPreconditionerFactory");
259 
260  // grab library and preconditioner name
261  Teuchos::RCP<const InverseLibrary> il = getInverseLibrary();
262  std::string precName = settings.get<std::string>(str_precType);
263 
264  // build preconditioner factory
265  precFactory_ = il->getInverseFactory(precName);
266  TEUCHOS_TEST_FOR_EXCEPTION(
267  precFactory_ == Teuchos::null, std::runtime_error,
268  "ERROR: \"" << str_precType << "\" = " << precName << " could not be found");
269 
270  // post conditions required to be satisfied
272 
273  TEUCHOS_ASSERT(sweepCount_ > 0);
274  TEUCHOS_ASSERT(initialGuessType_ != Unspecified);
275  TEUCHOS_ASSERT(precFactory_ != Teuchos::null);
276 }
277 
278 } // end namespace Teko
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
virtual void setRequestHandler(const Teuchos::RCP< RequestHandler > &rh)
Set the request handler with pointers to the appropriate callbacks.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
virtual Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
virtual void implicitApply(const MultiVector &x, MultiVector &y, const double alpha=1.0, const double beta=0.0) const
Perform a matrix vector multiply with this implicitly defined blocked operator.