vector interfaces]

Base class for all linear operators. More...

#include <Thyra_LinearOpBaseDecl.hpp>

Inheritance diagram for Thyra::LinearOpBase< RangeScalar, DomainScalar >:


Public pure virtual functions (must be overridden by subclass)
virtual RCP< const VectorSpaceBase< RangeScalar > >	range () const =0
	Return a smart pointer for the range space for `this` operator.
virtual RCP< const VectorSpaceBase< DomainScalar > >	domain () const =0
	Return a smart pointer for the domain space for `this` operator.
virtual void	apply (const EConj conj, const MultiVectorBase< DomainScalar > &X, MultiVectorBase< RangeScalar > *Y, const RangeScalar alpha=Teuchos::ScalarTraits< RangeScalar >::one(), const RangeScalar beta=Teuchos::ScalarTraits< RangeScalar >::zero()) const =0
	Apply the forward non-conjugate or conjugate linear operator to a multi-vector : `Y = alphaMX + beta*Y`.
Public virtual functions with default implementations
virtual bool	applySupports (const EConj conj) const
	Determines if `apply()` supports this `conj` argument.
virtual bool	applyTransposeSupports (const EConj conj) const
	Determines if `applyTranspose()` supports this `conj` argument.
virtual void	applyTranspose (const EConj conj, const MultiVectorBase< RangeScalar > &X, MultiVectorBase< DomainScalar > *Y, const DomainScalar alpha=Teuchos::ScalarTraits< DomainScalar >::one(), const DomainScalar beta=Teuchos::ScalarTraits< DomainScalar >::zero()) const
	Apply the non-conjugate or conjugate transposed linear operator to a multi-vector : `Y = alphatrans(M)X + beta*Y`.
virtual RCP< const LinearOpBase< RangeScalar, DomainScalar > >	clone () const
	Clone the linear operator object (if supported).
Related Functions
(Note that these are not member functions.)
bool	opSupported (const LinearOpBase< Scalar > &M, ETransp M_trans)
	Determines if an operation is supported for a single scalar type.
void	apply (const LinearOpBase< RangeScalar, DomainScalar > &M, const EConj conj, const MultiVectorBase< DomainScalar > &X, MultiVectorBase< RangeScalar > *Y, const RangeScalar alpha=ScalarTraits< RangeScalar >::one(), const RangeScalar beta=ScalarTraits< RangeScalar >::zero())
	Call `LinearOpBase::apply()` as a global function call.
void	applyTranspose (const LinearOpBase< RangeScalar, DomainScalar > &M, const EConj conj, const MultiVectorBase< RangeScalar > &X, MultiVectorBase< DomainScalar > *Y, const DomainScalar alpha=ScalarTraits< DomainScalar >::one(), const DomainScalar beta=ScalarTraits< DomainScalar >::zero())
	Call `LinearOpBase::applyTranspose()` as a global function call.
void	apply (const LinearOpBase< Scalar > &M, const ETransp M_trans, const MultiVectorBase< Scalar > &X, MultiVectorBase< Scalar > *Y, const Scalar alpha=ScalarTraits< Scalar >::one(), const Scalar beta=ScalarTraits< Scalar >::zero())
	Call `LinearOpBase::apply()` or `LinearOpBase::applyTranspose()` as a global function call (for a single scalar type).

Detailed Description

template<class RangeScalar, class DomainScalar>
class Thyra::LinearOpBase< RangeScalar, DomainScalar >

Base class for all linear operators.

Introduction

A linear operator can optionally perform the forward operations

Forward non-conjugate apply
```
Y = alpha*M*X + beta*Y 
```
Forward conjugate apply
```
Y = alpha*conjugate(M)*X + beta*Y 
```

through the apply() function and the operations

Transpose non-conjugate apply
```
Y = alpha*transpose(M)*X + beta*Y 
```
Transpose conjugate (i.e. adjoint) apply
```
Y = alpha*adjoint(M)*X + beta*Y 
```

through the applyTranspose() function where Y and X are MultiVectorBase objects. The reason for the exact form of the above operations is that there are direct BLAS and equivalent versions of these operations and performing a sum-into multiplication is more efficient in general.

Range and domain spaces

A linear operator has vector spaces associated with it for the vectors x and y that lie in the domain and the range spaces of the non-transposed linear operator y = M*x and these spaces are returned by domain() and range().

Support for different range and domain scalar types

This interface allows for different scalar types for the range and domain spaces, RangeScalar and DomainScalar respectively. This is needed to support such things as real-to-complex FFTs (see the example RealToComplex1DFFTLinearOp for instance) and real to extended-precision linear operators.

Single scalar type linear operators

