Thyra::LinearOpWithSolveFactoryBase< Scalar > Class Template Reference

Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects. More...

#include <Thyra_LinearOpWithSolveFactoryBase_decl.hpp>

Inheritance diagram for Thyra::LinearOpWithSolveFactoryBase< Scalar >:

[legend]List of all members.

Preconditioner Factory Management
virtual bool	acceptsPreconditionerFactory () const
	Determines if *this accepts external preconditioner factories.
virtual void	setPreconditionerFactory (const RCP< PreconditionerFactoryBase< Scalar > > &precFactory, const std::string &precFactoryName)
	Set a preconditioner factory object.
virtual RCP< PreconditionerFactoryBase< Scalar > >	getPreconditionerFactory () const
	Get a preconditioner factory object.
virtual void	unsetPreconditionerFactory (RCP< PreconditionerFactoryBase< Scalar > > *precFactory=NULL, std::string *precFactoryName=NULL)
	Unset the preconditioner factory (if one is set).
Creation/Initialization of basic LinearOpWithSolveBase objects
virtual bool	isCompatible (const LinearOpSourceBase< Scalar > &fwdOpSrc) const =0
	Check that a LinearOpBase object is compatible with *this factory object.
virtual RCP< LinearOpWithSolveBase< Scalar > >	createOp () const =0
	Create an (uninitialized) LinearOpWithSolveBase object to be initialized later in this->initializeOp().
virtual void	initializeOp (const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED) const =0
	Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object.
virtual void	initializeAndReuseOp (const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op) const
	Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object but allow for reuse of any preprocessing that is in *Op.
virtual void	uninitializeOp (LinearOpWithSolveBase< Scalar > *Op, RCP< const LinearOpSourceBase< Scalar > > *fwdOpSrc=NULL, RCP< const PreconditionerBase< Scalar > > *prec=NULL, RCP< const LinearOpSourceBase< Scalar > > *approxFwdOpSrc=NULL, ESupportSolveUse *supportSolveUse=NULL) const =0
	Uninitialize a LinearOpWithSolveBase object and return its remembered forward linear operator and potentially also its externally generated preconditioner.
Creation/Initialization of Preconditioned LinearOpWithSolveBase objects
virtual bool	supportsPreconditionerInputType (const EPreconditionerInputType precOpType) const
	Determines if *this supports given preconditioner type.
virtual void	initializePreconditionedOp (const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const RCP< const PreconditionerBase< Scalar > > &prec, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED) const
	Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object and an optional PreconditionerBase object.
virtual void	initializeApproxPreconditionedOp (const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED) const
	Initialize a pre-created LinearOpWithSolveBase object given a "compatible" forward LinearOpBase object and an approximate forward LinearOpBase object.

---

Detailed Description

template<class Scalar>
class Thyra::LinearOpWithSolveFactoryBase< Scalar >

Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects.

Outline

• Introduction
• Use cases
• Output and verbosity
• Notes to subclass developers

Introduction

This strategy interface allows a client to take one or more "compatible" LinearOpBase objects and then create one or more LinearOpWithSolveBase objects that can then be used to solve for linear systems. This interface carefully separates the construction from the initialization of a LinearOpWithSolveBase object.

Note that the non-member functions defined here provide for simpler use cases and are demonstrated in the example code below.

This interface can be implemented by both direct and iterative linear solvers.

This interface supports the concept of an external preconditioner that can be created by the client an passed into the function initializePreconditionedOp()

Use cases

The following use cases demonstrate and specify the behavior of this interface in several different situations. These use cases don't cover very possible variation but implementors and clients should get a pretty good idea of the desired behavior from what is specified below.

Note: All of the code fragments shown in the below use cases is code that is compiled and run automatically by the test harness and the code fragments are automatically pulled into this HTML documentation whenever the documentation is rebuilt. Therefore, one can have a very high degree of confidence that the code shown in these examples is correct.

The following use cases are described below:

