Thyra_LinearOpWithSolveFactoryExamples.hpp

// @HEADER
// ***********************************************************************
// 
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//                 Copyright (2004) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP
#define THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP


#include "Thyra_LinearOpWithSolveFactoryBase.hpp"
#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
#include "Thyra_PreconditionerFactoryHelpers.hpp"
#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
#include "Thyra_DefaultPreconditioner.hpp"


namespace Thyra {


//
// Helper code
//

template<class Scalar>
class LinearOpChanger {
public:
  virtual ~LinearOpChanger() {}
  virtual void changeOp( LinearOpBase<Scalar> *op ) const = 0;
};

template<class Scalar>
class NullLinearOpChanger : public LinearOpChanger<Scalar> {
public:
  void changeOp( LinearOpBase<Scalar> *op ) const {}
};


} // namespace Thyra


//
// Individual non-externally preconditioned use cases
//


// begin singleLinearSolve
template<class Scalar>
void singleLinearSolve(
  const Thyra::LinearOpBase<Scalar> &A,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::VectorBase<Scalar> &b,
  Thyra::VectorBase<Scalar> *x,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nPerforming a single linear solve ...\n";
  // Create the LOWSB object that will be used to solve the linear system
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = Thyra::linearOpWithSolve(lowsFactory,Teuchos::rcp(&A,false));
  // Solve the system using a default solve criteria using a non-member helper function 
  assign(&*x,Teuchos::ScalarTraits<Scalar>::zero()); // Must initialize to a guess before solve!
  Thyra::SolveStatus<Scalar>
    status = Thyra::solve(*invertibleA,Thyra::NOTRANS,b,x);
  out << "\nSolve status:\n" << status;
} // end singleLinearSolve


// begin createScaledAdjointLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createScaledAdjointLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Scalar &scalar,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating a scaled adjoint LinearOpWithSolveBase object ...\n";
  // Create the LOWSB object that will be used to solve the linear system
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleAdjointA
    = Thyra::linearOpWithSolve(lowsFactory,scale(scalar,adjoint(A)));
  out << "\nCreated LOWSB object:\n" << describe(*invertibleAdjointA,Teuchos::VERB_MEDIUM);
  return invertibleAdjointA;
} // end createScaledAdjointLinearOpWithSolve


// begin solveNumericalChangeSolve
template<class Scalar>
void solveNumericalChangeSolve(
  Thyra::LinearOpBase<Scalar> *A,
  const Thyra::LinearOpChanger<Scalar> &opChanger,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::VectorBase<Scalar> &b1,
  Thyra::VectorBase<Scalar> *x1,
  const Thyra::VectorBase<Scalar> &b2,
  Thyra::VectorBase<Scalar> *x2,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nPerforming a solve, changing the operator, then performing another solve ...\n";
  // Get a local non-owned RCP to A to be used by lowsFactory
  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
    rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
  // Create the LOWSB object that will be used to solve the linear system
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  // Initialize the invertible linear operator given the forward operator
  Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Solve the system using a default solve criteria using a non-member helper function
  Thyra::assign(x1,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b1,x1);
  // Before the forward operator A is changed it is recommended that you
  // uninitialize *invertibleA first to avoid accidental use of *invertiableA
  // while it may be in an inconsistent state from the time between *A changes
  // and *invertibleA is explicitly updated. However, this step is not
  // required!
  Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
  // Change the operator and reinitialize the invertible operator
  opChanger.changeOp(A);
  Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Note that above will reuse any factorization structures that may have been
  // created in the first call to initializeOp(...).
  // Finally, solve another linear system with new values of A
  Thyra::assign(x2,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b2,x2);
} // end solveNumericalChangeSolve