While this interface supports the notion of different range and domain scalar types, many different LinearOpBase implementations and many different ANAs will only support a single scalar type. There is a lot of support for single-scalar-type linear operators.

First, note that there is a forward declaration for this class of the form

  template<class RangeScalar, class DomainScalar = RangeScalar> class LinearOpBase;

that allows the class to be refereed to using just one Scalar type as LinearOpBase<Scalar>. This is useful for both clients and subclass implementations.

When a client ANA can only support a single-scalar type of linear operator, it may be more convenient to use some of the wrapper functions such as Thyra::apply() and Thyra::opSupported() that are described here.

Scalar products and the adjoint relation

Note that the vector spaces returned from domain() and range() may have specialized implementations of the scalar product $<u,w>$

(i.e. $<u,w> \neq u^H w$ in general). As a result, the operator and adjoint operator must obey the defined scalar products. Specifically, for any two vectors $w\in\mathcal{D}$ (in the domain space of A) and $u\in\mathcal{R}$ (in the range space of A), the adjoint operation must obey the adjoint property

$<u,A v>_{\mathcal{R}} =\!= <A^H u, v>_{\mathcal{D}}$

where $<.,.>_{\mathcal{R}}$ is the scalar product defined by this->range()->scalarProd() and $<.,.>_{\mathcal{D}}$ is the scalar product defined by this->domain()->scalarProd(). This property of the adjoint can be checked numerically, if adjoints are supported, using the testing utility class LinearOpTester.

Aliasing policy

It is strictly forbidden to alias the input/output object Y with the input object X in apply() or applyTranspose(). Allowing aliasing would greatly complicate the development of concrete subclasses.

Optional support for specific types of operator applications

This interface does not require that a linear operator implementation support all of the different types of operator applications defined in the introduction above. If a LinearOpBase object can not support a particular type of operator application, then this is determined by the functions applySupports() and applyTransposeSupports().

Testing LinearOpBase objects

The concrete class LinearOpTester provides a full featured set of tests for any LinearOpBase object. This testing class can check if the operator is truly "linear", and/or if the adjoint relationship holds, and/or if an operator is symmetric. All of the tests are controlled by the client, can be turned on and off, and pass/failure is determined by tolerances that the client can specify. In addition, this testing class can also check if two linear operators are approximately equal.

Notes for subclass developers

There are only three functions that a concrete subclass is required to override: domain(), range() and apply(). Note that the functions domain() and range() should simply return VectorSpaceBase objects for subclasses that are already defined for the vectors that the linear operator interacts with through the function apply(). Therefore, given that appropriate VectorSpaceBase and MultiVectorBase (and/or VectorBase) subclasses exist, the only real work involved in implementing a LinearOpBase subclass is in defining a single function apply().

This interface provides default implementations for the functions applyTranspose() and applyTransposeSupports() where it is assumed that the operator does not support transpose (or adjoint) operator applications. If transpose (and/or adjoint) operator applications can be supported, then the functions applyTranspose() and applyTransposeSupports() should be overridden as well.

If possible, the subclass should also override the clone() function which allows clients to create copies of a LinearOpBase object. This functionality is useful in some circumstances. However, this functionality is not required and the default clone() implementation returns a null smart pointer object.

If a concrete subclass can only support a single scalar type, then the concrete subclass should perhaps inherit form the SingleScalarLinearOpBase node subclass. This node subclass provides just a single function SingleScalarLinearOpBase::apply() that will support both forward and transpose (adjoint) operator applications.

If, in addition to only supporting a single scalar type, a concrete subclass can only support single RHS operator applications, then perhaps the node subclass SingleRhsLinearOpBase should be inherited from. This node subclass provides a version of SingleRhsLinearOpBase::apply() that takes VectorBase objects instead of MultiVectorBase objects.

Definition at line 202 of file Thyra_LinearOpBaseDecl.hpp.

Member Function Documentation

template<class RangeScalar, class DomainScalar>

virtual RCP< const VectorSpaceBase<RangeScalar> > Thyra::LinearOpBase< RangeScalar, DomainScalar >::range ( ) const [pure virtual]

Return a smart pointer for the range space for this operator.

Note that a return value of return.get()==NULL is a flag that *this is not fully initialized.

If return.get()!=NULL, it is required that the object referenced by *return.get() must have lifetime that extends past the lifetime of the returned smart pointer object. However, the object referenced by *return.get() may change if *this modified so this reference should not be maintained for too long.

New Behavior! It is required that the VectorSpaceBase object embedded in return must be valid past the lifetime of *this linear operator object.

template<class RangeScalar, class DomainScalar>

virtual RCP< const VectorSpaceBase<DomainScalar> > Thyra::LinearOpBase< RangeScalar, DomainScalar >::domain ( ) const [pure virtual]

