Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_IfpackPreconditionerFactory.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stratimikos: Thyra-based strategies for linear solvers
5 // Copyright (2006) 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 (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
44 #include "Thyra_EpetraLinearOp.hpp"
45 #include "Thyra_DefaultPreconditioner.hpp"
46 #include "Ifpack_ValidParameters.h"
47 #include "Ifpack_Preconditioner.h"
48 #include "Ifpack.h"
49 #include "Epetra_RowMatrix.h"
50 #include "Teuchos_TimeMonitor.hpp"
51 #include "Teuchos_dyn_cast.hpp"
57 
58 
59 namespace {
60 
61 Teuchos::RCP<Teuchos::Time> overallTimer, creationTimer, factorizationTimer;
62 
63 const std::string Ifpack_name = "Ifpack";
64 
65 const std::string IfpackSettings_name = "Ifpack Settings";
66 
67 const std::string PrecType_name = "Prec Type";
69 precTypeValidator; // Will be setup below!
70 const Ifpack::EPrecType PrecType_default = Ifpack::ILU;
71 const std::string PrecTypeName_default = Ifpack::precTypeNames[PrecType_default];
72 
73 const std::string Overlap_name = "Overlap";
74 const int Overlap_default = 0;
75 
76 
78 {
80 }
81 
82 
83 } // namespace
84 
85 namespace Thyra {
86 
87 // Constructors/initializers/accessors
88 
90  :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
91  ,precType_(PrecType_default)
92  ,overlap_(Overlap_default)
93 {
95  getValidParameters(); // Make sure validators get created!
96 }
97 
98 // Overridden from PreconditionerFactoryBase
99 
101  const LinearOpSourceBase<double> &fwdOpSrc
102  ) const
103 {
104  using Teuchos::outArg;
106  EOpTransp epetraFwdOpTransp;
107  EApplyEpetraOpAs epetraFwdOpApplyAs;
108  EAdjointEpetraOp epetraFwdOpAdjointSupport;
109  double epetraFwdOpScalar;
110  epetraFwdOpViewExtractor_->getEpetraOpView(
111  fwdOpSrc.getOp(),
112  outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
113  outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
114  outArg(epetraFwdOpScalar)
115  );
116  if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
117  return false;
118  return true;
119 }
120 
122 {
123  return true;
124 }
125 
127 {
128  return false; // See comment below
129 }
130 
133 {
134  return Teuchos::rcp(new DefaultPreconditioner<double>());
135 }
136 
138  const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc
139  ,PreconditionerBase<double> *prec
140  ,const ESupportSolveUse /* supportSolveUse */
141  ) const
142 {
143  using Teuchos::outArg;
144  using Teuchos::OSTab;
145  using Teuchos::dyn_cast;
146  using Teuchos::RCP;
147  using Teuchos::null;
148  using Teuchos::rcp;
149  using Teuchos::rcp_dynamic_cast;
150  using Teuchos::rcp_const_cast;
151  using Teuchos::set_extra_data;
152  using Teuchos::get_optional_extra_data;
154  Teuchos::Time totalTimer(""), timer("");
155  totalTimer.start(true);
156 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
157  Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
158 #endif
159  const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
160  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
161  Teuchos::OSTab tab(out);
162  if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
163  *out << "\nEntering Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
164 #ifdef TEUCHOS_DEBUG
165  TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
166  TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
167 #endif
169  fwdOp = fwdOpSrc->getOp();
170 #ifdef TEUCHOS_DEBUG
171  TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
172 #endif
173  //
174  // Unwrap and get the forward Epetra_Operator object
175  //
177  EOpTransp epetraFwdOpTransp;
178  EApplyEpetraOpAs epetraFwdOpApplyAs;
179  EAdjointEpetraOp epetraFwdOpAdjointSupport;
180  double epetraFwdOpScalar;
181  epetraFwdOpViewExtractor_->getEpetraOpView(
182  fwdOp,
183  outArg(epetraFwdOp), outArg(epetraFwdOpTransp),
184  outArg(epetraFwdOpApplyAs), outArg(epetraFwdOpAdjointSupport),
185  outArg(epetraFwdOpScalar)
186  );
187  // Validate what we get is what we need
189  epetraFwdRowMat = rcp_dynamic_cast<const Epetra_RowMatrix>(epetraFwdOp,true);
191  epetraFwdOpApplyAs != EPETRA_OP_APPLY_APPLY, std::logic_error
192  ,"Error, incorrect apply mode for an Epetra_RowMatrix"
193  );
194  //
195  // Get the concrete preconditioner object
196  //
197  DefaultPreconditioner<double>
198  *defaultPrec = &Teuchos::dyn_cast<DefaultPreconditioner<double> >(*prec);
199  //
200  // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
201  //
203  epetra_precOp = rcp_dynamic_cast<EpetraLinearOp>(defaultPrec->getNonconstUnspecifiedPrecOp(),true);
204  //
205  // Get the embedded Ifpack_Preconditioner object if it exists
206  //
208  ifpack_precOp;
209  if(epetra_precOp.get())
210  ifpack_precOp = rcp_dynamic_cast<Ifpack_Preconditioner>(epetra_precOp->epetra_op(),true);
211  //
212  // Get the attached forward operator if it exists and make sure that it matches
213  //
214  if(ifpack_precOp.get()) {
215  // ToDo: Get the forward operator and make sure that it matches what is
216  // already being used!
217  }
218  //
219  // Permform initialization if needed
220  //
221  //const bool startingOver = (ifpack_precOp.get() == NULL);
222  const bool startingOver = true;
223  // ToDo: Comment back in the above original version of startingOver to allow
224  // for resuse. Rob H. just pointed out to me that a memory leak is being
225  // created when you just call Ifpack_ILU::Compute() over and over again.
226  // Rob H. said that he will check in a fix the the development branch when
227  // he can.
228  if(startingOver) {
229  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
230  *out << "\nCreating the initial Ifpack_Preconditioner object of type \'"<<Ifpack::toString(precType_)<<"\' ...\n";
231  timer.start(true);
232 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
233  Teuchos::TimeMonitor creationTimeMonitor(*creationTimer);
234 #endif
235  // Create the initial preconditioner
236  ifpack_precOp = rcp(
237  ::Ifpack::Create(
238  precType_
239  ,const_cast<Epetra_RowMatrix*>(&*epetraFwdRowMat)
240  ,overlap_
241  )
242  );
243  timer.stop();
244  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
245  OSTab(out).o() <<"=> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
246  // Set parameters if the list exists
247  if(paramList_.get()) {
249  &ifpackSettingsPL = paramList_->sublist(IfpackSettings_name);
250  // Above will create new sublist if it does not exist!
251  TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->SetParameters(ifpackSettingsPL));
252  // Above, I have not idea how any error messages for a mistake will be
253  // reported back to the user!
254  }
255  // Initialize the structure for the preconditioner
256  TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Initialize());
257  }
258  //
259  // Attach the epetraFwdOp to the ifpack_precOp to guarantee that it will not go away
260  //
261  set_extra_data(epetraFwdOp, "IFPF::epetraFwdOp", Teuchos::inOutArg(ifpack_precOp),
262  Teuchos::POST_DESTROY, false);
263  //
264  // Update the factorization
265  //
266  {
267  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
268  *out << "\nComputing the preconditioner ...\n";
269 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
270  Teuchos::TimeMonitor factorizationTimeMonitor(*factorizationTimer);
271 #endif
272  timer.start(true);
273  TEUCHOS_TEST_FOR_EXCEPT(0!=ifpack_precOp->Compute());
274  timer.stop();
275  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
276  OSTab(out).o() <<"=> Setup time = "<<timer.totalElapsedTime()<<" sec\n";
277  }
278  //
279  // Compute the conditioner number estimate if asked
280  //
281 
282  // ToDo: Implement
283 
284  //
285  // Attach fwdOp to the ifpack_precOp
286  //
287  set_extra_data(fwdOpSrc, "IFPF::fwdOpSrc", Teuchos::inOutArg(ifpack_precOp),
288  Teuchos::POST_DESTROY, false);
289  //
290  // Initialize the output EpetraLinearOp
291  //
292  if(startingOver) {
293  epetra_precOp = rcp(new EpetraLinearOp);
294  }
295  epetra_precOp->initialize(
296  ifpack_precOp
297  ,epetraFwdOpTransp
300  );
301  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_MEDIUM)) {
302  *out << "\nDescription of created preconditioner:\n";
303  OSTab tab2(out);
304  ifpack_precOp->Print(*out);
305  }
306 
307  //
308  // Initialize the preconditioner
309  //
310  defaultPrec->initializeUnspecified(
311  Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp)
312  );
313  totalTimer.stop();
314  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
315  *out << "\nTotal time in IfpackPreconditionerFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
316  if(out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW))
317  *out << "\nLeaving Thyra::IfpackPreconditionerFactory::initializePrec(...) ...\n";
318 }
319 
321  PreconditionerBase<double> * /* prec */
322  ,Teuchos::RCP<const LinearOpSourceBase<double> > * /* fwdOpSrc */
323  ,ESupportSolveUse * /* supportSolveUse */
324  ) const
325 {
326  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
327 }
328 
329 // Overridden from ParameterListAcceptor
330 
332 {
333  TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
334  paramList->validateParameters(*this->getValidParameters(),2);
335  // Note: The above validation will validate right down into the the sublist
336  // named IfpackSettings_name!
337  paramList_ = paramList;
338  overlap_ = paramList_->get(Overlap_name,Overlap_default);
339  std::ostringstream oss;
340  oss << "(sub)list \""<<paramList->name()<<"\"parameter \"Prec Type\"";
341  precType_ =
342  ( paramList_.get()
343  ? precTypeValidator->getIntegralValue(*paramList_,PrecType_name,PrecTypeName_default)
344  : PrecType_default
345  );
346  Teuchos::readVerboseObjectSublist(&*paramList_,this);
347 #ifdef TEUCHOS_DEBUG
348  // Validate my use of the parameters!
349  paramList->validateParameters(*this->getValidParameters(),1);
350 #endif
351 }
352 
355 {
356  return paramList_;
357 }
358 
361 {
364  return _paramList;
365 }
366 
369 {
370  return paramList_;
371 }
372 
375 {
376  using Teuchos::rcp;
377  using Teuchos::rcp_implicit_cast;
379  static Teuchos::RCP<Teuchos::ParameterList> validParamList;
380  if(validParamList.get()==NULL) {
381  validParamList = Teuchos::rcp(new Teuchos::ParameterList(Ifpack_name));
382  {
383  // Create the validator for the preconditioner type!
385  precTypeNames;
386  precTypeNames.insert(
387  precTypeNames.begin(),
388  &Ifpack::precTypeNames[0],
389  &Ifpack::precTypeNames[0] + Ifpack::numPrecTypes
390  );
392  precTypeValues;
393  precTypeValues.insert(
394  precTypeValues.begin(),
395  &Ifpack::precTypeValues[0],
396  &Ifpack::precTypeValues[0] + Ifpack::numPrecTypes
397  );
398  precTypeValidator = rcp(
400  precTypeNames,precTypeValues,PrecType_name
401  )
402  );
403  }
404  validParamList->set(
405  PrecType_name, PrecTypeName_default,
406  "Type of Ifpack preconditioner to use.",
407  rcp_implicit_cast<const PEV>(precTypeValidator)
408  );
409  validParamList->set(
410  Overlap_name, Overlap_default,
411  "Number of rows/columns overlapped between subdomains in different"
412  "\nprocesses in the additive Schwarz-type domain-decomposition preconditioners."
413  );
414  validParamList->sublist(
415  IfpackSettings_name, false,
416  "Preconditioner settings that are passed onto the Ifpack preconditioners themselves."
417  ).setParameters(Ifpack_GetValidParameters());
418  // Note that in the above setParameterList(...) function that we actually
419  // validate down into the first level of this sublist. Really the
420  // Ifpack_Preconditioner objects themselves should do validation but we do
421  // it ourselves taking the return from the Ifpack_GetValidParameters()
422  // function as gospel!
423  Teuchos::setupVerboseObjectSublist(&*validParamList);
424  }
425  return validParamList;
426 }
427 
428 // Public functions overridden from Teuchos::Describable
429 
431 {
432  std::ostringstream oss;
433  oss << "Thyra::IfpackPreconditionerFactory{";
434  oss << "precType=\"" << ::Ifpack::toString(precType_) << "\"";
435  oss << ",overlap=" << overlap_;
436  oss << "}";
437  return oss.str();
438 }
439 
440 // private
441 
443 {
444  if(!overallTimer.get()) {
445 #ifdef STRATIMIKOS_TEUCHOS_TIME_MONITOR
446  overallTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF");
447  creationTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Creation");
448  factorizationTimer = Teuchos::TimeMonitor::getNewTimer("Stratimikos: IfpackPF:Factorization");
449 #endif
450  }
451 }
452 
453 } // namespace Thyra
const std::string & name() const
const std::string toString(const EAdjointEpetraOp adjointEpetraOp)
Teuchos::RCP< Teuchos::ParameterList > paramList_
Teuchos::RCP< PreconditionerBase< double > > createPrec() const
T & get(ParameterList &l, const std::string &name)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool isCompatible(const LinearOpSourceBase< double > &fwdOpSrc) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
T * get() const
void uninitializePrec(PreconditionerBase< double > *prec, Teuchos::RCP< const LinearOpSourceBase< double > > *fwdOpSrc, ESupportSolveUse *supportSolveUse) const
basic_OSTab< char > OSTab
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
static RCP< Time > getNewTimer(const std::string &name)
IntegralType getIntegralValue(ParameterList const &paramList, std::string const &paramName)
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double stop()
EAdjointEpetraOp
T_To & dyn_cast(T_From &from)
TypeTo implicit_cast(const TypeFrom &t)
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
#define TEUCHOS_ADD_STRINGTOINTEGRALVALIDATOR_CONVERTER(INTEGRALTYPE)
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< double > > &fwdOpSrc, PreconditionerBase< double > *prec, const ESupportSolveUse supportSolveUse) const
double totalElapsedTime(bool readCurrentTime=false) const
EApplyEpetraOpAs
iterator insert(iterator position, const value_type &x)
iterator begin()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
TEUCHOS_STATIC_SETUP()