• Performing a single linear solve given a forward operator
• Creating invertible operators for scaled and/or adjoint forward operators
• Updates of linear operator between linear solves
• Reuse of factorizations for small changes in the forward operator
• Updating the linear solver for major changes in the structure of the forward operator
• Externally defined preconditioners
  • Use of externally defined preconditioner operators
  • The use of externally defined preconditioner matrices
  • Reuse of externally defined preconditioner operators
  • Reuse of externally defined preconditioner matrices

Performing a single linear solve given a forward operator

Performing a single linear solve is at minimum a two step process. The following example function shows how to take a LinearOpBase object, a compatible LinearOpWithSolveFactoryBase object, a RHS and a LHS and then solve the linear system using the default solve criteria:

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

See the documentation for the LinearOpWithSolveBase interface for how to specify solve criteria, how to perform adjoint (i.e. transpose) solve, and how to solve systems with multiple RHSs.

Note that the forward operator A is passed into this function singleLinearSolve(...) as a raw object reference and not as a Teuchos::RCP wrapped object since no persisting relationship with this object will be last after this function exists, even if an exception is thrown.

Also note that once the above function exits that all memory of the invertible operator created as part of the invertibleA object will be erased. The following use cases show more sophisticated uses of these interfaces.

Creating invertible operators for scaled and/or adjoint forward operators

This interface requires that all good implementations support implicitly scaled and adjoint (or transposed) forward operators. The following example function shows how this looks:

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

In the above example, the functions adjoint() and scale() create an implicitly scaled adjoint operator of type DefaultScaledAdjointLinearOp which is then unwrapped by the lowsFactory implementation. The idea is that any operation that works with a particular forward operator A should automatically work with an implicitly scaled and/or adjoint view of that forward operator. The specification of this interface actually says that all LinearOpBase objects that support the abstract mix-in interface ScaledAdjointLinearOpBase will allow this feature.

Above also note that the forward operator A is passed in as a Teuchos::RCP wrapped object since it will be used to create a persisting relationship with the returned Thyra::LinearOpWithSolveBase object. Also note that the lowsFactory object is still passed in as a raw object reference since no persisting relationship with lowsFactory is created as a side effect of calling this function. Remember, the specification of this interface requires that the returned LinearOpWithSolveBase objects be independent from the lowsFactory object that creates them. In other words, once a LinearOpWithSolveFactoryBase object creates a LinearOpWithSolveBase object, then the LinearOpWithSolveFactoryBase object can be deleted and the LinearOpWithSolveBase it created must still be valid.

Updates of linear operator between linear solves

In this use case is comprised of:

• Creating a LinearOpWithSolveBase object from a LinearOpBase object
• Using the LinearOpWithSolveBase object to solve one or more linear systems
• Modifying the LinearOpBase object and reinitializing the LinearOpWithSolveBase object
• Using the updated LinearOpWithSolveBase object to solve one or more linear systems

The below code fragment shows how this looks.

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

In the above code fragment the call to the function lowsFactory.uninitializeOp(&*invertibleA) may not fully uninitialize the *invertibleA object as it may contain memory of a factorization structure, or a communication pattern, or whatever may have been computed in the first call to lowsFactory.initializeOp(rcpA,&*invertibleA). This allows for a more efficient reinitialization on the second call of lowsFactory.initializeOp(rcpA,&*invertibleA).

Reuse of factorizations for small changes in the forward operator

This interface supports the notion of the reuse of all factorizations or other expensive preprocessing used in the last initialization of the LinearOpWithSolveBase object. The primary use for this is the reuse of a preconditioner for small changes in the matrix values. While this approach generally is not very effective in many cases, there are some cases, such as in transient solvers, where this has been shown to be effective in some problems.

The below example function shows what this looks like:

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

Note that the *invertibleA object reinitialized in the second call to lowsFactory.initializeAndReuseOp(rcpA,&*invertibleA) must be able to correctly solve the linear systems to an appropriate accuracy. For example, a preconditioned iterative linear solver might keep the same preconditioner but would use the updated forward operator to define the linear system. A direct solver can not reuse the factorization unless it has an iterative refinement feature (which uses the updated forward operator) or some other device. Or a direct solver might, for instance, reuse the same factorization structure and not re-pivot where it might re-pivot generally.

