Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_LinearOpWithSolveFactoryExamples.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP
00030 #define THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP
00031 
00032 
00033 #include "Thyra_LinearOpWithSolveFactoryBase.hpp"
00034 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00035 #include "Thyra_LinearOpWithSolveBase.hpp"
00036 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00037 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00038 #include "Thyra_DefaultPreconditioner.hpp"
00039 #include "Thyra_MultiVectorStdOps.hpp"
00040 #include "Thyra_VectorStdOps.hpp"
00041 #include "Thyra_VectorBase.hpp"
00042 
00043 
00044 namespace Thyra {
00045 
00046 
00047 //
00048 // Helper code
00049 //
00050 
00061 template<class Scalar>
00062 class LinearOpChanger {
00063 public:
00065   virtual ~LinearOpChanger() {}
00067   virtual void changeOp( const Teuchos::Ptr<LinearOpBase<Scalar> > &op ) const = 0;
00068 };
00069 
00073 template<class Scalar>
00074 class NullLinearOpChanger : public LinearOpChanger<Scalar> {
00075 public:
00077   void changeOp( const Teuchos::Ptr<LinearOpBase<Scalar> > &op ) const {}
00078 };
00079 
00080 
00081 } // namespace Thyra
00082 
00083 
00084 //
00085 // Individual non-externally preconditioned use cases
00086 //
00087 
00088 
00089 // begin singleLinearSolve
00093 template<class Scalar>
00094 void singleLinearSolve(
00095   const Thyra::LinearOpBase<Scalar> &A,
00096   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00097   const Thyra::VectorBase<Scalar> &b,
00098   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x,
00099   Teuchos::FancyOStream &out
00100   )
00101 {
00102   typedef Teuchos::ScalarTraits<Scalar> ST;
00103   using Teuchos::rcpFromRef;
00104   Teuchos::OSTab tab(out);
00105   out << "\nPerforming a single linear solve ...\n";
00106   // Create the LOWSB object that will be used to solve the linear system
00107   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00108     Thyra::linearOpWithSolve(lowsFactory, rcpFromRef(A));
00109   // Solve the system using a default solve criteria using a non-member helper function 
00110   assign(x, ST::zero()); // Must initialize to a guess before solve!
00111   Thyra::SolveStatus<Scalar> 
00112     status = Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b, x);
00113   out << "\nSolve status:\n" << status;
00114 } // end singleLinearSolve
00115 
00116 
00117 // begin createScaledAdjointLinearOpWithSolve
00122 template<class Scalar>
00123 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00124 createScaledAdjointLinearOpWithSolve(
00125   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00126   const Scalar &scalar,
00127   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00128   Teuchos::FancyOStream &out
00129   )
00130 {
00131   Teuchos::OSTab tab(out);
00132   out << "\nCreating a scaled adjoint LinearOpWithSolveBase object ...\n";
00133   // Create the LOWSB object that will be used to solve the linear system
00134   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleAdjointA =
00135     Thyra::linearOpWithSolve(lowsFactory,scale(scalar,adjoint(A)));
00136   out << "\nCreated LOWSB object:\n" << describe(*invertibleAdjointA,
00137     Teuchos::VERB_MEDIUM);
00138   return invertibleAdjointA;
00139 } // end createScaledAdjointLinearOpWithSolve
00140 
00141 
00142 // begin solveNumericalChangeSolve
00147 template<class Scalar>
00148 void solveNumericalChangeSolve(
00149   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > &A,
00150   const Thyra::LinearOpChanger<Scalar> &opChanger,
00151   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00152   const Thyra::VectorBase<Scalar> &b1,
00153   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x1,
00154   const Thyra::VectorBase<Scalar> &b2,
00155   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x2,
00156   Teuchos::FancyOStream &out
00157   )
00158 {
00159   using Teuchos::as; using Teuchos::ptr; using Teuchos::rcpFromPtr;
00160   Teuchos::OSTab tab(out);
00161   out << "\nPerforming a solve, changing the operator, then performing another"
00162       << " solve ...\n";
00163   // Get a local non-owned RCP to A to be used by lowsFactory
00164   Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > rcpA = rcpFromPtr(A);
00165   // Create the LOWSB object that will be used to solve the linear system
00166   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00167     lowsFactory.createOp();
00168   // Initialize the invertible linear operator given the forward operator
00169   Thyra::initializeOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00170   // Solve the system using a default solve criteria using a non-member helper function
00171   Thyra::assign(x1, as<Scalar>(0.0));
00172   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b1, x1);
00173   // Before the forward operator A is changed it is recommended that you
00174   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00175   // while it may be in an inconsistent state from the time between *A changes
00176   // and *invertibleA is explicitly updated. However, this step is not
00177   // required!
00178   Thyra::uninitializeOp<Scalar>(lowsFactory, invertibleA.ptr());
00179   // Change the operator and reinitialize the invertible operator
00180   opChanger.changeOp(A);
00181   Thyra::initializeOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00182   // Note that above will reuse any factorization structures that may have been
00183   // created in the first call to initializeOp(...).
00184   // Finally, solve another linear system with new values of A
00185   Thyra::assign<Scalar>(x2, as<Scalar>(0.0));
00186   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b2, x2);
00187 } // end solveNumericalChangeSolve
00188 
00189 
00190 // begin solveSmallNumericalChangeSolve
00196 template<class Scalar>
00197 void solveSmallNumericalChangeSolve(
00198   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > &A,
00199   const Thyra::LinearOpChanger<Scalar> &opSmallChanger,
00200   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00201   const Thyra::VectorBase<Scalar> &b1,
00202   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x1,
00203   const Thyra::VectorBase<Scalar> &b2,
00204   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x2,
00205   Teuchos::FancyOStream &out
00206   )
00207 {
00208   using Teuchos::ptr; using Teuchos::as; using Teuchos::rcpFromPtr;
00209   Teuchos::OSTab tab(out);
00210   out << "\nPerforming a solve, changing the operator in a very small way,"
00211       << " then performing another solve ...\n";
00212   // Get a local non-owned RCP to A to be used by lowsFactory
00213   Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > rcpA = rcpFromPtr(A);
00214   // Create the LOWSB object that will be used to solve the linear system
00215   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00216     lowsFactory.createOp();
00217   // Initialize the invertible linear operator given the forward operator
00218   Thyra::initializeOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00219   // Solve the system using a default solve criteria using a non-member helper function
00220   Thyra::assign(x1, as<Scalar>(0.0));
00221   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b1, x1);
00222   // Before the forward operator A is changed it is recommended that you
00223   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00224   // while it may be in an inconsistent state from the time between *A changes
00225   // and *invertibleA is explicitly updated. However, this step is not
00226   // required!
00227   Thyra::uninitializeOp<Scalar>(lowsFactory, invertibleA.ptr());
00228   // Change the operator and reinitialize the invertible operator
00229   opSmallChanger.changeOp(A);
00230   Thyra::initializeAndReuseOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00231   // Note that above a maximum amount of reuse will be achieved, such as
00232   // keeping the same precondtioner.
00233   Thyra::assign(x2, as<Scalar>(0.0));
00234   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b2, x2);
00235 } // end solveSmallNumericalChangeSolve
00236 
00237 
00238 // begin solveMajorChangeSolve
00244 template<class Scalar>
00245 void solveMajorChangeSolve(
00246   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > &A,
00247   const Thyra::LinearOpChanger<Scalar> &opMajorChanger,
00248   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00249   const Thyra::VectorBase<Scalar> &b1,
00250   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x1,
00251   const Thyra::VectorBase<Scalar> &b2,
00252   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x2,
00253   Teuchos::FancyOStream &out
00254   )
00255 {
00256   using Teuchos::as; using Teuchos::rcpFromPtr;
00257   Teuchos::OSTab tab(out);
00258   out << "\nPerforming a solve, changing the operator in a major way, then performing"
00259       << " another solve ...\n";
00260   // Get a local non-owned RCP to A to be used by lowsFactory
00261   Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > rcpA = rcpFromPtr(A);
00262   // Create the LOWSB object that will be used to solve the linear system
00263   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00264     lowsFactory.createOp();
00265   // Initialize the invertible linear operator given the forward operator
00266   Thyra::initializeOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00267   // Solve the system using a default solve criteria using a non-member helper function
00268   Thyra::assign(x1, as<Scalar>(0.0));
00269   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b1, x1);
00270   // Before the forward operator A is changed it is recommended that you
00271   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00272   // while it may be in an inconsistent state from the time between *A changes
00273   // and *invertibleA is explicitly updated. However, this step is not
00274   // required!
00275   Thyra::uninitializeOp<Scalar>(lowsFactory, invertibleA.ptr());
00276   // Change the operator in some major way (perhaps even changing its structure)
00277   opMajorChanger.changeOp(A);
00278   // Recreate the LOWSB object and initialize it from scratch
00279   invertibleA = lowsFactory.createOp();
00280   Thyra::initializeOp<Scalar>(lowsFactory, rcpA, invertibleA.ptr());
00281   // Solve another set of linear systems
00282   Thyra::assign(x2, as<Scalar>(0.0));
00283   Thyra::solve<Scalar>(*invertibleA, Thyra::NOTRANS, b2, x2);
00284 } // end solveMajorChangeSolve
00285 
00286 
00287 //
00288 // Individual externally preconditioned use cases
00289 //
00290 
00291 
00292 // begin createGeneralPreconditionedLinearOpWithSolve
00298 template<class Scalar>
00299 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00300 createGeneralPreconditionedLinearOpWithSolve(
00301   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00302   const Teuchos::RCP<const Thyra::PreconditionerBase<Scalar> > &P,
00303   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00304   Teuchos::FancyOStream &out
00305   )
00306 {
00307   Teuchos::OSTab tab(out);
00308   out << "\nCreating an externally preconditioned LinearOpWithSolveBase object ...\n";
00309   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00310     lowsFactory.createOp();
00311   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A, P, invertibleA.ptr());
00312   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00313   return invertibleA;
00314 } // end createGeneralPreconditionedLinearOpWithSolve
00315 
00316 
00317 // begin createUnspecifiedPreconditionedLinearOpWithSolve
00323 template<class Scalar>
00324 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00325 createUnspecifiedPreconditionedLinearOpWithSolve(
00326   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00327   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op,
00328   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00329   Teuchos::FancyOStream &out
00330   )
00331 {
00332   Teuchos::OSTab tab(out);
00333   out << "\nCreating an LinearOpWithSolveBase object given a preconditioner operator"
00334       << " not targeted to the left or right ...\n";
00335   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00336     lowsFactory.createOp();
00337   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A,
00338     Thyra::unspecifiedPrec<Scalar>(P_op), invertibleA.ptr());
00339   // Above, the lowsFactory object will decide whether to apply the single
00340   // preconditioner operator on the left or on the right.
00341   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00342   return invertibleA;
00343 } // end createUnspecifiedPreconditionedLinearOpWithSolve
00344 
00345 
00346 // begin createLeftPreconditionedLinearOpWithSolve
00352 template<class Scalar>
00353 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00354 createLeftPreconditionedLinearOpWithSolve(
00355   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00356   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_left,
00357   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00358   Teuchos::FancyOStream &out
00359   )
00360 {
00361   Teuchos::OSTab tab(out);
00362   out << "\nCreating an LinearOpWithSolveBase object given a left preconditioner"
00363       << " operator ...\n";
00364   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00365     lowsFactory.createOp();
00366   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A,
00367     Thyra::leftPrec<Scalar>(P_op_left), invertibleA.ptr());
00368   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00369   return invertibleA;
00370 } // end createLeftPreconditionedLinearOpWithSolve
00371 
00372 
00373 // begin createRightPreconditionedLinearOpWithSolve
00379 template<class Scalar>
00380 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00381 createRightPreconditionedLinearOpWithSolve(
00382   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00383   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_right,
00384   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00385   Teuchos::FancyOStream &out
00386   )
00387 {
00388   Teuchos::OSTab tab(out);
00389   out << "\nCreating an LinearOpWithSolveBase object given a right"
00390       << " preconditioner operator ...\n";
00391   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00392     lowsFactory.createOp();
00393   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A,
00394     Thyra::rightPrec<Scalar>(P_op_right), invertibleA.ptr());
00395   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00396   return invertibleA;
00397 } // end createRightPreconditionedLinearOpWithSolve
00398 
00399 
00400 // begin createLeftRightPreconditionedLinearOpWithSolve
00406 template<class Scalar>
00407 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00408 createLeftRightPreconditionedLinearOpWithSolve(
00409   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00410   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_left,
00411   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_right,
00412   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00413   Teuchos::FancyOStream &out
00414   )
00415 {
00416   Teuchos::OSTab tab(out);
00417   out << "\nCreating an LinearOpWithSolveBase object given a left and"
00418       << "right preconditioner operator ...\n";
00419   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00420     lowsFactory.createOp();
00421   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A,
00422     Thyra::splitPrec<Scalar>(P_op_left, P_op_right), invertibleA.ptr());
00423   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00424   return invertibleA;
00425 } // end createLeftRightPreconditionedLinearOpWithSolve
00426 
00427 
00428 // begin createMatrixPreconditionedLinearOpWithSolve
00434 template<class Scalar>
00435 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00436 createMatrixPreconditionedLinearOpWithSolve(
00437   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
00438   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A_approx,
00439   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00440   Teuchos::FancyOStream &out
00441   )
00442 {
00443   Teuchos::OSTab tab(out);
00444   out << "\nCreating a LinearOpWithSolveBase object given an approximate forward"
00445       << " operator to define the preconditioner ...\n";
00446   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00447     lowsFactory.createOp();
00448   Thyra::initializeApproxPreconditionedOp<Scalar>(lowsFactory, A, A_approx,
00449     invertibleA.ptr());
00450   out << "\nCreated LOWSB object:\n" << describe(*invertibleA, Teuchos::VERB_MEDIUM);
00451   return invertibleA;
00452 } // end createMatrixPreconditionedLinearOpWithSolve
00453 
00454 
00455 // begin externalPreconditionerReuseWithSolves
00459 template<class Scalar>
00460 void externalPreconditionerReuseWithSolves(
00461   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > &A_inout,
00462   const Thyra::LinearOpChanger<Scalar> &opChanger,
00463   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00464   const Thyra::PreconditionerFactoryBase<Scalar> &precFactory,
00465   const Thyra::VectorBase<Scalar> &b1,
00466   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x1,
00467   const Thyra::VectorBase<Scalar> &b2,
00468   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x2,
00469   Teuchos::FancyOStream &out
00470   )
00471 {
00472   using Teuchos::tab; using Teuchos::rcpFromPtr;
00473   typedef Teuchos::ScalarTraits<Scalar> ST;
00474   Teuchos::OSTab tab2(out);
00475   out << "\nShowing resuse of the preconditioner ...\n";
00476   Teuchos::RCP<Thyra::LinearOpBase<Scalar> > A = rcpFromPtr(A_inout);
00477   // Create the initial preconditioner for the input forward operator
00478   Teuchos::RCP<Thyra::PreconditionerBase<Scalar> > P =
00479     precFactory.createPrec();
00480   Thyra::initializePrec<Scalar>(precFactory, A, P.ptr());
00481   // Create the invertible LOWS object given the preconditioner
00482   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > invertibleA =
00483     lowsFactory.createOp();
00484   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A, P, invertibleA.ptr());
00485   // Solve the first linear system
00486   assign(x1, ST::zero());
00487   Thyra::SolveStatus<Scalar> status1 = Thyra::solve<Scalar>(*invertibleA,
00488     Thyra::NOTRANS, b1, x1);
00489   out << "\nSolve status:\n" << status1;
00490   // Change the forward linear operator without changing the preconditioner
00491   opChanger.changeOp(A.ptr());
00492   // Warning! After the above change the integrity of the preconditioner
00493   // linear operators in P is undefined. For some implementations of the
00494   // preconditioner, its behavior will remain unchanged (e.g. ILU) which in
00495   // other cases the behavior will change but the preconditioner will still
00496   // work (e.g. Jacobi). However, there may be valid implementations where
00497   // the preconditioner will simply break if the forward operator that it is
00498   // based on breaks.
00499   //
00500   // Reinitialize the LOWS object given the updated forward operator A and the
00501   // old preconditioner P.
00502   Thyra::initializePreconditionedOp<Scalar>(lowsFactory, A, P, invertibleA.ptr());
00503   // Solve the second linear system
00504   assign(x2, ST::zero());
00505   Thyra::SolveStatus<Scalar>status2 = Thyra::solve<Scalar>(*invertibleA,
00506     Thyra::NOTRANS, b2, x2);
00507   out << "\nSolve status:\n" << status2;
00508 } // end externalPreconditionerReuseWithSolves
00509 
00510 
00511 //
00512 // Combined use cases
00513 //
00514 
00515 
00521 template<class Scalar>
00522 void nonExternallyPreconditionedLinearSolveUseCases(
00523   const Thyra::LinearOpBase<Scalar> &A,
00524   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00525   bool supportsAdjoints,
00526   Teuchos::FancyOStream &out
00527   )
00528 {
00529   using Teuchos::as;
00530   Teuchos::OSTab tab(out);
00531   out << "\nRunning example use cases for a LinearOpWithSolveFactoryBase object ...\n";
00532   // Create a non-const A object (don't worry, it will not be changed)
00533   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > A_nonconst =
00534     Teuchos::ptrFromRef(const_cast<Thyra::LinearOpBase<Scalar>&>(A));
00535   // Create the RHS (which is just a random set of coefficients)
00536   Teuchos::RCP<Thyra::VectorBase<Scalar> >
00537     b1 = Thyra::createMember(A.range()),
00538     b2 = Thyra::createMember(A.range());
00539   Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), b1.ptr() );
00540   Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), b2.ptr() );
00541   // Create the LHS for the linear solve
00542   Teuchos::RCP<Thyra::VectorBase<Scalar> >
00543     x1 = Thyra::createMember(A.domain()),
00544     x2 = Thyra::createMember(A.domain());
00545   // Perform a single, non-adjoint, linear solve
00546   singleLinearSolve(A, lowsFactory, *b1, x1.ptr(), out);
00547   // Creating a scaled adjoint LinearOpWithSolveBase object
00548   if(supportsAdjoints) {
00549     Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00550       invertibleAdjointA = createScaledAdjointLinearOpWithSolve(
00551         Teuchos::rcp(&A,false),as<Scalar>(2.0),lowsFactory,out);
00552   }
00553   // Perform a solve, change the operator, and then solve again.
00554   solveNumericalChangeSolve<Scalar>(
00555     A_nonconst,                            // Don't worry, it will not be changed!
00556     Thyra::NullLinearOpChanger<Scalar>(),  // This object will not really change A!
00557     lowsFactory, *b1, x1.ptr(), *b2, x1.ptr(), out );
00558   // Perform a solve, change the operator in a very small way, and then solve again.
00559   solveSmallNumericalChangeSolve<Scalar>(
00560     A_nonconst,                            // Don't worry, it will not be changed!
00561     Thyra::NullLinearOpChanger<Scalar>(),  // This object will not really change A!
00562     lowsFactory, *b1, x1.ptr(), *b2, x1.ptr(), out );
00563   // Perform a solve, change the operator in a major way, and then solve again.
00564   solveMajorChangeSolve<Scalar>(
00565     A_nonconst,                            // Don't worry, it will not be changed!
00566     Thyra::NullLinearOpChanger<Scalar>(),  // This object will not really change A!
00567     lowsFactory, *b1, x1.ptr(), *b2, x1.ptr(), out );
00568 }
00569 
00570 
00576 template<class Scalar>
00577 void externallyPreconditionedLinearSolveUseCases(
00578   const Thyra::LinearOpBase<Scalar> &A,
00579   const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
00580   const Thyra::PreconditionerFactoryBase<Scalar> &precFactory,
00581   const bool supportsLeftPrec,
00582   const bool supportsRightPrec,
00583   Teuchos::FancyOStream &out
00584   )
00585 {
00586   using Teuchos::rcpFromRef; using Teuchos::as;
00587   Teuchos::OSTab tab(out);
00588   out << "\nRunning example use cases with an externally defined"
00589       << " preconditioner with a LinearOpWithSolveFactoryBase object ...\n";
00590   // Create a non-const A object (don't worry, it will not be changed)
00591   const Teuchos::Ptr<Thyra::LinearOpBase<Scalar> > A_nonconst =
00592     Teuchos::ptrFromRef(const_cast<Thyra::LinearOpBase<Scalar>&>(A));
00593   // Create the RHS (which is just a random set of coefficients)
00594   Teuchos::RCP<Thyra::VectorBase<Scalar> >
00595     b1 = Thyra::createMember(A.range()),
00596     b2 = Thyra::createMember(A.range());
00597   Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), b1.ptr() );
00598   Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), b2.ptr() );
00599   // Create the LHS for the linear solve
00600   Teuchos::RCP<Thyra::VectorBase<Scalar> >
00601     x1 = Thyra::createMember(A.domain()),
00602     x2 = Thyra::createMember(A.domain());
00603   // Create a preconditioner for the input forward operator
00604   Teuchos::RCP<Thyra::PreconditionerBase<Scalar> >
00605     P = precFactory.createPrec();
00606   Thyra::initializePrec<Scalar>(precFactory, rcpFromRef(A), P.ptr());
00607   // Above, we don't really know the nature of the preconditioner. It could a
00608   // single linear operator to be applied on the left or the right or it could
00609   // be a split preconditioner with different linear operators to be applied
00610   // on the right or left. Or, it could be a single linear operator that is
00611   // not targeted to the left or the right.
00612   //
00613   // Create a LOWSB object given the created preconditioner
00614   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00615     invertibleA = createGeneralPreconditionedLinearOpWithSolve<Scalar>(
00616       Teuchos::rcp(&A, false), P.getConst(), lowsFactory, out);
00617   // Grab a preconditioner operator out of the preconditioner object
00618   Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > P_op;
00619   if (nonnull(P_op=P->getUnspecifiedPrecOp()));
00620   else if (nonnull(P_op=P->getLeftPrecOp()));
00621   else if (nonnull(P_op=P->getRightPrecOp()));
00622   // Create a LOWSB object given an unspecified preconditioner operator
00623   invertibleA = createUnspecifiedPreconditionedLinearOpWithSolve(
00624     rcpFromRef(A), P_op, lowsFactory, out);
00625   // Create a LOWSB object given a left preconditioner operator
00626   if(supportsLeftPrec) {
00627     invertibleA = createLeftPreconditionedLinearOpWithSolve(
00628       rcpFromRef(A), P_op, lowsFactory,out);
00629   }
00630   // Create a LOWSB object given a right preconditioner operator
00631   if(supportsRightPrec) {
00632     invertibleA = createRightPreconditionedLinearOpWithSolve(
00633       rcpFromRef(A), P_op, lowsFactory, out);
00634   }
00635   // Create a LOWSB object given (bad set of) left and right preconditioner
00636   // operators
00637   if( supportsLeftPrec && supportsRightPrec ) {
00638     invertibleA = createLeftRightPreconditionedLinearOpWithSolve(
00639       rcpFromRef(A), P_op, P_op, lowsFactory, out);
00640   }
00641   // Create a LOWSB object given a (very good) approximate forward linear
00642   // operator to construct the preconditoner from..
00643   invertibleA = createMatrixPreconditionedLinearOpWithSolve<Scalar>(
00644     rcpFromRef(A), rcpFromRef(A), lowsFactory,out);
00645   // Preconditioner reuse example
00646   externalPreconditionerReuseWithSolves<Scalar>(
00647     A_nonconst,                            // Don't worry, it will not be changed!
00648     Thyra::NullLinearOpChanger<Scalar>(),  // This object will not really change A!
00649     lowsFactory, precFactory,
00650     *b1, x1.ptr(), *b2, x2.ptr(), out );
00651 }
00652 
00653 
00654 #endif // THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines