Teko  Version of the Day
Teko_StratimikosFactory.cpp
1 #include "Teko_StratimikosFactory.hpp"
2 
3 #include "Teuchos_Time.hpp"
4 #include "Teuchos_AbstractFactoryStd.hpp"
5 
6 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
7 #include "Thyra_EpetraLinearOp.hpp"
8 #include "Thyra_DefaultPreconditioner.hpp"
9 
10 #include "Teko_InverseLibrary.hpp"
11 #include "Teko_Preconditioner.hpp"
12 #include "Teko_InverseFactoryOperator.hpp"
13 #include "Teko_Utilities.hpp"
14 #include "Teko_InverseLibrary.hpp"
15 #include "Teko_StridedEpetraOperator.hpp"
16 #include "Teko_BlockedEpetraOperator.hpp"
17 #include "Teko_ReorderedLinearOp.hpp"
18 
19 #include "EpetraExt_RowMatrixOut.h"
20 
21 namespace Teko {
22 
23 using Teuchos::RCP;
24 using Teuchos::ParameterList;
25 
26 // hide stuff
27 namespace {
28  // Simple preconditioner class that adds a counter
29  class StratimikosFactoryPreconditioner : public Thyra::DefaultPreconditioner<double>{
30  public:
31  StratimikosFactoryPreconditioner() : iter_(0) {}
32 
33  inline void incrIter() { iter_++; }
34  inline std::size_t getIter() { return iter_; }
35 
36  private:
37  StratimikosFactoryPreconditioner(const StratimikosFactoryPreconditioner &);
38 
39  std::size_t iter_;
40  };
41 
42  // factory used to initialize the Teko::StratimikosFactory
43  // user data
44  class TekoFactoryBuilder
45  : public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
46  public:
47  TekoFactoryBuilder(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
48  const Teuchos::RCP<Teko::RequestHandler> & rh) : builder_(builder), requestHandler_(rh) {}
49  Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create() const
50  { return Teuchos::rcp(new StratimikosFactory(builder_,requestHandler_)); }
51 
52  private:
53  Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builder_;
54  Teuchos::RCP<Teko::RequestHandler> requestHandler_;
55  };
56 }
57 
58 // Constructors/initializers/accessors
60  :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
61 {}
62 
63 // Constructors/initializers/accessors
64 StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Teko::RequestHandler> & rh)
65  :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
66 {
68 }
69 
70 StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
71  const Teuchos::RCP<Teko::RequestHandler> & rh)
72  :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())), builder_(builder)
73 {
75 }
76 
77 
78 // Overridden from PreconditionerFactoryBase
79 bool StratimikosFactory::isCompatible(const Thyra::LinearOpSourceBase<double> &fwdOpSrc) const
80 {
81  using Teuchos::outArg;
82 
83  TEUCHOS_ASSERT(false); // what you doing?
84 
85  Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
86  Thyra::EOpTransp epetraFwdOpTransp;
87  Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
88  Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
89  double epetraFwdOpScalar;
90 
91  // check to make sure this is an epetra CrsMatrix
92  Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc.getOp();
93  epetraFwdOpViewExtractor_->getEpetraOpView(
94  fwdOp,
95  outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
96  outArg(epetraFwdOpApplyAs),
97  outArg(epetraFwdOpAdjointSupport),
98  outArg(epetraFwdOpScalar));
99 
100  if( !dynamic_cast<const Epetra_CrsMatrix*>(&*epetraFwdOp) )
101  return false;
102 
103  return true;
104 }
105 
106 
107 bool StratimikosFactory::applySupportsConj(Thyra::EConj conj) const
108 {
109  return false;
110 }
111 
112 
114 {
115  return false; // See comment below
116 }
117 
118 
119 Teuchos::RCP<Thyra::PreconditionerBase<double> >
121 {
122  return Teuchos::rcp(new StratimikosFactoryPreconditioner());
123 }
124 
125 
127  const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
128  Thyra::PreconditionerBase<double> *prec,
129  const Thyra::ESupportSolveUse supportSolveUse
130  ) const
131 {
132  Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
133 
134  if(epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
135  initializePrec_Epetra(fwdOpSrc,prec,supportSolveUse);
136  else
137  initializePrec_Thyra(fwdOpSrc,prec,supportSolveUse);
138 }
139 
141  const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
142  Thyra::PreconditionerBase<double> *prec,
143  const Thyra::ESupportSolveUse supportSolveUse
144  ) const
145 {
146  using Teuchos::RCP;
147  using Teuchos::rcp;
148  using Teuchos::rcp_dynamic_cast;
149  using Teuchos::implicit_cast;
150  using Thyra::LinearOpBase;
151 
152  Teuchos::Time totalTimer(""), timer("");
153  totalTimer.start(true);
154 
155  const RCP<Teuchos::FancyOStream> out = this->getOStream();
156  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
157 
158  bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
159 
160  Teuchos::OSTab tab(out);
161  if(mediumVerbosity)
162  *out << "\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
163 
164  Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
165 
166  // Get the concrete preconditioner object
167  StratimikosFactoryPreconditioner & defaultPrec
168  = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
169  Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
170 
171  // Permform initialization if needed
172  const bool startingOver = (prec_Op == Teuchos::null);
173  if(startingOver) {
174  // start over
175  invLib_ = Teuchos::null;
176  invFactory_ = Teuchos::null;
177 
178  if(mediumVerbosity)
179  *out << "\nCreating the initial Teko Operator object...\n";
180 
181  timer.start(true);
182 
183  // build library, and set request handler (user defined!)
184  invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"),builder_);
185  invLib_->setRequestHandler(reqHandler_);
186 
187  // build preconditioner factory
188  invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
189 
190  timer.stop();
191  if(mediumVerbosity)
192  Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
193  }
194 
195  if(mediumVerbosity)
196  *out << "\nComputing the preconditioner ...\n";
197 
198  // setup reordering if required
199  std::string reorderType = paramList_->get<std::string>("Reorder Type");
200  if(reorderType!="") {
201 
202  Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
203  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(fwdOp,true);
204  RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
205  Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm,blkFwdOp);
206 
207  if(prec_Op==Teuchos::null) {
208  Teko::ModifiableLinearOp reorderedPrec = Teko::buildInverse(*invFactory_,blockedFwdOp);
209  prec_Op = Teuchos::rcp(new ReorderedLinearOp(brm,reorderedPrec));
210  }
211  else {
212  Teko::ModifiableLinearOp reorderedPrec = Teuchos::rcp_dynamic_cast<ReorderedLinearOp>(prec_Op,true)->getBlockedOp();
213  Teko::rebuildInverse(*invFactory_,blockedFwdOp,reorderedPrec);
214  }
215  }
216  else {
217  // no reordering required
218  if(prec_Op==Teuchos::null)
219  prec_Op = Teko::buildInverse(*invFactory_,fwdOp);
220  else
221  Teko::rebuildInverse(*invFactory_,fwdOp,prec_Op);
222  }
223 
224  // construct preconditioner
225  timer.stop();
226 
227  if(mediumVerbosity)
228  Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
229 
230 
231  // Initialize the preconditioner
232  defaultPrec.initializeUnspecified( prec_Op);
233  defaultPrec.incrIter();
234  totalTimer.stop();
235 
236  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
237  *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
238  if(mediumVerbosity)
239  *out << "\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
240 }
241 
243  const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
244  Thyra::PreconditionerBase<double> *prec,
245  const Thyra::ESupportSolveUse supportSolveUse
246  ) const
247 {
248  using Teuchos::outArg;
249  using Teuchos::RCP;
250  using Teuchos::rcp;
251  using Teuchos::rcp_dynamic_cast;
252  using Teuchos::implicit_cast;
253 
254  Teuchos::Time totalTimer(""), timer("");
255  totalTimer.start(true);
256 
257  const RCP<Teuchos::FancyOStream> out = this->getOStream();
258  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
259 
260  bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
261 
262  Teuchos::OSTab tab(out);
263  if(mediumVerbosity)
264  *out << "\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
265 
266  Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
267 #ifdef _DEBUG
268  TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
269  TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
270 #endif
271 
272  //
273  // Unwrap and get the forward Epetra_Operator object
274  //
275  Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
276  Thyra::EOpTransp epetraFwdOpTransp;
277  Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
278  Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
279  double epetraFwdOpScalar;
280  epetraFwdOpViewExtractor_->getEpetraOpView(
281  fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
282  outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
283  );
284  // Get the concrete preconditioner object
285  StratimikosFactoryPreconditioner & defaultPrec
286  = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
287 
288  // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
289  RCP<Thyra::EpetraLinearOp> epetra_precOp
290  = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),true);
291 
292  // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
293  Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
294  if(epetra_precOp!=Teuchos::null)
295  teko_precOp = rcp_dynamic_cast<Teko::Epetra::InverseFactoryOperator>(epetra_precOp->epetra_op(),true);
296 
297  // Permform initialization if needed
298  const bool startingOver = (teko_precOp == Teuchos::null);
299  if(startingOver) {
300  // start over
301  invLib_ = Teuchos::null;
302  invFactory_ = Teuchos::null;
303  decomp_.clear();
304 
305  if(mediumVerbosity)
306  *out << "\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
307 
308  timer.start(true);
309 
310  std::stringstream ss;
311  ss << paramList_->get<std::string>("Strided Blocking");
312 
313  // figure out the decomposition requested by the string
314  while(not ss.eof()) {
315  int num=0;
316  ss >> num;
317 
318  TEUCHOS_ASSERT(num>0);
319 
320  decomp_.push_back(num);
321  }
322 
323  // build library, and set request handler (user defined!)
324  invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"),builder_);
325  invLib_->setRequestHandler(reqHandler_);
326 
327  // build preconditioner factory
328  invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
329 
330  // Create the initial preconditioner: DO NOT compute it yet
331  teko_precOp = rcp( new Teko::Epetra::InverseFactoryOperator(invFactory_));
332 
333  timer.stop();
334  if(mediumVerbosity)
335  Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
336  }
337 
338  if(mediumVerbosity)
339  *out << "\nComputing the preconditioner ...\n";
340 
341  // construct preconditioner
342  {
343  bool writeBlockOps = paramList_->get<bool>("Write Block Operator");
344  timer.start(true);
345  teko_precOp->initInverse();
346 
347  if(decomp_.size()==1) {
348  teko_precOp->rebuildInverseOperator(epetraFwdOp);
349 
350  // write out to disk
351  if(writeBlockOps) {
352  std::stringstream ss;
353  ss << "block-" << defaultPrec.getIter() << "_00.mm";
354  EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*rcp_dynamic_cast<const Epetra_CrsMatrix>(epetraFwdOp));
355  }
356  }
357  else {
358  Teuchos::RCP<Epetra_Operator> wrappedFwdOp
359  = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
360 
361  // write out to disk
362  if(writeBlockOps) {
363  std::stringstream ss;
364  ss << "block-" << defaultPrec.getIter();
365  // Teuchos::RCP<Teko::Epetra::StridedEpetraOperator> stridedJac
366  // = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedFwdOp);
367  Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac
368  = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedFwdOp);
369  if(stridedJac!=Teuchos::null) {
370  // write out blocks of strided operator
371  stridedJac->WriteBlocks(ss.str());
372  }
373  else TEUCHOS_ASSERT(false);
374  }
375 
376  teko_precOp->rebuildInverseOperator(wrappedFwdOp);
377  }
378 
379  timer.stop();
380  }
381 
382  if(mediumVerbosity)
383  Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
384 
385  // Initialize the output EpetraLinearOp
386  if(startingOver) {
387  epetra_precOp = rcp(new Thyra::EpetraLinearOp);
388  }
389  epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
390  ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
391  ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
392 
393  // Initialize the preconditioner
394  defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
395  defaultPrec.incrIter();
396  totalTimer.stop();
397 
398  if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
399  *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
400  if(mediumVerbosity)
401  *out << "\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
402 }
403 
404 
406  Thyra::PreconditionerBase<double> *prec,
407  Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > *fwdOp,
408  Thyra::ESupportSolveUse *supportSolveUse
409  ) const
410 {
411  TEUCHOS_TEST_FOR_EXCEPT(true);
412 }
413 
414 
415 // Overridden from ParameterListAcceptor
416 
417 
419  Teuchos::RCP<Teuchos::ParameterList> const& paramList
420  )
421 {
422  TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
423 
424  paramList->validateParametersAndSetDefaults(*this->getValidParameters(),0);
425  paramList_ = paramList;
426 
427 }
428 
429 
430 Teuchos::RCP<Teuchos::ParameterList>
432 {
433  return paramList_;
434 }
435 
436 
437 Teuchos::RCP<Teuchos::ParameterList>
439 {
440  Teuchos::RCP<ParameterList> _paramList = paramList_;
441  paramList_ = Teuchos::null;
442  return _paramList;
443 }
444 
445 
446 Teuchos::RCP<const Teuchos::ParameterList>
448 {
449  return paramList_;
450 }
451 
452 
453 Teuchos::RCP<const Teuchos::ParameterList>
455 {
456 
457  using Teuchos::rcp;
458  using Teuchos::tuple;
459  using Teuchos::implicit_cast;
460  using Teuchos::rcp_implicit_cast;
461 
462  static RCP<const ParameterList> validPL;
463 
464  if(is_null(validPL)) {
465 
466  Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
467 
468  pl->set("Test Block Operator",false,
469  "If Stratiikos/Teko is used to break an operator into its parts,\n"
470  "then setting this parameter to true will compare applications of the\n"
471  "segregated operator to the original operator.");
472  pl->set("Write Block Operator",false,
473  "Write out the segregated operator to disk with the name \"block-?_xx\"");
474  pl->set("Strided Blocking","1",
475  "Assuming that the user wants Strided blocking, break the operating into\n"
476  "blocks. The syntax can be thought to be associated with the solution\n"
477  "vector. For example if your variables are [u v w p T], and we want [u v w]\n"
478  "blocked together, and p and T seperate then the relevant string is \"3 1 1\".\n"
479  "Meaning put the first 3 unknowns per node together and sperate the v and w\n"
480  "components.");
481  pl->set("Reorder Type","",
482  "This specifies how the blocks are reordered for use in the preconditioner.\n"
483  "For example, assume the linear system is generated from 3D Navier-Stokes\n"
484  "with an energy equation, yielding the unknowns [u v w p T]. If the\n"
485  "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n"
486  "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n"
487  "velocity and pressure forming an inner two-by-two block, and then the\n"
488  "temperature unknowns forming a two-by-two system with the velocity-pressure\n"
489  "block.");
490  pl->set("Inverse Type","Amesos",
491  "The type of inverse operator the user wants. This can be one of the defaults\n"
492  "from Stratimikos, or a Teko preconditioner defined in the\n"
493  "\"Inverse Factory Library\".");
494  pl->sublist("Inverse Factory Library",false,"Definition of Teko preconditioners.");
495 
496  validPL = pl;
497  }
498 
499  return validPL;
500 }
501 
505 Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
506  const Teuchos::RCP<const Epetra_Operator> & Jac,
507  const Teuchos::RCP<Epetra_Operator> & wrapInput,
508  std::ostream & out
509  ) const
510 {
511  Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
512 
513 // // initialize jacobian
514 // if(wrappedOp==Teuchos::null)
515 // {
516 // wrappedOp = Teuchos::rcp(new Teko::Epetra::StridedEpetraOperator(decomp_,Jac));
517 //
518 // // reorder the blocks if requested
519 // std::string reorderType = paramList_->get<std::string>("Reorder Type");
520 // if(reorderType!="") {
521 // RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
522 //
523 // // out << "Teko: Reordering = " << brm->toString() << std::endl;
524 // Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->Reorder(*brm);
525 // }
526 // }
527 // else {
528 // Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->RebuildOps();
529 // }
530 //
531 // // test blocked operator for correctness
532 // if(paramList_->get<bool>("Test Block Operator")) {
533 // bool result
534 // = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
535 //
536 // out << "Teko: Tested operator correctness: " << (result ? "passed" : "FAILED!") << std::endl;
537 // }
538 
539  // initialize jacobian
540  if(wrappedOp==Teuchos::null)
541  {
542  // build strided vector
543  std::vector<std::vector<int> > vars;
544  buildStridedVectors(*Jac,decomp_,vars);
545  wrappedOp = Teuchos::rcp(new Teko::Epetra::BlockedEpetraOperator(vars,Jac));
546 
547  // reorder the blocks if requested
548  std::string reorderType = paramList_->get<std::string>("Reorder Type");
549  if(reorderType!="") {
550  RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
551 
552  // out << "Teko: Reordering = " << brm->toString() << std::endl;
553  Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->Reorder(*brm);
554  }
555  }
556  else {
557  Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->RebuildOps();
558  }
559 
560  // test blocked operator for correctness
561  if(paramList_->get<bool>("Test Block Operator")) {
562  bool result
563  = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
564 
565  out << "Teko: Tested operator correctness: " << (result ? "passed" : "FAILED!") << std::endl;
566  }
567 
568  return wrappedOp;
569 }
570 
572 {
573  std::ostringstream oss;
574  oss << "Teko::StratimikosFactory";
575  return oss.str();
576 }
577 
578 void StratimikosFactory::buildStridedVectors(const Epetra_Operator & Jac,
579  const std::vector<int> & decomp,
580  std::vector<std::vector<int> > & vars) const
581 {
582  const Epetra_Map & rangeMap = Jac.OperatorRangeMap();
583 
584  // compute total number of variables
585  int numVars = 0;
586  for(std::size_t i=0;i<decomp.size();i++)
587  numVars += decomp[i];
588 
589  // verify that the decomposition is appropriate for this matrix
590  TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars)==0);
591  TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars)==0);
592 
593  int * globalIds = rangeMap.MyGlobalElements();
594 
595  vars.resize(decomp.size());
596  for(int i=0;i<rangeMap.NumMyElements();) {
597 
598  // for each "node" copy global ids to vectors
599  for(std::size_t d=0;d<decomp.size();d++) {
600  // for this variable copy global ids to variable arrays
601  int current = decomp[d];
602  for(int v=0;v<current;v++,i++)
603  vars[d].push_back(globalIds[i]);
604  }
605  }
606 }
607 
608 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
609  const std::string & stratName)
610 {
611  TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
612  "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
613 
614  Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
615  = Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder(builder));
616 
617  // use default constructor to add Teko::StratimikosFactory
618  Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(builderCopy,Teuchos::null));
619  builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
620 }
621 
622 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
623  const Teuchos::RCP<Teko::RequestHandler> & rh,
624  const std::string & stratName)
625 {
626  TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
627  "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
628 
629  Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
630  = Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder(builder));
631 
632  // build an instance of a Teuchos::AbsractFactory<Thyra::PFB> so request handler is passed onto
633  // the resulting StratimikosFactory
634  Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(builderCopy,rh));
635  builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
636 }
637 
638 } // namespace Thyra
bool applySupportsConj(Thyra::EConj conj) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void setRequestHandler(const Teuchos::RCP< Teko::RequestHandler > &rh)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void initializePrec(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
void uninitializePrec(Thyra::PreconditionerBase< double > *prec, Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > *fwdOp, Thyra::ESupportSolveUse *supportSolveUse) const
bool applyTransposeSupportsConj(Thyra::EConj conj) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
void initializePrec_Thyra(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
A single Epetra wrapper for all operators constructed from an inverse operator.
bool isCompatible(const Thyra::LinearOpSourceBase< double > &fwdOp) const
Teuchos::RCP< Thyra::PreconditionerBase< double > > createPrec() const
This class takes a blocked linear op and represents it in a flattened form.
void initializePrec_Epetra(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Tear about a user specified Epetra_Operator (CrsMatrix) using a vector of vectors of GIDs for each bl...
virtual void WriteBlocks(const std::string &prefix) const
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.