Implements Epetra_Operator as a product of one or more Epetra_Operator objects. More...

#include <EpetraExt_ProductOperator.h>

Inheritance diagram for EpetraExt::ProductOperator:

Public types
enum	EApplyMode { APPLY_MODE_APPLY, APPLY_MODE_APPLY_INVERSE }

Constructors / initializers / accessors
	ProductOperator ()
	Construct to uninitialized. More...

	ProductOperator (const int num_Op, const Teuchos::RCP< const Epetra_Operator > Op[], const Teuchos::ETransp Op_trans[], const EApplyMode Op_inverse[])
	Calls `initialize()`. More...

void	initialize (const int num_Op, const Teuchos::RCP< const Epetra_Operator > Op[], const Teuchos::ETransp Op_trans[], const EApplyMode Op_inverse[])
	Setup with constituent operators. More...

void	uninitialize (int *num_Op, Teuchos::RCP< const Epetra_Operator > Op[], Teuchos::ETransp Op_trans[], EApplyMode p_inverse[])
	Set to an uninitialized state and wipe out memory. More...

void	applyConstituent (const int k, Teuchos::ETransp Op_trans, EApplyMode Op_inverse, const Epetra_MultiVector &X_k, Epetra_MultiVector *Y_k) const
	Apply the kth aggregate operator M[k] correctly. More...

int	num_Op () const
	Return the number of aggregate opeators. More...

Teuchos::RCP< const Epetra_Operator >	Op (int k) const
	Access the kth operator (zero-based). More...

Teuchos::ETransp	Op_trans (int k) const
	Access the transpose mode of the kth operator (zero-based). More...

EApplyMode	Op_inverse (int k) const
	Access the inverse mode of the kth operator (zero-based). More...

Overridden from Epetra_Operator
int	SetUseTranspose (bool UseTranspose)

int	Apply (const Epetra_MultiVector &X, Epetra_MultiVector &Y) const

int	ApplyInverse (const Epetra_MultiVector &X, Epetra_MultiVector &Y) const

double	NormInf () const

const char *	Label () const

bool	UseTranspose () const

bool	HasNormInf () const

const Epetra_Comm &	Comm () const

const Epetra_Map &	OperatorDomainMap () const

const Epetra_Map &	OperatorRangeMap () const

Additional Inherited Members

Detailed Description

Implements Epetra_Operator as a product of one or more Epetra_Operator objects.

This class implements a product operator of the form:

M = M[0]*M[1]*...*M[num_Op-1]

and operator applications are performed one constituent operator at a time as:

Forward Mat-vec: Y = M * X

  T[k-1] = M[k]*T[k]

    for k = num_Op-1...0

      where: T[num_Op-1] = X (input vector)
      where: T[-1]       = Y (output vector)

Adjoint Mat-vec: Y = M' * X

  T[k] = M[k]'*T[k-1]

    for k = 0...num_Op-1

      where: T[-1]       = X (input vector)
      where: T[num_Op-1] = Y (output vector)

Likewise, the inverse can also be applied (if all of the constituent operators support the inverse operation) as:

Forward Inverse Mat-vec: Y = inv(M) * X

  T[k] = inv(M[k])*T[k-1]

    for k = 0...num_Op-1

    for k = 0...num_Op-1

      where: T[-1]       = X (input vector)
      where: T[num_Op-1] = Y (output vector)

