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 #include "Thyra_LinearOpWithSolveFactoryBase.hpp"
00033 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00034 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00035 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00036 #include "Thyra_DefaultPreconditioner.hpp"
00037 
00038 namespace Thyra {
00039 
00040 //
00041 // Helper code
00042 //
00043 
00054 template<class Scalar>
00055 class LinearOpChanger {
00056 public:
00058   virtual ~LinearOpChanger() {}
00060   virtual void changeOp( LinearOpBase<Scalar> *op ) const = 0;
00061 };
00062 
00066 template<class Scalar>
00067 class NullLinearOpChanger : public LinearOpChanger<Scalar> {
00068 public:
00070   void changeOp( LinearOpBase<Scalar> *op ) const {}
00071 };
00072 
00073 } // namespace Thyra
00074 
00075 //
00076 // Individual non-externally preconditioned use cases
00077 //
00078 
00079 // begin singleLinearSolve
00083 template<class Scalar>
00084 void singleLinearSolve(
00085   const Thyra::LinearOpBase<Scalar>                    &A
00086   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00087   ,const Thyra::VectorBase<Scalar>                     &b
00088   ,Thyra::VectorBase<Scalar>                           *x
00089   ,Teuchos::FancyOStream                               &out
00090   )
00091 {
00092   Teuchos::OSTab tab(out);
00093   out << "\nPerforming a single linear solve ...\n";
00094   // Create the LOWSB object that will be used to solve the linear system
00095   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00096     invertibleA = Thyra::linearOpWithSolve(lowsFactory,Teuchos::rcp(&A,false));
00097   // Solve the system using a default solve criteria using a non-member helper function 
00098   assign(&*x,Teuchos::ScalarTraits<Scalar>::zero()); // Must initialize to a guess before solve!
00099   Thyra::SolveStatus<Scalar>
00100     status = Thyra::solve(*invertibleA,Thyra::NOTRANS,b,x);
00101   out << "\nSolve status:\n" << status;
00102 } // end singleLinearSolve
00103 
00104 // begin createScaledAdjointLinearOpWithSolve
00109 template<class Scalar>
00110 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00111 createScaledAdjointLinearOpWithSolve(
00112   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >   &A
00113   ,const Scalar                                                    &scalar
00114   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>               &lowsFactory
00115   ,Teuchos::FancyOStream                                           &out
00116   )
00117 {
00118   Teuchos::OSTab tab(out);
00119   out << "\nCreating a scaled adjoint LinearOpWithSolveBase object ...\n";
00120   // Create the LOWSB object that will be used to solve the linear system
00121   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00122     invertibleAdjointA
00123     = Thyra::linearOpWithSolve(lowsFactory,scale(scalar,adjoint(A)));
00124   out << "\nCreated LOWSB object:\n" << describe(*invertibleAdjointA,Teuchos::VERB_MEDIUM);
00125   return invertibleAdjointA;
00126 } // end createScaledAdjointLinearOpWithSolve
00127 
00128 // begin solveNumericalChangeSolve
00133 template<class Scalar>
00134 void solveNumericalChangeSolve(
00135   Thyra::LinearOpBase<Scalar>                          *A
00136   ,const Thyra::LinearOpChanger<Scalar>                &opChanger
00137   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00138   ,const Thyra::VectorBase<Scalar>                     &b1
00139   ,Thyra::VectorBase<Scalar>                           *x1
00140   ,const Thyra::VectorBase<Scalar>                     &b2
00141   ,Thyra::VectorBase<Scalar>                           *x2
00142   ,Teuchos::FancyOStream                               &out
00143   )
00144 {
00145   Teuchos::OSTab tab(out);
00146   out << "\nPerforming a solve, changing the operator, then performing another solve ...\n";
00147   // Get a local non-owned RCP to A to be used by lowsFactory
00148   Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00149     rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
00150   // Create the LOWSB object that will be used to solve the linear system
00151   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00152     invertibleA = lowsFactory.createOp();
00153   // Initialize the invertible linear operator given the forward operator
00154   Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00155   // Solve the system using a default solve criteria using a non-member helper function
00156   Thyra::assign(x1,Scalar(0.0));
00157   solve(*invertibleA,Thyra::NOTRANS,b1,x1);
00158   // Before the forward operator A is changed it is recommended that you
00159   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00160   // while it may be in an inconsistent state from the time between *A changes
00161   // and *invertibleA is explicitly updated.  However, this step is not
00162   // required!
00163   Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
00164   // Change the operator and reinitialize the invertible operator
00165   opChanger.changeOp(A);
00166   Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00167   // Note that above will reuse any factorization structures that may have been
00168   // created in the first call to initializeOp(...).
00169   // Finally, solve another linear system with new values of A
00170   Thyra::assign(x2,Scalar(0.0));
00171   solve(*invertibleA,Thyra::NOTRANS,b2,x2);
00172 } // end solveNumericalChangeSolve
00173 
00174 // begin solveSmallNumericalChangeSolve
00180 template<class Scalar>
00181 void solveSmallNumericalChangeSolve(
00182   Thyra::LinearOpBase<Scalar>                          *A
00183   ,const Thyra::LinearOpChanger<Scalar>                &opSmallChanger
00184   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00185   ,const Thyra::VectorBase<Scalar>                     &b1
00186   ,Thyra::VectorBase<Scalar>                           *x1
00187   ,const Thyra::VectorBase<Scalar>                     &b2
00188   ,Thyra::VectorBase<Scalar>                           *x2
00189   ,Teuchos::FancyOStream                               &out
00190   )
00191 {
00192   Teuchos::OSTab tab(out);
00193   out << "\nPerforming a solve, changing the operator in a very small way, then performing another solve ...\n";
00194   // Get a local non-owned RCP to A to be used by lowsFactory
00195   Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00196     rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
00197   // Create the LOWSB object that will be used to solve the linear system
00198   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00199     invertibleA = lowsFactory.createOp();
00200   // Initialize the invertible linear operator given the forward operator
00201   Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00202   // Solve the system using a default solve criteria using a non-member helper function
00203   Thyra::assign(x1,Scalar(0.0));
00204   solve(*invertibleA,Thyra::NOTRANS,b1,x1);
00205   // Before the forward operator A is changed it is recommended that you
00206   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00207   // while it may be in an inconsistent state from the time between *A changes
00208   // and *invertibleA is explicitly updated.  However, this step is not
00209   // required!
00210   Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
00211   // Change the operator and reinitialize the invertible operator
00212   opSmallChanger.changeOp(A);
00213   Thyra::initializeAndReuseOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00214   // Note that above a maximum amount of reuse will be achieved, such as
00215   // keeping the same precondtioner.
00216   Thyra::assign(x2,Scalar(0.0));
00217   solve(*invertibleA,Thyra::NOTRANS,b2,x2);
00218 } // end solveSmallNumericalChangeSolve
00219 
00220 // begin solveMajorChangeSolve
00226 template<class Scalar>
00227 void solveMajorChangeSolve(
00228   Thyra::LinearOpBase<Scalar>                          *A
00229   ,const Thyra::LinearOpChanger<Scalar>                &opMajorChanger
00230   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00231   ,const Thyra::VectorBase<Scalar>                     &b1
00232   ,Thyra::VectorBase<Scalar>                           *x1
00233   ,const Thyra::VectorBase<Scalar>                     &b2
00234   ,Thyra::VectorBase<Scalar>                           *x2
00235   ,Teuchos::FancyOStream                               &out
00236   )
00237 {
00238   Teuchos::OSTab tab(out);
00239   out << "\nPerforming a solve, changing the operator in a major way, then performing another solve ...\n";
00240   // Get a local non-owned RCP to A to be used by lowsFactory
00241   Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00242     rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
00243   // Create the LOWSB object that will be used to solve the linear system
00244   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00245     invertibleA = lowsFactory.createOp();
00246   // Initialize the invertible linear operator given the forward operator
00247   Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00248   // Solve the system using a default solve criteria using a non-member helper function
00249   Thyra::assign(x1,Scalar(0.0));
00250   solve(*invertibleA,Thyra::NOTRANS,b1,x1);
00251   // Before the forward operator A is changed it is recommended that you
00252   // uninitialize *invertibleA first to avoid accidental use of *invertiableA
00253   // while it may be in an inconsistent state from the time between *A changes
00254   // and *invertibleA is explicitly updated.  However, this step is not
00255   // required!
00256   Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
00257   // Change the operator in some major way (perhaps even changing its structure)
00258   opMajorChanger.changeOp(A);
00259   // Recreate the LOWSB object and initialize it from scratch
00260   invertibleA = lowsFactory.createOp();
00261   Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
00262   // Solve another set of linear systems
00263   Thyra::assign(x2,Scalar(0.0));
00264   solve(*invertibleA,Thyra::NOTRANS,b2,x2);
00265 } // end solveMajorChangeSolve
00266 
00267 //
00268 // Individual externally preconditioned use cases
00269 //
00270 
00271 // begin createGeneralPreconditionedLinearOpWithSolve
00277 template<class Scalar>
00278 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00279 createGeneralPreconditionedLinearOpWithSolve(
00280   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >        &A
00281   ,const Teuchos::RefCountPtr<const Thyra::PreconditionerBase<Scalar> > &P
00282   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                    &lowsFactory
00283   ,Teuchos::FancyOStream                                                &out
00284   )
00285 {
00286   Teuchos::OSTab tab(out);
00287   out << "\nCreating an externally preconditioned LinearOpWithSolveBase object ...\n";
00288   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00289     invertibleA = lowsFactory.createOp();
00290   Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
00291   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00292   return invertibleA;
00293 } // end createGeneralPreconditionedLinearOpWithSolve
00294 
00295 // begin createUnspecifiedPreconditionedLinearOpWithSolve
00301 template<class Scalar>
00302 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00303 createUnspecifiedPreconditionedLinearOpWithSolve(
00304   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >      &A
00305   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &P_op
00306   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00307   ,Teuchos::FancyOStream                                              &out
00308   )
00309 {
00310   Teuchos::OSTab tab(out);
00311   out << "\nCreating an LinearOpWithSolveBase object given a preconditioner operator not targeted to the left or right ...\n";
00312   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00313     invertibleA = lowsFactory.createOp();
00314   Thyra::initializePreconditionedOp<Scalar>(
00315     lowsFactory,A,Thyra::unspecifiedPrec<Scalar>(P_op)
00316     ,&*invertibleA
00317     );
00318   // Above, the lowsFactory object will decide whether to apply the single
00319   // preconditioner operator on the left or on the right.
00320   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00321   return invertibleA;
00322 } // end createUnspecifiedPreconditionedLinearOpWithSolve
00323 
00324 // begin createLeftPreconditionedLinearOpWithSolve
00330 template<class Scalar>
00331 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00332 createLeftPreconditionedLinearOpWithSolve(
00333   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >      &A
00334   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &P_op_left
00335   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00336   ,Teuchos::FancyOStream                                              &out
00337   )
00338 {
00339   Teuchos::OSTab tab(out);
00340   out << "\nCreating an LinearOpWithSolveBase object given a left preconditioner operator ...\n";
00341   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00342     invertibleA = lowsFactory.createOp();
00343   Thyra::initializePreconditionedOp<Scalar>(
00344     lowsFactory,A,Thyra::leftPrec<Scalar>(P_op_left)
00345     ,&*invertibleA
00346     );
00347   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00348   return invertibleA;
00349 } // end createLeftPreconditionedLinearOpWithSolve
00350 
00351 // begin createRightPreconditionedLinearOpWithSolve
00357 template<class Scalar>
00358 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00359 createRightPreconditionedLinearOpWithSolve(
00360   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >      &A
00361   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &P_op_right
00362   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00363   ,Teuchos::FancyOStream                                              &out
00364   )
00365 {
00366   Teuchos::OSTab tab(out);
00367   out << "\nCreating an LinearOpWithSolveBase object given a right preconditioner operator ...\n";
00368   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00369     invertibleA = lowsFactory.createOp();
00370   Thyra::initializePreconditionedOp<Scalar>(
00371     lowsFactory,A,Thyra::rightPrec<Scalar>(P_op_right)
00372     ,&*invertibleA
00373     );
00374   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00375   return invertibleA;
00376 } // end createRightPreconditionedLinearOpWithSolve
00377 
00378 // begin createLeftRightPreconditionedLinearOpWithSolve
00384 template<class Scalar>
00385 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00386 createLeftRightPreconditionedLinearOpWithSolve(
00387   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >      &A
00388   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &P_op_left
00389   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &P_op_right
00390   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00391   ,Teuchos::FancyOStream                                              &out
00392   )
00393 {
00394   Teuchos::OSTab tab(out);
00395   out << "\nCreating an LinearOpWithSolveBase object given a left and right preconditioner operator ...\n";
00396   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00397     invertibleA = lowsFactory.createOp();
00398   Thyra::initializePreconditionedOp<Scalar>(
00399     lowsFactory,A,Thyra::splitPrec<Scalar>(P_op_left,P_op_right)
00400     ,&*invertibleA
00401     );
00402   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00403   return invertibleA;
00404 } // end createLeftRightPreconditionedLinearOpWithSolve
00405 
00406 // begin createMatrixPreconditionedLinearOpWithSolve
00412 template<class Scalar>
00413 Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00414 createMatrixPreconditionedLinearOpWithSolve(
00415   const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >      &A
00416   ,const Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> >     &A_approx
00417   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00418   ,Teuchos::FancyOStream                                              &out
00419   )
00420 {
00421   Teuchos::OSTab tab(out);
00422   out << "\nCreating a LinearOpWithSolveBase object given an approximate forward"
00423       << " operator to define the preconditioner ...\n";
00424   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00425     invertibleA = lowsFactory.createOp();
00426   Thyra::initializeApproxPreconditionedOp<Scalar>(lowsFactory,A,A_approx,&*invertibleA);
00427   out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
00428   return invertibleA;
00429 } // end createMatrixPreconditionedLinearOpWithSolve
00430 
00431 // begin externalPreconditionerReuseWithSolves
00435 template<class Scalar>
00436 void externalPreconditionerReuseWithSolves(
00437   Thyra::LinearOpBase<Scalar>                                         *A_inout
00438   ,const Thyra::LinearOpChanger<Scalar>                               &opChanger
00439   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>                  &lowsFactory
00440   ,const Thyra::PreconditionerFactoryBase<Scalar>                     &precFactory
00441   ,const Thyra::VectorBase<Scalar>                                    &b1
00442   ,Thyra::VectorBase<Scalar>                                          *x1
00443   ,const Thyra::VectorBase<Scalar>                                    &b2
00444   ,Thyra::VectorBase<Scalar>                                          *x2
00445   ,Teuchos::FancyOStream                                              &out
00446   )
00447 {
00448   Teuchos::OSTab tab(out);
00449   out << "\nShowing resuse of the preconditioner ...\n";
00450   Teuchos::RefCountPtr<Thyra::LinearOpBase<Scalar> >
00451     A = rcp(A_inout,false);
00452   // Create the initial preconditioner for the input forward operator
00453   Teuchos::RefCountPtr<Thyra::PreconditionerBase<Scalar> >
00454     P = precFactory.createPrec();
00455   Thyra::initializePrec<Scalar>(precFactory,A,&*P);
00456   // Create the invertible LOWS object given the preconditioner
00457   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00458     invertibleA = lowsFactory.createOp();
00459   Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
00460   // Solve the first linear system
00461   assign(&*x1,Teuchos::ScalarTraits<Scalar>::zero());
00462   Thyra::SolveStatus<Scalar>
00463     status1 = Thyra::solve(*invertibleA,Thyra::NOTRANS,b1,x1);
00464   out << "\nSolve status:\n" << status1;
00465   // Change the forward linear operator without changing the preconditioner
00466   opChanger.changeOp(&*A);
00467   // Warning!  After the above change the integrity of the preconditioner
00468   // linear operators in P is undefined.  For some implementations of the
00469   // preconditioner, its behavior will remain unchanged (e.g. ILU) which in
00470   // other cases the behavior will change but the preconditioner will still
00471   // work (e.g. Jacobi).  However, there may be valid implementations where
00472   // the preconditioner will simply break if the forward operator that it is
00473   // based on breaks.
00474   //
00475   // Reinitialize the LOWS object given the updated forward operator A and the
00476   // old preconditioner P.
00477   Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
00478   // Solve the second linear system
00479   assign(&*x2,Teuchos::ScalarTraits<Scalar>::zero());
00480   Thyra::SolveStatus<Scalar>
00481     status2 = Thyra::solve(*invertibleA,Thyra::NOTRANS,b2,x2);
00482   out << "\nSolve status:\n" << status2;
00483 } // end externalPreconditionerReuseWithSolves
00484 
00485 //
00486 // Combined use cases
00487 //
00488 
00494 template<class Scalar>
00495 void nonExternallyPreconditionedLinearSolveUseCases(
00496   const Thyra::LinearOpBase<Scalar>                    &A
00497   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00498   ,bool                                                supportsAdjoints
00499   ,Teuchos::FancyOStream                               &out
00500   )
00501 {
00502   Teuchos::OSTab tab(out);
00503   out << "\nRunning example use cases for a LinearOpWithSolveFactoryBase object ...\n";
00504   // Create the RHS (which is just a random set of coefficients)
00505   Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
00506     b1 = Thyra::createMember(A.range()),
00507     b2 = Thyra::createMember(A.range());
00508   Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b1 );
00509   Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b2 );
00510   // Create the LHS for the linear solve
00511   Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
00512     x1 = Thyra::createMember(A.domain()),
00513     x2 = Thyra::createMember(A.domain());
00514   // Perform a single, non-adjoint, linear solve
00515   singleLinearSolve(A,lowsFactory,*b1,&*x1,out);
00516   // Creating a scaled adjoint LinearOpWithSolveBase object
00517   if(supportsAdjoints) {
00518     Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00519       invertibleAdjointA = createScaledAdjointLinearOpWithSolve(
00520         Teuchos::rcp(&A,false),Scalar(2.0),lowsFactory,out);
00521   }
00522   // Perform a solve, change the operator, and then solve again.
00523   solveNumericalChangeSolve(
00524     const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
00525     ,Thyra::NullLinearOpChanger<Scalar>()        // This object will not really change A!
00526     ,lowsFactory,*b1,&*x1,*b2,&*x1,out
00527     );
00528   // Perform a solve, change the operator in a very small way, and then solve again.
00529   solveSmallNumericalChangeSolve(
00530     const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
00531     ,Thyra::NullLinearOpChanger<Scalar>()        // This object will not really change A!
00532     ,lowsFactory,*b1,&*x1,*b2,&*x1,out
00533     );
00534   // Perform a solve, change the operator in a major way, and then solve again.
00535   solveMajorChangeSolve(
00536     const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
00537     ,Thyra::NullLinearOpChanger<Scalar>()        // This object will not really change A!
00538     ,lowsFactory,*b1,&*x1,*b2,&*x1,out
00539     );
00540 }
00541 
00547 template<class Scalar>
00548 void externallyPreconditionedLinearSolveUseCases(
00549   const Thyra::LinearOpBase<Scalar>                    &A
00550   ,const Thyra::LinearOpWithSolveFactoryBase<Scalar>   &lowsFactory
00551   ,const Thyra::PreconditionerFactoryBase<Scalar>      &precFactory
00552   ,const bool                                          supportsLeftPrec
00553   ,const bool                                          supportsRightPrec
00554   ,Teuchos::FancyOStream                               &out
00555   )
00556 {
00557   Teuchos::OSTab tab(out);
00558   out << "\nRunning example use cases with an externally defined"
00559       << " preconditioner with a LinearOpWithSolveFactoryBase object ...\n";
00560   // Create the RHS (which is just a random set of coefficients)
00561   Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
00562     b1 = Thyra::createMember(A.range()),
00563     b2 = Thyra::createMember(A.range());
00564   Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b1 );
00565   Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b2 );
00566   // Create the LHS for the linear solve
00567   Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
00568     x1 = Thyra::createMember(A.domain()),
00569     x2 = Thyra::createMember(A.domain());
00570   // Create a preconditioner for the input forward operator
00571   Teuchos::RefCountPtr<Thyra::PreconditionerBase<Scalar> >
00572     P = precFactory.createPrec();
00573   Thyra::initializePrec<Scalar>(precFactory,Teuchos::rcp(&A,false),&*P);
00574   // Above, we don't really know the nature of the preconditioner.  It could a
00575   // single linear operator to be applied on the left or the right or it could
00576   // be a split preconditioner with different linear operators to be applied
00577   // on the right or left.  Or, it could be a single linear operator that is
00578   // not targeted to the left or the right.
00579   //
00580   // Create a LOWSB object given the created preconditioner
00581   Teuchos::RefCountPtr<Thyra::LinearOpWithSolveBase<Scalar> >
00582     invertibleA = createGeneralPreconditionedLinearOpWithSolve<Scalar>(
00583       Teuchos::rcp(&A,false),P,lowsFactory,out);
00584   // Grab a preconditioner operator out of the preconditioner object
00585   Teuchos::RefCountPtr<const Thyra::LinearOpBase<Scalar> > P_op;
00586   if((P_op=P->getUnspecifiedPrecOp()).get());
00587   else if((P_op=P->getLeftPrecOp()).get());
00588   else if((P_op=P->getRightPrecOp()).get());
00589   // Create a LOWSB object given an unspecified preconditioner operator
00590   invertibleA = createUnspecifiedPreconditionedLinearOpWithSolve(
00591     Teuchos::rcp(&A,false),P_op,lowsFactory,out);
00592   // Create a LOWSB object given a left preconditioner operator
00593   if(supportsLeftPrec) {
00594     invertibleA = createLeftPreconditionedLinearOpWithSolve(
00595       Teuchos::rcp(&A,false),P_op,lowsFactory,out);
00596   }
00597   // Create a LOWSB object given a right preconditioner operator
00598   if(supportsRightPrec) {
00599     invertibleA = createRightPreconditionedLinearOpWithSolve(
00600       Teuchos::rcp(&A,false),P_op,lowsFactory,out);
00601   }
00602   // Create a LOWSB object given (bad set of) left and right preconditioner
00603   // operators
00604   if( supportsLeftPrec && supportsRightPrec ) {
00605     invertibleA = createLeftRightPreconditionedLinearOpWithSolve(
00606       Teuchos::rcp(&A,false),P_op,P_op,lowsFactory,out);
00607   }
00608   // Create a LOWSB object given a (very good) approximate forward linear
00609   // operator to construct the preconditoner from..
00610   invertibleA = createMatrixPreconditionedLinearOpWithSolve<Scalar>(
00611     Teuchos::rcp(&A,false),Teuchos::rcp(&A,false),lowsFactory,out);
00612   // Preconditioner reuse example
00613   externalPreconditionerReuseWithSolves(
00614     const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
00615     ,Thyra::NullLinearOpChanger<Scalar>()        // This object will not really change A!
00616     ,lowsFactory,precFactory
00617     ,*b1,&*x1,*b2,&*x2,out
00618     );
00619 }
00620 
00621 #endif // THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1