Return a smart pointer for the domain space for this operator.

Note that a return value of return.get()==NULL is a flag that *this is not fully initialized.

New Behavior! It is required that the VectorSpaceBase object embedded in return must be valid past the lifetime of *this linear operator object.

template<class RangeScalar, class DomainScalar>

virtual void Thyra::LinearOpBase< RangeScalar, DomainScalar >::apply	(	const EConj	conj,
		const MultiVectorBase< DomainScalar > &	X,
		MultiVectorBase< RangeScalar > *	Y,
		const RangeScalar	alpha = `Teuchos::ScalarTraits< RangeScalar >::one()`,
		const RangeScalar	beta = `Teuchos::ScalarTraits< RangeScalar >::zero()`
	)			const `[pure virtual]`

Apply the forward non-conjugate or conjugate linear operator to a multi-vector : Y = alpha*M*X + beta*Y.

Parameters:

	conj	[in] Determines whether the elements are non-conjugate or conjugate. The value `NONCONJ_ELE` gives the standard forward operator while the value of `CONJ_ELE` gives the forward operator with the complex conjugate matrix elements. For a real-valued operators, this argument is ignored and has no effect.
	X	[in] The right hand side multi-vector
	Y	[in/out] The target multi-vector being transformed
	alpha	[in] Scalar multiplying `M`, where `M==*this`. The default value of `alpha` is 1.0
	beta	[in] The multiplier for the target multi-vector `Y`. The default value of `beta` is `0.0`.

Preconditions:

this->applySupports(conj)==true (throw Exceptions::OpNotSupported)
this->domain().get()!=NULL && this->range().get()!=NULL (throw std::logic_error)
X.range()->isCompatible(this->domain()) == true (throw Exceptions::IncompatibleVectorSpaces)
Y->range()->isCompatible(*this->range()) == true (throw Exceptions::IncompatibleVectorSpaces)
Y->domain()->isCompatible(*X.domain()) == true (throw Exceptions::IncompatibleVectorSpaces)
Y can not alias X. It is up to the client to ensure that Y and X are distinct since in general this can not be verified by the implementation until, perhaps, it is too late. If possible, an exception will be thrown if aliasing is detected.

Postconditions:

Is it not obvious? After the function returns the multi-vector Y is transformed as indicated above.

template<class RangeScalar, class DomainScalar>

bool Thyra::LinearOpBase< RangeScalar, DomainScalar >::applySupports ( const EConj conj ) const [virtual]

Determines if apply() supports this conj argument.

The default implementation returns true for real valued scalar types or when conj==NONCONJ_ELE for complex valued types.

Definition at line 41 of file Thyra_LinearOpBase.hpp.

template<class RangeScalar, class DomainScalar>

bool Thyra::LinearOpBase< RangeScalar, DomainScalar >::applyTransposeSupports ( const EConj conj ) const [virtual]

Determines if applyTranspose() supports this conj argument.

The default implementation returns false which is consistent with the below default implementation for applyTranspose().

Definition at line 51 of file Thyra_LinearOpBase.hpp.

template<class RangeScalar, class DomainScalar>

void Thyra::LinearOpBase< RangeScalar, DomainScalar >::applyTranspose	(	const EConj	conj,
		const MultiVectorBase< RangeScalar > &	X,
		MultiVectorBase< DomainScalar > *	Y,
		const DomainScalar	alpha = `Teuchos::ScalarTraits< DomainScalar >::one()`,
		const DomainScalar	beta = `Teuchos::ScalarTraits< DomainScalar >::zero()`
	)			const `[virtual]`

Apply the non-conjugate or conjugate transposed linear operator to a multi-vector : Y = alpha*trans(M)*X + beta*Y.

Parameters:

	conj	[in] Determines whether the elements are non-conjugate or conjugate. The value `NONCONJ_ELE` gives the standard transposed operator with non-conjugate transposed matrix elements while the value of `CONJ_ELE` gives the standard adjoint operator with the complex conjugate transposed matrix elements. For a real-valued operators, this argument is ignored and has no effect.
	X	[in] The right hand side multi-vector
	Y	[in/out] The target multi-vector being transformed
	alpha	[in] Scalar multiplying `M`, where `M==*this`. The default value of `alpha` is 1.0
	beta	[in] The multiplier for the target multi-vector `Y`. The default value of `beta` is `0.0`.

Preconditions:

