Thyra  Version of the Day
Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
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 // 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 Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
43 #define THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
44 
45 
46 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory_decl.hpp"
47 #include "Thyra_LinearOpWithSolveBase.hpp"
48 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
49 #include "Thyra_PhysicallyBlockedLinearOpBase.hpp"
50 #include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
51 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
52 #include "Thyra_DefaultLinearOpSource.hpp"
53 
54 
55 namespace Thyra {
56 
57 
58 // Overridden from Constructors/Initializers/Accessors
59 
60 
61 template<class Scalar>
63  const RCP<LinearOpWithSolveFactoryBase<Scalar> > &lowsf
64  )
65 {
66 #ifdef TEUCHOS_DEBUG
67  TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf));
68 #endif
69  lowsf_.initialize(lowsf);
70 }
71 
72 
73 template<class Scalar>
75  const RCP<const LinearOpWithSolveFactoryBase<Scalar> > &lowsf
76  )
77 {
78 #ifdef TEUCHOS_DEBUG
79  TEUCHOS_TEST_FOR_EXCEPT(is_null(lowsf));
80 #endif
81  lowsf_.initialize(lowsf);
82 }
83 
84 
85 template<class Scalar>
86 RCP<LinearOpWithSolveFactoryBase<Scalar> >
88 {
89  return lowsf_.getNonconstObj();
90 }
91 
92 
93 template<class Scalar>
94 RCP<const LinearOpWithSolveFactoryBase<Scalar> >
96 {
97  return lowsf_.getConstObj();
98 }
99 
100 
101 // Overridden from Teuchos::Describable
102 
103 
104 template<class Scalar>
105 std::string
107 {
108  std::ostringstream oss;
109  oss << this->Teuchos::Describable::description()
110  << "{"
111  << "lowsf=";
112  if (!is_null(lowsf_.getConstObj()))
113  oss << lowsf_.getConstObj()->description();
114  else
115  oss << "NULL";
116  oss << "}";
117  return oss.str();
118 }
119 
120 
121 // Overridden from ParameterListAcceptor
122 
123 
124 template<class Scalar>
125 void
127  RCP<ParameterList> const& paramList
128  )
129 {
130  lowsf_.getNonconstObj()->setParameterList(paramList);
131 }
132 
133 
134 template<class Scalar>
135 RCP<ParameterList>
137 {
138  return lowsf_.getNonconstObj()->getNonconstParameterList();
139 }
140 
141 
142 template<class Scalar>
143 RCP<ParameterList>
145 {
146  return lowsf_.getNonconstObj()->unsetParameterList();
147 }
148 
149 
150 template<class Scalar>
151 RCP<const ParameterList>
153 {
154  return lowsf_.getConstObj()->getParameterList();
155 }
156 
157 
158 template<class Scalar>
159 RCP<const ParameterList>
161 {
162  return lowsf_.getConstObj()->getValidParameters();
163 }
164 
165 
166 // Overridden from LinearOpWithSolveFactoyBase
167 
168 
169 template<class Scalar>
170 bool
172 {
173  return false;
174 }
175 
176 
177 template<class Scalar>
178 void
180  const RCP<PreconditionerFactoryBase<Scalar> > &precFactory,
181  const std::string &precFactoryName
182  )
183 {
184  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
185  "Error, we don't support a preconditioner factory!");
186 }
187 
188 
189 template<class Scalar>
190 RCP<PreconditionerFactoryBase<Scalar> >
192 {
193  return Teuchos::null;
194 }
195 
196 
197 template<class Scalar>
199  RCP<PreconditionerFactoryBase<Scalar> > *precFactory,
200  std::string *precFactoryName
201  )
202 {
203  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
204  "Error, we don't support a preconditioner factory!");
205 }
206 
207 
208 template<class Scalar>
209 bool
211  const LinearOpSourceBase<Scalar> &fwdOpSrc
212  ) const
213 {
214  TEUCHOS_TEST_FOR_EXCEPT(true);
215  return false;
216 }
217 
218 
219 template<class Scalar>
220 RCP<LinearOpWithSolveBase<Scalar> >
222 {
223  return defaultBlockedTriangularLinearOpWithSolve<Scalar>();
224 }
225 
226 
227 template<class Scalar>
228 void
230  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
232  const ESupportSolveUse supportSolveUse
233  ) const
234 {
235 
236  using Teuchos::dyn_cast;
237  using Teuchos::rcp_dynamic_cast;
238 
239 #ifdef TEUCHOS_DEBUG
240  TEUCHOS_TEST_FOR_EXCEPT(0==Op);
241 #endif
242 
243  // Set the verbosity settings for the wrapped LOWSF object!
244  lowsf_.getConstObj()->setOStream(this->getOStream());
245  lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel());
246 
247  // Get the block interface to get at the blocks
249  const RCP<const PBLOB> blo =
250  rcp_dynamic_cast<const PBLOB>(fwdOpSrc->getOp().assert_not_null());
251 
252  // Dynamic cast to get the DefaultBlockedTriangularLinearOpWithSolveBase
253  // interface that we will fill.
254 
256  DBTLOWS &btlows = dyn_cast<DBTLOWS>(*Op);
257 
258  // Determine if this is the first time through or if we have already
259  // initialized before. This will be needed to allow efficient reuse of the
260  // LOWSB objects for the diagonal blocks.
261  const bool firstTime = is_null(btlows.range());
262 
263  // If this is the first time through, we need to fill and create the block
264  // structure
265  if (firstTime)
266  btlows.beginBlockFill(blo->productRange(),blo->productDomain());
267 
268  const int N = blo->productRange()->numBlocks();
269  for ( int k = 0; k < N; ++k ) {
270  const RCP<const LinearOpBase<Scalar> > fwdOp_k =
271  blo->getBlock(k,k).assert_not_null();
272  if (firstTime) {
273  // This is the first time through so reate and initialize a new LOWSB
274  // object for each block
275  btlows.setNonconstLOWSBlock( k, k,
276  linearOpWithSolve<Scalar>(*lowsf_.getConstObj(),fwdOp_k)
277  );
278  }
279  else {
280  // This is not the first time through so we need to just reinitiallize
281  // the object that is already created. This allows us to efficiently
282  // reuse precreated structure and storage.
283  RCP<LinearOpWithSolveBase<Scalar> >
284  invOp_k = btlows.getNonconstLOWSBlock(k,k).assert_not_null();
285  Thyra::initializeOp<Scalar>(*lowsf_.getConstObj(), fwdOp_k, invOp_k.ptr());
286  }
287  }
288 
289  // If this is the first time through, then we need to finalize the block
290  // structure.
291  if (firstTime)
292  btlows.endBlockFill();
293 
294  // After the block structure has been setup, set the off-diagonal blocks.
295  // Note that this also sets the diagonal blocks but these are ignored since
296  // the LOWSB blocks created above override these.
297  btlows.setBlocks(blo);
298 
299  // Set the verbosity settings
300  btlows.setOStream(this->getOStream());
301  btlows.setVerbLevel(this->getVerbLevel());
302 
303 }
304 
305 
306 template<class Scalar>
307 void
309  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
311  ) const
312 {
313  TEUCHOS_TEST_FOR_EXCEPT(true);
314 }
315 
316 
317 template<class Scalar>
318 void
321  RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
322  RCP<const PreconditionerBase<Scalar> > *prec,
323  RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
324  ESupportSolveUse *supportSolveUse
325  ) const
326 {
327  using Teuchos::dyn_cast;
328  using Teuchos::rcp_implicit_cast;
329  using Teuchos::rcp_dynamic_cast;
331  TEUCHOS_TEST_FOR_EXCEPT(0==Op);
332  DBTLOWS &btlowsOp = dyn_cast<DBTLOWS>(*Op);
333  if (fwdOpSrc) {
334  const RCP<const LinearOpBase<Scalar> > fwdOp = btlowsOp.getBlocks();
335  if (!is_null(fwdOp))
336  *fwdOpSrc = defaultLinearOpSource<Scalar>(fwdOp);
337  else
338  *fwdOpSrc = Teuchos::null;
339  }
340  if (prec) *prec = Teuchos::null;
341  if (approxFwdOpSrc) *approxFwdOpSrc = Teuchos::null;
342 }
343 
344 
345 template<class Scalar>
346 bool
348  const EPreconditionerInputType precOpType
349  ) const
350 {
351  // We don't support any external preconditioners!
352  return false;
353  // 20071006: rabartl: Note: We could support external preconditioners but it
354  // will take some work. We would have to extract out the individual
355  // preconditioners from each block. This would be pretty easy to do but I
356  // am not going to do this until we have to.
357 }
358 
359 
360 template<class Scalar>
361 void
363  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
364  const RCP<const PreconditionerBase<Scalar> > &prec,
366  const ESupportSolveUse supportSolveUse
367  ) const
368 {
369  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
370  "Error, we don't support an external preconditioner!");
371 }
372 
373 
374 template<class Scalar>
375 void
377  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
378  const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
380  const ESupportSolveUse supportSolveUse
381  ) const
382 {
383  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
384  "Error, we don't support an external preconditioner!");
385 }
386 
387 
388 // protected
389 
390 
391 template<class Scalar>
392 void
394 {
395  lowsf_.getConstObj()->setVerbLevel(this->getVerbLevel());
396  lowsf_.getConstObj()->setOStream(this->getOStream());
397 }
398 
399 
400 } // namespace Thyra
401 
402 
403 #endif // THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
virtual RCP< PreconditionerFactoryBase< Scalar > > getPreconditionerFactory() const
Returns null .
Base class for all linear operators that can support a high-level solve operation.
virtual void initializeAndReuseOp(const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op) const
Concrete composite LinearOpWithSolveBase subclass that creates single upper or lower block triangular...
virtual void unsetPreconditionerFactory(RCP< PreconditionerFactoryBase< Scalar > > *precFactory, std::string *precFactoryName)
Throws exception.
Simple interface class to access a precreated preconditioner as one or more linear operators objects ...
Implicit subclass that takes a blocked triangular LOWB object and turns it into a LOWSB object...
virtual void initializePreconditionedOp(const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const RCP< const PreconditionerBase< Scalar > > &prec, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects...
Base interface for physically blocked linear operators.
Factory interface for creating preconditioner objects from LinearOpBase objects.
Base interface for objects that can return a linear operator.
ESupportSolveUse
Enum that specifies how a LinearOpWithSolveBase object will be used for solves after it is constructe...
virtual void setPreconditionerFactory(const RCP< PreconditionerFactoryBase< Scalar > > &precFactory, const std::string &precFactoryName)
Throws exception.
virtual void uninitializeOp(LinearOpWithSolveBase< Scalar > *Op, RCP< const LinearOpSourceBase< Scalar > > *fwdOpSrc, RCP< const PreconditionerBase< Scalar > > *prec, RCP< const LinearOpSourceBase< Scalar > > *approxFwdOpSrc, ESupportSolveUse *supportSolveUse) const
virtual void initializeApproxPreconditionedOp(const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
EPreconditionerInputType
Enum defining the status of a preconditioner object.
virtual bool supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
virtual void initializeOp(const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const