Updating the linear solver for major changes in the structure of the forward operator

A major change in the shape or structure of a forward operator generally eliminates the possibility of reusing any computations between different calls to initializeOp(). In these instances, the client might as well just recreate the LinearOpWithSolveBase using createOp() before the reinitialization.

The below example function shows what this looks like:

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

Externally defined preconditioners

This interface also supports the use of externally defined preconditioners that are created and controlled by the client. Client-created and client-controlled preconditioners can be passed along with the forward operator through the function initializePreconditionedOp(). Only LinearOpWithSolveFactoryBase objects that return supportsPreconditionerType()==true will support externally defined preconditioners. In general, iterative linear solver implementations will support externally defined preconditioners and direct linear solver implementations will not.

Externally defined preconditioners can be passed in as operators or as matrices and these two sub use cases are described below.

Use of externally defined preconditioner operators

The client can use externally defined preconditioner linear operator(s) for LinearOpWithSolveFactoryBase objects that accept preconditioners. Preconditioner operator(s) are specified through the PreconditionerBase interface. A PreconditionerBase object is essentially an encapsulation of a left, right, or split left/right preconditioner.

Using an externally defined PreconditionerBase object

Given an externally defined PreconditionerBase object, the client can pass it allow with a forward linear operator to initialize a LinearOpWithSolveBase object as shown in the following code fragment:

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

Using a single externally defined LinearOpBase object as an unspecified preconditioner

A client can easily use any compatible appropriately defined LinearOpBase object as a preconditioner operator. To do so, it can just be wrapped by the default PreconditionerBase implementation class DefaultPreconditioner.

The following code fragment shows show to use an externally defined LinearOpBase object as a preconditioner operator without specifying whether it should be applied on the right or on the left.

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

Using a single externally defined LinearOpBase object as a left preconditioner

An externally defined LinearOpBase object can be used as a preconditioner operator targeted as a left preconditioner as shown in the following code fragment:

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

Using a single externally defined LinearOpBase object as a right preconditioner

An externally defined LinearOpBase object can be used as a preconditioner operator targeted as a right preconditioner as shown in the following code fragement:

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

Using two single externally defined LinearOpBase objects as a split left/right preconditioner

Two externally defined LinearOpBase objects can be used as preconditioner operators for form a split preconditioner targeted to the left and right as shown in the following code fragement:

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

The use of externally defined preconditioner matrices

A different linear operator can be used to build a preconditioner internally when initializing a LinearOpWithSolveBase object. The following code fragement shows how to use an approximate forward linear from which to form a preconditioner internally:

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

Reuse of externally defined preconditioner operators

Reusing an externally defined preconditioner is a very straightforward matter conceptually. Since the client controls the creation and initialization of the preconditioner, the client can control which preconditioners are paired with which forward operators in order to form LinearOpWithSolveBase objects. However, when an abstract PreconditionerFactoryBase object is used to create and maintain the preconditioner, the exact nature of the preconditioner after the operator is updated is undefined. The following code fragment show a use case where the operator is reused without explicitly updating the preconditioner between solves:

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

Reuse of externally defined preconditioner matrices

This interface does not guarantee the correct reuse of externally defined preconditioner matrices. The problem is that the internally generated preconditioner may be dependent on the input preconditioner matrix and it is not clear for the current interface design how to cleanly handle this use case. However, when the client knows that the internally generated preconditioner operator is independent from the externally defined preconditioner input matrix, then the function initializeAndReuseOp() can be called to reuse the internal preconditioner but this is implementation defined.

Output and verbosity

TODO: Provide some guidance on how clients and subclasses should interpret Teuchos::EVerbLevel.

Notes to subclass developers

This interface assumes a minimal default set of functionality that is appropriate for direct and simple iterative linear solver implementations. The pure virtual functions that must be overridden by a subclass are isCompatible(), createOp(), initializeOp(), and uninitializeOp(). By far the most complex function to implement is initializeOp() which is where the real guts of the factory method is defined.