// begin solveSmallNumericalChangeSolve
template<class Scalar>
void solveSmallNumericalChangeSolve(
  Thyra::LinearOpBase<Scalar> *A,
  const Thyra::LinearOpChanger<Scalar> &opSmallChanger,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::VectorBase<Scalar> &b1,
  Thyra::VectorBase<Scalar> *x1,
  const Thyra::VectorBase<Scalar> &b2,
  Thyra::VectorBase<Scalar> *x2,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nPerforming a solve, changing the operator in a very small way, then performing another solve ...\n";
  // Get a local non-owned RCP to A to be used by lowsFactory
  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
    rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
  // Create the LOWSB object that will be used to solve the linear system
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  // Initialize the invertible linear operator given the forward operator
  Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Solve the system using a default solve criteria using a non-member helper function
  Thyra::assign(x1,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b1,x1);
  // Before the forward operator A is changed it is recommended that you
  // uninitialize *invertibleA first to avoid accidental use of *invertiableA
  // while it may be in an inconsistent state from the time between *A changes
  // and *invertibleA is explicitly updated. However, this step is not
  // required!
  Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
  // Change the operator and reinitialize the invertible operator
  opSmallChanger.changeOp(A);
  Thyra::initializeAndReuseOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Note that above a maximum amount of reuse will be achieved, such as
  // keeping the same precondtioner.
  Thyra::assign(x2,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b2,x2);
} // end solveSmallNumericalChangeSolve


// begin solveMajorChangeSolve
template<class Scalar>
void solveMajorChangeSolve(
  Thyra::LinearOpBase<Scalar> *A,
  const Thyra::LinearOpChanger<Scalar> &opMajorChanger,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::VectorBase<Scalar> &b1,
  Thyra::VectorBase<Scalar> *x1,
  const Thyra::VectorBase<Scalar> &b2,
  Thyra::VectorBase<Scalar> *x2,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nPerforming a solve, changing the operator in a major way, then performing another solve ...\n";
  // Get a local non-owned RCP to A to be used by lowsFactory
  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
    rcpA = Teuchos::rcp(A,false); // Note: This is the right way to do this!
  // Create the LOWSB object that will be used to solve the linear system
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  // Initialize the invertible linear operator given the forward operator
  Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Solve the system using a default solve criteria using a non-member helper function
  Thyra::assign(x1,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b1,x1);
  // Before the forward operator A is changed it is recommended that you
  // uninitialize *invertibleA first to avoid accidental use of *invertiableA
  // while it may be in an inconsistent state from the time between *A changes
  // and *invertibleA is explicitly updated. However, this step is not
  // required!
  Thyra::uninitializeOp<Scalar>(lowsFactory,&*invertibleA);
  // Change the operator in some major way (perhaps even changing its structure)
  opMajorChanger.changeOp(A);
  // Recreate the LOWSB object and initialize it from scratch
  invertibleA = lowsFactory.createOp();
  Thyra::initializeOp<Scalar>(lowsFactory,rcpA,&*invertibleA);
  // Solve another set of linear systems
  Thyra::assign(x2,Scalar(0.0));
  solve(*invertibleA,Thyra::NOTRANS,b2,x2);
} // end solveMajorChangeSolve


//
// Individual externally preconditioned use cases
//


// begin createGeneralPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createGeneralPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::PreconditionerBase<Scalar> > &P,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating an externally preconditioned LinearOpWithSolveBase object ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createGeneralPreconditionedLinearOpWithSolve


// begin createUnspecifiedPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createUnspecifiedPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating an LinearOpWithSolveBase object given a preconditioner operator not targeted to the left or right ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(
    lowsFactory,A,Thyra::unspecifiedPrec<Scalar>(P_op)
    ,&*invertibleA
    );
  // Above, the lowsFactory object will decide whether to apply the single
  // preconditioner operator on the left or on the right.
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createUnspecifiedPreconditionedLinearOpWithSolve


// begin createLeftPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createLeftPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_left,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating an LinearOpWithSolveBase object given a left preconditioner"
      << " operator ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(
    lowsFactory,A,Thyra::leftPrec<Scalar>(P_op_left)
    ,&*invertibleA
    );
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createLeftPreconditionedLinearOpWithSolve


// begin createRightPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createRightPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_right,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating an LinearOpWithSolveBase object given a right preconditioner operator ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(
    lowsFactory,A,Thyra::rightPrec<Scalar>(P_op_right)
    ,&*invertibleA
    );
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createRightPreconditionedLinearOpWithSolve


// begin createLeftRightPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createLeftRightPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_left,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &P_op_right,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating an LinearOpWithSolveBase object given a left and"
      << "right preconditioner operator ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(
    lowsFactory,A,Thyra::splitPrec<Scalar>(P_op_left,P_op_right)
    ,&*invertibleA
    );
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createLeftRightPreconditionedLinearOpWithSolve