Adjoint Inverse Mat-vec: Y = inv(M') * X

  T[k] = inv(M[k]')*T[k-1]

    for k = num_Op-1...0

      where: T[num_Op-1] = X (input vector)
      where: T[-1]       = Y (output vector)

Note that maps for the result of the inverse of an operator is the same as the result of the adjoint of the operator and the map for the result of the inverse of the adjoint is the same as for the result the non-inverse forward opeator.

The client constructs this object with a list of Epetra_Operator objects an how the non-transposed operator is to be viewed and if it is to be views as its inverse or not (see initialize()).

Note: The Epetra_Map objects returned from OperatorDomainMap() and OperatorRangeMap() must always be with respect to the non-transposed operator! This is very strange behavior and is totally undocumented in the Epetra_Operator interface but it seems to be the case.

Definition at line 139 of file EpetraExt_ProductOperator.h.

Member Enumeration Documentation

enum EpetraExt::ProductOperator::EApplyMode

Enumerator
APPLY_MODE_APPLY
APPLY_MODE_APPLY_INVERSE

Definition at line 146 of file EpetraExt_ProductOperator.h.

Constructor & Destructor Documentation

EpetraExt::ProductOperator::ProductOperator ( )

Construct to uninitialized.

Definition at line 50 of file EpetraExt_ProductOperator.cpp.

EpetraExt::ProductOperator::ProductOperator	(	const int	num_Op,
		const Teuchos::RCP< const Epetra_Operator >	Op[],
		const Teuchos::ETransp	Op_trans[],
		const EApplyMode	Op_inverse[]
	)

Calls initialize().

Definition at line 53 of file EpetraExt_ProductOperator.cpp.

Member Function Documentation

void EpetraExt::ProductOperator::initialize	(	const int	num_Op,
		const Teuchos::RCP< const Epetra_Operator >	Op[],
		const Teuchos::ETransp	Op_trans[],
		const EApplyMode	Op_inverse[]
	)

Setup with constituent operators.

    @param  num_Op  [in] Number of constinuent operators.
    @param  Op      [in] Array (length <tt>num_Op</tt>) of smart pointers
                    to <tt>Epetra_Operator objects that implement each
                    constituent operator.
    @param  Op_trans
                    [in] Array (length <tt>num_Op</tt>) that determines
                    the transpose mode of the constituent operator (see below).
    @param  Op_inverse
                    [in] Array (length <tt>num_Op</tt>) that determines
                    the inverse apply mode of the constituent operator (see below).

    Preconditions:<ul>
    <li><tt>num_Op > 0</tt>
    <li><tt>Op!=NULL</tt>
    <li><tt>Op[k].get()!=NULL</tt>, for <tt>k=0...num_Op-1</tt>
    <li><tt>Op_trans!=NULL</tt>
    <li><tt>Op_inverse!=NULL</tt>
    </ul>

    Postconditions:<ul>
    <li><tt>this->num_Op()==num_Op</tt>
    <li><tt>this->Op(k).get()==Op[k].get()</tt>, for <tt>k=0...num_Op-1</tt>
    <li><tt>this->Op_trans(k)==Op_trans[k]</tt>, for <tt>k=0...num_Op-1</tt>
    <li><tt>this->Op_inverse(k)==Op_inverse[k]</tt>, for <tt>k=0...num_Op-1</tt>
    </ul>

    The forward constituent operator <tt>T[k-1] = M[k]*T[k]</tt> described in the
    main documenatation above is defined as follows:

Op[k]->SetUseTranspose( Op_trans[k]!=Teuchos::NO_TRANS );
if( Op_inverse[k]==APPLY_MODE_APPLY )
  Op[k]->Apply( T[k], T[k-1] );
      else
  Op[k]->ApplyInverse( T[k], T[k-1] );

The inverse constituent operator T[k] = inv(M[k])*T[k-1] described in the main documenatation above is defined as follows:

Op[k]->SetUseTranspose( Op_trans[k]!=Teuchos::NO_TRANS );
if( Op_inverse[k]==APPLY_MODE_APPLY )
  Op[k]->ApplyInverse( T[k-1], T[k] );
      else
  Op[k]->Apply( T[k-1], T[k] );

The other transposed constituent operators M[k]' and inv(M[k]') are defined by simply changing the value of the transpose as Op[k]->SetUseTranspose( Op_trans[k]==Teuchos::NO_TRANS );. Note, that Op[k]->SetUseTranspose(...) is called immediately before Op[k]->Apply(...) or Op[k]->ApplyInverse(...) is called to avoid tricky problems that could occur with multiple uses of the same operator.

Definition at line 63 of file EpetraExt_ProductOperator.cpp.

void EpetraExt::ProductOperator::uninitialize	(	int *	num_Op,
		Teuchos::RCP< const Epetra_Operator >	Op[],
		Teuchos::ETransp	Op_trans[],
		EApplyMode	p_inverse[]
	)

Set to an uninitialized state and wipe out memory.

Postconditions:

this->num_Op()==0

Definition at line 89 of file EpetraExt_ProductOperator.cpp.

void EpetraExt::ProductOperator::applyConstituent	(	const int	k,
		Teuchos::ETransp	Op_trans,
		EApplyMode	Op_inverse,
		const Epetra_MultiVector &	X_k,
		Epetra_MultiVector *	Y_k
	)		const

Apply the kth aggregate operator M[k] correctly.

Parameters

k	[in] Gives the index (zero-based) of the constituent operator to apply.
Op_trans	[in] Determines if the transpose of the constituent operator is to be applied.
Op_inverse	[in] Determines if the operator or its inverse (if supported) should be applied.
X_k	[in] The input vector.
Y_k	[out] The output vector.