If the concrete subclass can support some significant type of preprocessing reuse then it may override to function initializeAndReuseOp(). The most common use of this function it to support the reuse of preconditioners between small changes in forward operator matrix values.

If the concrete subclass can utilize a preconditioner, then it should override the functions supportsPreconditionerInputType() and initializePreconditionedOp(). The subclass implementation can decide if preconditioner operators and/or matrices are supported or not and this is determined by the return value from supportsPreconditionerInputType().

Definition at line 391 of file Thyra_LinearOpWithSolveFactoryBase_decl.hpp.

---

Member Function Documentation

template<class Scalar>

bool Thyra::LinearOpWithSolveFactoryBase< Scalar >::acceptsPreconditionerFactory

(

)

const [virtual]

Determines if *this accepts external preconditioner factories.

The default implementation returns false.

Definition at line 40 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

void Thyra::LinearOpWithSolveFactoryBase< Scalar >::setPreconditionerFactory	(	const RCP< PreconditionerFactoryBase< Scalar > > &	precFactory,
		const std::string &	precFactoryName
	)			[virtual]

Set a preconditioner factory object.

Parameters:

	precFactory	[in] The preconditioner factory to be used internally to create preconditioners.
	precFactoryName	[in] The name to give to the preconditioner factory internally. This name is used when setting parameters in the parameter list.

Preconditions:

• this->acceptsPreconditionerFactory()==true
• precFactory.get()!=NULL

Postconditions:

• this->getPreconditionerFactory().get()==precFactory.get()

The default implementation thrown an exception which is consistent with the default implementation acceptsPreconditionerFactory()==false.

Definition at line 47 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

RCP< PreconditionerFactoryBase< Scalar > > Thyra::LinearOpWithSolveFactoryBase< Scalar >::getPreconditionerFactory

(

)

const [virtual]

Get a preconditioner factory object.

The default implementation returns Teuchos::null.

Definition at line 62 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

void Thyra::LinearOpWithSolveFactoryBase< Scalar >::unsetPreconditionerFactory	(	RCP< PreconditionerFactoryBase< Scalar > > *	precFactory = NULL,
		std::string *	precFactoryName = NULL
	)			[virtual]

Unset the preconditioner factory (if one is set).

Postconditions:

• this->getPreconditionerFactory().get()==NULL

The default implementation returns Teuchos::null.

Definition at line 69 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

virtual bool Thyra::LinearOpWithSolveFactoryBase< Scalar >::isCompatible

(

const LinearOpSourceBase< Scalar > &

fwdOpSrc

)

const [pure virtual]

Check that a LinearOpBase object is compatible with *this factory object.

template<class Scalar>

virtual RCP<LinearOpWithSolveBase<Scalar> > Thyra::LinearOpWithSolveFactoryBase< Scalar >::createOp

(

)

const [pure virtual]

Create an (uninitialized) LinearOpWithSolveBase object to be initialized later in this->initializeOp().

Note that on output return->domain().get()==NULL may be true which means that the operator is not fully initialized. In fact, the output operator object is not guaranteed to be fully initialized until after it is passed through this->initializeOp().

template<class Scalar>

virtual void Thyra::LinearOpWithSolveFactoryBase< Scalar >::initializeOp	(	const RCP< const LinearOpSourceBase< Scalar > > &	fwdOpSrc,
		LinearOpWithSolveBase< Scalar > *	Op,
		const ESupportSolveUse	supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED
	)			const [pure virtual]

Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object.

Parameters:

	fwdOpSrc	[in] The forward linear operator that will be used to create the output LinearOpWithSolveBase object. Note that this object is remembered by the *Op object on output.
	Op	[in/out] The output LinearOpWithSolveBase object. This object must have be created first by this->createOp(). The object may have also already been passed through this function several times. Note that subclasses should always first strip off the transpose and scaling by calling unwrap() before attempting to dynamic cast the object.
	supportSolveUse	[in] Determines if Op->solve(...) or Op->solveTranspose(...) will be called. This allows *this factory object to determine how to best initialize the *Op object. Default supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED

Preconditions:

• fwdOpSrc.get()!=NULL

• this->isCompatible(*fwdOpSrc)==true

• Op!=NULL

• *Op must have been created by this->createOp() prior to calling this function.

• [supportSolveUse==SUPPORT_SOLVE_FORWARD_ONLY] this->solveSupportsConj(conj)==true for any value of conj

• [supportSolveUse==SUPPORT_SOLVE_TRANSPOSE_ONLY] this->solveTransposeSupportsConj(conj)==true for any value of conj

• [supportSolveUse==SUPPORT_SOLVE_FORWARD_AND_TRANSPOSE] this->solveSupportsConj(conj)==true && this->solveTransposeSupportsConj(conj)==true for any value of conj

Postconditions:

• Throws CatastrophicSolveFailure if the underlying linear solver could not be created successfully (e.g. due to a factorization failure or some other cause).

• Op->range()->isCompatible(*fwdOpSrc->range())==true

• Op->domain()->isCompatible(*fwdOpSrc->domain())==true

• Op->apply() and Op->applyTranspose() must behave exactly the same as fwdOpSrc->apply() and fwdOpSrc->applyTranspose()

• Op->solveSupportsConj(conj)==this->solveSupportsConj(conj)

• Op->solveTransposeSupportsConj(conj)==this->solveTransposeSupportsConj(conj)

• fwdOpSrc.count() after output is greater than fwdOpSrc.count() just before this call and therefore the client can assume that the *fwdOpSrc object will be remembered by the *Op object. The client must be careful not to modify the *fwdOpSrc object or else the *Op object may also be modified and become invalid.

template<class Scalar>

void Thyra::LinearOpWithSolveFactoryBase< Scalar >::initializeAndReuseOp	(	const RCP< const LinearOpSourceBase< Scalar > > &	fwdOpSrc,
		LinearOpWithSolveBase< Scalar > *	Op
	)			const [virtual]

Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object but allow for reuse of any preprocessing that is in *Op.

Parameters:

	fwdOpSrc	[in] The forward linear operator that will be used to create the output LinearOpWithSolveBase object.
	Op	[in/out] The output LinearOpWithSolveBase object. This object must have be created first by this->createOp() and may have already been through at least one previous set of calls to this->initializeOp() and this->uninitializeOp(). Note that subclasses should always first strip off the transpose and scaling by calling unwrap() before attempting to dynamic cast the object.

Preconditions:

• fwdOpSrc.get()!=NULL

• this->isCompatible(*fwdOpSrc)==true

• Op!=NULL

• *Op must have been created by this->createOp() prior to calling this function.

Postconditions:

• Throws CatastrophicSolveFailure if the underlying linear solver could not be created successfully (e.g. due to a factorization failure or some other cause).

• Op->range()->isCompatible(*fwdOpSrc->range())==true

• Op->domain()->isCompatible(*fwdOpSrc->domain())==true

• Op->apply() and Op->applyTranspose() must behave exactly the same as fwdOpSrc->apply() and fwdOpSrc->applyTranspose()

• Op->solveSupportsConj(conj)==this->solveSupportsConj(conj)

• Op->solveTransposeSupportsConj(conj)==this->solveTransposeSupportsConj(conj)

• fwdOpSrc.count() after output is greater than fwdOpSrc.count() just before this call and therefore the client can assume that the *fwdOpSrc object will be remembered by the *Op object. The client must be careful not to modify the *fwdOpSrc object or else the *Op object may also be modified.

The purpose of this function is to allow the reuse of old factorizations and/or preconditioners that may go into the initialization of the *Op objects. Note that by calling this function, the performance Op->solve(...) may not be as good as when calling the function this->initializeOp(...,Op) to initialize *Op.

The default implementation of this function just calls this->initializeOp(fwdOpSrc,Op) which does the default implementation.

Definition at line 80 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