// begin createMatrixPreconditionedLinearOpWithSolve
template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
createMatrixPreconditionedLinearOpWithSolve(
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &A_approx,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nCreating a LinearOpWithSolveBase object given an approximate forward"
      << " operator to define the preconditioner ...\n";
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializeApproxPreconditionedOp<Scalar>(lowsFactory,A,A_approx,&*invertibleA);
  out << "\nCreated LOWSB object:\n" << describe(*invertibleA,Teuchos::VERB_MEDIUM);
  return invertibleA;
} // end createMatrixPreconditionedLinearOpWithSolve


// begin externalPreconditionerReuseWithSolves
template<class Scalar>
void externalPreconditionerReuseWithSolves(
  Thyra::LinearOpBase<Scalar> *A_inout,
  const Thyra::LinearOpChanger<Scalar> &opChanger,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::PreconditionerFactoryBase<Scalar> &precFactory,
  const Thyra::VectorBase<Scalar> &b1,
  Thyra::VectorBase<Scalar> *x1,
  const Thyra::VectorBase<Scalar> &b2,
  Thyra::VectorBase<Scalar> *x2,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nShowing resuse of the preconditioner ...\n";
  Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
    A = rcp(A_inout,false);
  // Create the initial preconditioner for the input forward operator
  Teuchos::RCP<Thyra::PreconditionerBase<Scalar> >
    P = precFactory.createPrec();
  Thyra::initializePrec<Scalar>(precFactory,A,&*P);
  // Create the invertible LOWS object given the preconditioner
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = lowsFactory.createOp();
  Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
  // Solve the first linear system
  assign(&*x1,Teuchos::ScalarTraits<Scalar>::zero());
  Thyra::SolveStatus<Scalar>
    status1 = Thyra::solve(*invertibleA,Thyra::NOTRANS,b1,x1);
  out << "\nSolve status:\n" << status1;
  // Change the forward linear operator without changing the preconditioner
  opChanger.changeOp(&*A);
  // Warning! After the above change the integrity of the preconditioner
  // linear operators in P is undefined. For some implementations of the
  // preconditioner, its behavior will remain unchanged (e.g. ILU) which in
  // other cases the behavior will change but the preconditioner will still
  // work (e.g. Jacobi). However, there may be valid implementations where
  // the preconditioner will simply break if the forward operator that it is
  // based on breaks.
  //
  // Reinitialize the LOWS object given the updated forward operator A and the
  // old preconditioner P.
  Thyra::initializePreconditionedOp<Scalar>(lowsFactory,A,P,&*invertibleA);
  // Solve the second linear system
  assign(&*x2,Teuchos::ScalarTraits<Scalar>::zero());
  Thyra::SolveStatus<Scalar>
    status2 = Thyra::solve(*invertibleA,Thyra::NOTRANS,b2,x2);
  out << "\nSolve status:\n" << status2;
} // end externalPreconditionerReuseWithSolves


//
// Combined use cases
//


template<class Scalar>
void nonExternallyPreconditionedLinearSolveUseCases(
  const Thyra::LinearOpBase<Scalar> &A,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  bool supportsAdjoints,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nRunning example use cases for a LinearOpWithSolveFactoryBase object ...\n";
  // Create the RHS (which is just a random set of coefficients)
  Teuchos::RCP<Thyra::VectorBase<Scalar> >
    b1 = Thyra::createMember(A.range()),
    b2 = Thyra::createMember(A.range());
  Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b1 );
  Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b2 );
  // Create the LHS for the linear solve
  Teuchos::RCP<Thyra::VectorBase<Scalar> >
    x1 = Thyra::createMember(A.domain()),
    x2 = Thyra::createMember(A.domain());
  // Perform a single, non-adjoint, linear solve
  singleLinearSolve(A,lowsFactory,*b1,&*x1,out);
  // Creating a scaled adjoint LinearOpWithSolveBase object
  if(supportsAdjoints) {
    Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
      invertibleAdjointA = createScaledAdjointLinearOpWithSolve(
        Teuchos::rcp(&A,false),Scalar(2.0),lowsFactory,out);
  }
  // Perform a solve, change the operator, and then solve again.
  solveNumericalChangeSolve(
    const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
    ,Thyra::NullLinearOpChanger<Scalar>() // This object will not really change A!
    ,lowsFactory,*b1,&*x1,*b2,&*x1,out
    );
  // Perform a solve, change the operator in a very small way, and then solve again.
  solveSmallNumericalChangeSolve(
    const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
    ,Thyra::NullLinearOpChanger<Scalar>() // This object will not really change A!
    ,lowsFactory,*b1,&*x1,*b2,&*x1,out
    );
  // Perform a solve, change the operator in a major way, and then solve again.
  solveMajorChangeSolve(
    const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
    ,Thyra::NullLinearOpChanger<Scalar>() // This object will not really change A!
    ,lowsFactory,*b1,&*x1,*b2,&*x1,out
    );
}