this->applyTransposeSupports(conj)==true (throw Exceptions::OpNotSupported)
this->domain().get()!=NULL && this->range().get()!=NULL (throw std::logic_error)
X.range()->isCompatible(this->range()) == true (throw Exceptions::IncompatibleVectorSpaces)
Y->range()->isCompatible(*this->domain()) == true (throw Exceptions::IncompatibleVectorSpaces)
Y->domain()->isCompatible(*X.domain()) == true (throw Exceptions::IncompatibleVectorSpaces) and X are distinct since in general this can not be verified by the implementation until, perhaps, it is too late. If possible, an exception will be thrown if aliasing is detected.
Y can not alias X. It is up to the client to ensure that Y

Postconditions:

Is it not obvious? After the function returns the multi-vector Y is transformed as indicated above.

The default implementation throws an exception but gives a very good error message. The assumption here is that most linear operators will not be able to support an transpose apply and that is why this default implementation is provided.

Definition at line 59 of file Thyra_LinearOpBase.hpp.

template<class RangeScalar, class DomainScalar>

RCP< const LinearOpBase< RangeScalar, DomainScalar > > Thyra::LinearOpBase< RangeScalar, DomainScalar >::clone ( ) const [virtual]

Clone the linear operator object (if supported).

The primary purpose for this function is to allow a client to capture the current state of a linear operator object and be guaranteed that some other client will not alter its behavior. A smart implementation will use reference counting and lazy evaluation internally and will not actually copy any large amount of data unless it has to.

The default implementation returns return.get()==NULL which is allowable. A linear operator object is not required to return a non-NULL value but many good matrix-based linear operator implementations will.

Reimplemented in Thyra::MultiVectorBase< Scalar >.

Definition at line 80 of file Thyra_LinearOpBase.hpp.

Friends And Related Function Documentation

template<class RangeScalar, class DomainScalar>

bool opSupported	(	const LinearOpBase< Scalar > &	M,
		ETransp	M_trans
	)			`[related]`

Determines if an operation is supported for a single scalar type.

Definition at line 378 of file Thyra_LinearOpBaseDecl.hpp.

template<class RangeScalar, class DomainScalar>

void apply	(	const LinearOpBase< RangeScalar, DomainScalar > &	M,
		const EConj	conj,
		const MultiVectorBase< DomainScalar > &	X,
		MultiVectorBase< RangeScalar > *	Y,
		const RangeScalar	alpha = `ScalarTraits<RangeScalar>::one()`,
		const RangeScalar	beta = `ScalarTraits<RangeScalar>::zero()`
	)			`[related]`

Call LinearOpBase::apply() as a global function call.

Calls M.apply(conj,X,Y,alpha,beta).

Definition at line 392 of file Thyra_LinearOpBaseDecl.hpp.

template<class RangeScalar, class DomainScalar>

void applyTranspose	(	const LinearOpBase< RangeScalar, DomainScalar > &	M,
		const EConj	conj,
		const MultiVectorBase< RangeScalar > &	X,
		MultiVectorBase< DomainScalar > *	Y,
		const DomainScalar	alpha = `ScalarTraits<DomainScalar>::one()`,
		const DomainScalar	beta = `ScalarTraits<DomainScalar>::zero()`
	)			`[related]`

Call LinearOpBase::applyTranspose() as a global function call.

Calls M.applyTranspose(conj,X,Y,alpha,beta).

Definition at line 412 of file Thyra_LinearOpBaseDecl.hpp.

template<class RangeScalar, class DomainScalar>

void apply	(	const LinearOpBase< Scalar > &	M,
		const ETransp	M_trans,
		const MultiVectorBase< Scalar > &	X,
		MultiVectorBase< Scalar > *	Y,
		const Scalar	alpha = `ScalarTraits<Scalar>::one()`,
		const Scalar	beta = `ScalarTraits<Scalar>::zero()`
	)			`[related]`

Call LinearOpBase::apply() or LinearOpBase::applyTranspose() as a global function call (for a single scalar type).

Calls M.apply(...,X,Y,alpha,beta) or M.applyTranspose(...,X,Y,alpha,beta).

Definition at line 435 of file Thyra_LinearOpBaseDecl.hpp.

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

Generated on Tue Oct 20 12:46:45 2009 for Fundamental Thyra ANA Operator/Vector Interfaces by

1.4.7

Thyra::LinearOpBase< RangeScalar, DomainScalar > Class Template Reference [C++ code for foundational Thyra operator/vector interfaces]

Public pure virtual functions (must be overridden by subclass)

Public virtual functions with default implementations

Related Functions

Detailed Description

template<class RangeScalar, class DomainScalar> class Thyra::LinearOpBase< RangeScalar, DomainScalar >

Member Function Documentation

Friends And Related Function Documentation

Thyra::LinearOpBase< RangeScalar, DomainScalar > Class Template Reference
[C++ code for foundational Thyra operator/vector interfaces]

template<class RangeScalar, class DomainScalar>
class Thyra::LinearOpBase< RangeScalar, DomainScalar >