virtual void Thyra::LinearOpWithSolveFactoryBase< Scalar >::uninitializeOp	(	LinearOpWithSolveBase< Scalar > *	Op,
		RCP< const LinearOpSourceBase< Scalar > > *	fwdOpSrc = NULL,
		RCP< const PreconditionerBase< Scalar > > *	prec = NULL,
		RCP< const LinearOpSourceBase< Scalar > > *	approxFwdOpSrc = NULL,
		ESupportSolveUse *	supportSolveUse = NULL
	)			const [pure virtual]

Uninitialize a LinearOpWithSolveBase object and return its remembered forward linear operator and potentially also its externally generated preconditioner.

Parameters:

	Op	[in/out] On input, *Op is an initialized or uninitialized object and on output is uninitialized. Note that "uninitialized" does not mean that Op is completely stateless. It may still remember some aspect of the matrix fwdOpSrc that will allow for a more efficient initialization next time through this->initializeOp().
	fwdOpSrc	[in/out] If fwdOpSrc!=NULL on input, then on output this is set to the same forward operator passed into this->initializeOp().
	prec	[in/out] If prep!=NULL on input, then on output, this this is set to same preconditioner that was passed into this->initializePreconditionedOp().
	approxFwdOpSrc	[in/out] If approxFwdOpSrc!=NULL on input, then on output, this is set to same approximate forward operator that was passed into this->initializePreconditionedOp().
	ESupportSolveUse	[in/out] If fwdOpSrc!=NULL on input, then on output this is set to same option value passed to this->initializeOp().

Preconditions:

• *Op must have been created by this->createOp() prior to calling this function.

• Op may or may not have been passed through a call to this->initializeOp() or this->initializePreconditionedOp().

Postconditions:

• If *Op on input was initialized through a call to this->initializeOp() and if fwdOpSrc!=NULL then (*fwdOpSrc).get()!=NULL.

• If *Op was uninitialized on input and fwdOpSrc!=NULL then fwdOpSrc->get()==NULL out output.

• On output, *Op can be considered to be uninitialized and it is safe to modify the forward operator object <tt>*(*fwdOpSrc) returned in fwdOpSrc. The default is fwdOpSrc==NULL in which case the forward operator will not be returned in *fwdOpSrc.

This function should be called before the forward operator passed in to this->initializeOp() is modified. Otherwise, *this could be left in an inconsistent state. However, this is not required.

template<class Scalar>

bool Thyra::LinearOpWithSolveFactoryBase< Scalar >::supportsPreconditionerInputType

(

const EPreconditionerInputType

precOpType

)

const [virtual]

Determines if *this supports given preconditioner type.

The default implementation returns false.

Definition at line 90 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

void Thyra::LinearOpWithSolveFactoryBase< Scalar >::initializePreconditionedOp	(	const RCP< const LinearOpSourceBase< Scalar > > &	fwdOpSrc,
		const RCP< const PreconditionerBase< Scalar > > &	prec,
		LinearOpWithSolveBase< Scalar > *	Op,
		const ESupportSolveUse	supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED
	)			const [virtual]

Initialize a pre-created LinearOpWithSolveBase object given a "compatible" LinearOpBase object and an optional PreconditionerBase object.

Parameters:

	fwdOpSrc	[in] The forward linear operator that will be used to create the output LinearOpWithSolveBase object.
	prec	[in] The preconditioner that will be used to create the output LinearOpWithSolveBase object if preconditioners are supported.
	Op	[in/out] The output LinearOpWithSolveBase object. This object must have be created first by this->createOp(). The object may have also already been passed through this function several times. Note that subclasses should always first strip off the transpose and scaling by calling unwrap() before attempting to dynamic cast the object.
	supportSolveUse	[in] Determines if Op->solve(...) or Op->solveTranspose(...) will be called. This allows *this factory object determine how to best initialize the *Op object. Default supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED

Preconditions:

• fwdOpSrc.get()!=NULL

• prec.get()!=NULL

• this->isCompatible(*fwdOpSrc)==true

• Op!=NULL

• *Op must have been created by this->createOp() prior to calling this function.

• It is allowed for an implementation to throw an exception if this->supportsPreconditionerInputType(PRECONDITIONER_INPUT_TYPE_AS_OPERATOR)==false but this is not required.