template<class Scalar>
void externallyPreconditionedLinearSolveUseCases(
  const Thyra::LinearOpBase<Scalar> &A,
  const Thyra::LinearOpWithSolveFactoryBase<Scalar> &lowsFactory,
  const Thyra::PreconditionerFactoryBase<Scalar> &precFactory,
  const bool supportsLeftPrec,
  const bool supportsRightPrec,
  Teuchos::FancyOStream &out
  )
{
  Teuchos::OSTab tab(out);
  out << "\nRunning example use cases with an externally defined"
      << " preconditioner with a LinearOpWithSolveFactoryBase object ...\n";
  // Create the RHS (which is just a random set of coefficients)
  Teuchos::RCP<Thyra::VectorBase<Scalar> >
    b1 = Thyra::createMember(A.range()),
    b2 = Thyra::createMember(A.range());
  Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b1 );
  Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*b2 );
  // Create the LHS for the linear solve
  Teuchos::RCP<Thyra::VectorBase<Scalar> >
    x1 = Thyra::createMember(A.domain()),
    x2 = Thyra::createMember(A.domain());
  // Create a preconditioner for the input forward operator
  Teuchos::RCP<Thyra::PreconditionerBase<Scalar> >
    P = precFactory.createPrec();
  Thyra::initializePrec<Scalar>(precFactory,Teuchos::rcp(&A,false),&*P);
  // Above, we don't really know the nature of the preconditioner. It could a
  // single linear operator to be applied on the left or the right or it could
  // be a split preconditioner with different linear operators to be applied
  // on the right or left. Or, it could be a single linear operator that is
  // not targeted to the left or the right.
  //
  // Create a LOWSB object given the created preconditioner
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
    invertibleA = createGeneralPreconditionedLinearOpWithSolve<Scalar>(
      Teuchos::rcp(&A,false),P,lowsFactory,out);
  // Grab a preconditioner operator out of the preconditioner object
  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > P_op;
  if((P_op=P->getUnspecifiedPrecOp()).get());
  else if((P_op=P->getLeftPrecOp()).get());
  else if((P_op=P->getRightPrecOp()).get());
  // Create a LOWSB object given an unspecified preconditioner operator
  invertibleA = createUnspecifiedPreconditionedLinearOpWithSolve(
    Teuchos::rcp(&A,false),P_op,lowsFactory,out);
  // Create a LOWSB object given a left preconditioner operator
  if(supportsLeftPrec) {
    invertibleA = createLeftPreconditionedLinearOpWithSolve(
      Teuchos::rcp(&A,false),P_op,lowsFactory,out);
  }
  // Create a LOWSB object given a right preconditioner operator
  if(supportsRightPrec) {
    invertibleA = createRightPreconditionedLinearOpWithSolve(
      Teuchos::rcp(&A,false),P_op,lowsFactory,out);
  }
  // Create a LOWSB object given (bad set of) left and right preconditioner
  // operators
  if( supportsLeftPrec && supportsRightPrec ) {
    invertibleA = createLeftRightPreconditionedLinearOpWithSolve(
      Teuchos::rcp(&A,false),P_op,P_op,lowsFactory,out);
  }
  // Create a LOWSB object given a (very good) approximate forward linear
  // operator to construct the preconditoner from..
  invertibleA = createMatrixPreconditionedLinearOpWithSolve<Scalar>(
    Teuchos::rcp(&A,false),Teuchos::rcp(&A,false),lowsFactory,out);
  // Preconditioner reuse example
  externalPreconditionerReuseWithSolves(
    const_cast<Thyra::LinearOpBase<Scalar>*>(&A) // Don't worry, it will not be changed!
    ,Thyra::NullLinearOpChanger<Scalar>() // This object will not really change A!
    ,lowsFactory,precFactory
    ,*b1,&*x1,*b2,&*x2,out
    );
}


#endif // THYRA_LINEAR_OP_WITH_SOLVE_EXAMPLES_HPP

---