• [supportSolveUse==SUPPORT_SOLVE_FORWARD_ONLY] this->solveSupportsConj(conj)==true for any value of conj

• [supportSolveUse==SUPPORT_SOLVE_TRANSPOSE_ONLY] this->solveTransposeSupportsConj(conj)==true for any value of conj

• [supportSolveUse==SUPPORT_SOLVE_FORWARD_AND_TRANSPOSE] this->solveSupportsConj(conj)==true && this->solveTransposeSupportsConj(conj)==true for any value of conj

Postconditions:

• Throws CatastrophicSolveFailure if the underlying linear solver could not be created successfully (e.g. due to a factorization failure or some other cause).

• Op->range()->isCompatible(*fwdOpSrc->range())==true

• Op->domain()->isCompatible(*fwdOpSrc->domain())==true

• Op->apply() and Op->applyTranspose() must behave exactly the same as fwdOpSrc->apply() and fwdOpSrc->applyTranspose()

• <t>Op->solveSupportsConj(conj)==this->solveSupportsConj(conj)

• <t>Op->solveTransposeSupportsConj(conj)==this->solveTransposeSupportsConj(conj)

• fwdOpSrc.count() after output is greater than fwdOpSrc.count() just before this call and therefore the client can assume that the *fwdOpSrc object will be remembered by the *Op object. The client must be careful not to modify the *fwdOpSrc object or else the *Op object may also be modified.

• If this->supportsPreconditionerInputType(PRECONDITIONER_INPUT_TYPE_AS_OPERATOR)==true then
  • prec.count() after output is greater than prec.count() just before this call and therefore the client can assume that the *prec object will be remembered by the *Op object. The client must be careful not to modify the *prec object or else the *Op object may also be modified.

• else if an exception is not thrown then
  • prec.count() after output is equal to prec.count() just before this call and therefore the *prec is ignored and is not remembered.

Warning! It is allowed for an implementation to throw an exception if this->supportsPreconditionerInputType(PRECONDITIONER_INPUT_TYPE_AS_OPERATOR)==false so therefore a client should not call this function if preconditioners are not supported! The mode of silently ignoring preconditioners is acceptable in some cases and is therefore allowed behavior.

The default implementation throws an exception which is consistent with the default implementation of this->supportsPreconditionerInputType() which assumes by default that preconditioners can not be supported..

Definition at line 99 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

template<class Scalar>

void Thyra::LinearOpWithSolveFactoryBase< Scalar >::initializeApproxPreconditionedOp	(	const RCP< const LinearOpSourceBase< Scalar > > &	fwdOpSrc,
		const RCP< const LinearOpSourceBase< Scalar > > &	approxFwdOpSrc,
		LinearOpWithSolveBase< Scalar > *	Op,
		const ESupportSolveUse	supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED
	)			const [virtual]

Initialize a pre-created LinearOpWithSolveBase object given a "compatible" forward LinearOpBase object and an approximate forward LinearOpBase object.

Parameters:

	fwdOpSrc	[in] The forward linear operator that will be used to create the output LinearOpWithSolveBase object.
	approxFwdOpSrc	[in] Approximation to fwdOpSrc from which a preconditioner will be created for.
	Op	[in/out] The output LinearOpWithSolveBase object. This object must have be created first by this->createOp(). The object may have also already been passed through this function several times. Note that subclasses should always first strip off the transpose and scaling by calling unwrap() before attempting to dynamic cast the object.
	supportSolveUse	[in] Determines if Op->solve(...) or Op->solveTranspose(...) will be called. This allows *this factory object determine how to best initialize the *Op object. Default supportSolveUse=SUPPORT_SOLVE_UNSPECIFIED

ToDo: finish documentation!

Definition at line 115 of file Thyra_LinearOpWithSolveFactoryBase_def.hpp.

---

The documentation for this class was generated from the following files:

• Thyra_LinearOpWithSolveFactoryBase_decl.hpp
• Thyra_LinearOpWithSolveFactoryBase_def.hpp

---