AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects: AbstractLinAlgPack::DirectSparseSolver Class Reference

This interface is designed to accommodate both unsymmetic and symmetric direct solvers. The motivation for using one common interface is that we would like to be able to easily use an unsymmetic solver on a symmetric system. Of couse this would also allow a client to attempt to solve an unsymmetric system with a symmetric solver which of course will not work. A symmetric implementation can of course check that A.rows()==A.cols() but it would be much more difficult to determine if the matrix was indeed symmetric. It is likely that the symmetric solver would only extract the upper or lower triangular region of an unsymmetric matrix and would never check the other triangular region for compatibility.

The interface is designed to allow for maximum memory sharing and recycling. For example, we would like for the client to be able to maintain more than one set of nonzero factors using the same sparsity and fill-in structure determined by the analyze_and_factor() method.

The basic idea is that given an m x n rectangular matrix A, where m <= n, objects of this type will find row and column permutations P and Q and a rank r such that C = (P'*A*Q)(1:r,1:r) is the largest square nonsingular (well conditioned) matrix that can be found. This basis matrix C is represented as a BasisMatrix object that supports the MatrixNonsingSerial interface.

All the information needed to solve for a linear system and to factor other matrices with the same structure is contained in the BasisMatrix object. Therefore, a DirectSparseSolver (DSS) object can be considered to be stateless if needed. The first point to make clear is that once the client has a reference to a BasisMatrix object in a smart pointer like C1_ptr, that BasisMatrix object will not be altered by any future operations on the DSS object that created it unless the client specifically requests an alteration. This behavior allows BasisMatrix objects to be completely decoupled from the DSS object that created them. This behavior is explained in more detail later.

It is through the MatrixNonsing interface of the BasisMatrix object that clients can solve for linear systems.

The usage of this interface will be illustrated by considering a few different scenarios:

  // Storage for the permutations and rank
  IVector row_perm, col_perm;
  size_type rank;
  
  // Analyze and factor A1
  typedef DirectSparseSolver::basis_matrix_ptr_t C_ptr_t;
  C_ptr_t  C1_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.analyze_and_factor( A1, &row_perm, &col_perm, &rank, C1_ptr.get() );
  // After the above method finishes direct_solver.get_fact_struc().get() points
  // to the same factorization structure as C1_ptr->get_fact_struc().get().

  // Solve for x1 = inv(C1)*b1
  VectorMutableDense x1(C1_ptr->cols());
  V_InvMtV( &x1, *C1_ptr, BLAS_Cpp::no_trans, b1 );

  // Factor another matrix with the same structure.
  // The current factorization structure in direct_solver.get_fact_struc() is used.
  C_ptr_t  C2_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.factor( A2, C2_ptr.get() );

  // Solve for x2 = inv(C2')*b2
  VectorMutableDense x2(C2_ptr->rows());
  V_InvMtV( &x2, *C2_ptr, BLAS_Cpp::trans, b2 );

2) Maintain the factorization of three unsymmetic matrices, two with the same structure and one with a different structure.

  // Storage for the permutations and rank
  IVector row_perm1, col_perm1;
  size_type rank1;
  IVector row_perm2, col_perm2;
  size_type rank2;

  // Analyze and factor A1
  typedef DirectSparseSolver::basis_matrix_ptr_t C_ptr_t;
  C_ptr_t  C1_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.analize_and_factor( A1, &row_perm1, &col_perm1, &rank1, C1_ptr.get() );

  // Solve for x1 = inv(C1)*b1
  VectorMutableDense x1(C1_ptr->cols());
  V_InvMtV( &x1, *C1_ptr, BLAS_Cpp::no_trans, b1 );

  // Factor another matrix A2 with a different structure from A1.
  C_ptr_t  C2_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.analize_and_factor( A2, &row_perm2, &col_perm2, &rank2, C2_ptr.get() );
  // In the above method, a new current factorization structure will be computed and
  // used for C2.  This is because the current factorization structure before the
  // method is called is being used by the basis matrix object *C1_ptr and must
  // be preserved.

  // Solve for x2 = inv(C2)*b2
  VectorMutableDense x2(C2_ptr->cols());
  V_InvMtV( &x2, *C2_ptr, BLAS_Cpp::no_trans, b2 );

  // Factor another matrix A3 with the same structure of as A1 but preserve
  // the factorization nonzeros of A1.
  C_ptr_t  C3_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.factor( A3, C3_ptr.get(), C1_ptr->get_fact_struc() );
  // Above, the current factorization structure direct_solver.get_fact_struc()
  // will be preserved before and after this method call.

  // Solve for x3 = inv(C3)*b3
  VectorMutableDense x3(C3_ptr->cols());
  V_InvMtV( &x3, *C3_ptr, BLAS_Cpp::no_trans, b3 );

3) Factor of two unsymmetic matrices with different structure and then recycle the storage of the first matrix to factor a third matrix with an all new structure.

  // Storage for the permutations and rank
  IVector row_perm1, col_perm1;
  size_type rank1;
  IVector row_perm2, col_perm2;
  size_type rank2;

  typedef DirectSparseSolver::basis_matrix_ptr_t  C_ptr_t;

  // Analyze and factor A1
  C_ptr_t  C1_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.analize_and_factor( A1, &row_perm1, &col_perm1, &rank1, C1_ptr.get() );

  // Solve for x1 = inv(C1)*b1
  VectorMutableDense x1(C1_ptr->cols());
  V_InvMtV( &x1, *C1_ptr, BLAS_Cpp::no_trans, b1 );

  // Factor another matrix A2 with a different structure from A1.
  C_ptr_t  C2_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.analize_and_factor( A2, &row_perm2, &col_perm2, &rank2, C2_ptr.get() );

  // Solve for x2 = inv(C2)*b2
  VectorMutableDense x2(C2_ptr->cols());
  V_InvMtV( &x2, *C2_ptr, BLAS_Cpp::no_trans, b2 );

  // Analyze and factor another matrix A3 with an all new structure and recycle storage of C1.
  C_ptr_t  C3_ptr = C1_ptr;
  C1_ptr = Teuchos::null;
  direct_solver.analyze_and_factor( A3, &row_perm1, &col_perm1, &rank1, C3_ptr.get() );

  // Solve for x3 = inv(C3)*b3
  VectorMutableDense x3(C1_ptr->cols());
  V_InvMtV( &x3, *C3_ptr, BLAS_Cpp::no_trans, b3 );

  // Factor another matrix A4 with the same structure as A2 but preserve
  // its factorization in the basis matrix C2.
  C_ptr_t  C4_ptr = direct_solver.basis_matrix_factory()->create();
  direct_solver.factor( A4, C4_ptr.get(), C2_ptr->get_fact_struc() );

  // Solve for x4 = inv(C4)*b4
  VectorMutableDense x4(C4_ptr->cols());
  V_InvMtV( &x4, *C4_ptr, BLAS_Cpp::no_trans, b4 );

  // Analyze and factor another matrix A5 with an all new structure
  // and recycle nonzero storage of C4.  Note that the behavoir of C2_ptr must
  // not change by this operation.  Therefore, a new factorization structure
  // must be allocated in this case.
  C_ptr_t  C5_ptr = C4_ptr;
  C4_ptr = Teuchos::null;
  direct_solver.analyze_and_factor( A5, &row_perm5, &col_perm5, &rank5, C4_ptr.get() );

If by some major mistake the client were to pass in basis_matrix or fact_struc objects of the incorrect type from non-compatible DSS objects then the InvalidObjectType exception would be thrown. If the client follows the protocal defined in this interface, this will never happen.

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

virtual void AbstractLinAlgPack::DirectSparseSolver::estimated_fillin_ratio ( value_type estimated_fillin_ratio ) [pure virtual]

Set the estimate of the fill-in ratio of the factorization.

This value is used to set the initial guess at the amount of storage needed for the factorization of the matrix plus some "elbow" room for the factorization algorithm.

The estimated amount of nonzeros for the factorization is:

num_factor_nonzeros = estimated_fillin_ratio * A.nz()

Implemented in AbstractLinAlgPack::DirectSparseSolverDense, and AbstractLinAlgPack::DirectSparseSolverMA28.

virtual void AbstractLinAlgPack::DirectSparseSolver::analyze_and_factor	(	const AbstractLinAlgPack::MatrixConvertToSparse &	A,
		DenseLinAlgPack::IVector *	row_perm,
		DenseLinAlgPack::IVector *	col_perm,
		size_type *	rank,
		BasisMatrix *	basis_matrix,
		std::ostream *	out = `NULL`
	)	`[pure virtual]`

Analyze and factor a matrix.

Parameters:

	A	[in] Matrix who's elements are extracted in order to form the factorization. Here `A.rows() <= A.cols()` must be `true`.
	row_perm	[out] On output holds an array (size `A.rows()`) of row permutations `P` used to select the basis matrix `C`. The ith row of `P'AQ` is row (*row_perm)(i) of `A`.
	col_perm	[out] On output holds an array (size `A.cols()`) of column permutations `Q` used to select the basis matrix `C`. The jth column of `P'AQ` is column (*col_perm)(j) of `A`. If the client is solving a symmetric system then `col_perm` may be set to `NULL` on input which forces the implementation to use symmetric permutations (see above discussion).
	rank	[out] On output holds the rank `r` (`r <= A.rows()`) of the basis matrix `C`.
	basis_matrix	[in/out] On input, this can be a previously initialized basis matrix object in which case its storage for the structure of the factorization and its nonzeros can be recycled if possible.
	out	[in/out] If `out!=NULL`, then some information about the operations performed internally may be printed to `*out`. The amount of this output should be very minimal and should not significantly increase with the size of the problem being solved.

Preconditions:

A.rows() <= A.cols() (throw std::invalid_argument)
basis_matrix != NULL (throw std::invalid_argument)

Postconditions:

*row_perm, *col_perm and *rank give the basis selection determined by the direct sparse solver implementation.
this->get_fact_struc().get() == basis_matrix->get_fact_struc().get() points to the factorizztion structure for this matrix.

This function extracts the nonzero elements from the matrix A (using the AbstractLinAlgPack::MatrixConvertToSparse interface) and then finds a nonsingular basis matrix and factors it (see intro).

Given the m x n matrix A (m <= n), this method finds the row and column permutations P and Q respectively and the rank r such that the basis matrix C = (P'*A*Q)(1:r,1:r) is nonsingular.

If a symmetric system is being solved and it is feared that A is not full rank and the client wants to use symmetric pivoting (P = Q) to find the basis C, then the client can pass in NULL for col_perm and force the solver to return a symmetric permutation. Note that the row and column permutations of row_perm(k) and col_perm(k) for k = 1...r does not neccesarly indicate the exact permutations used within the sparse solver. Instead, what is significant is the partitioning of rows and columns between k = 1..r and k = r+1..m for row permutations and k = r+1..n for column permutations. Therefore, if col_perm==NULL on input and an unsymmetric solver is used and the matrix is not full rank, then since in reality unsymmetric pivoting is being used internally, the implementation must throw an UnsymmetricRankDeficientException exception. This convention is needed to allow an unsymmetric solver to be used for a symmetric system but still meet the needs of users for symmetric systems. In summary, as long as the symmetric system is full rank, then an unsymmetic solver can always be used and the client need not care how the system is solved. On the other hand, if the client does not care if unsymmetric pivoting is used on a symmetric rank deficient system, then the client can set col_perm to point to a valid IVector object and deal with the unsymmetic permutations that are returned. Who knows what the use of this would be but there is no reason not to allow it.

This method allows for the careful sharing and recycling of the factorization structure and nonzeros of a factorized matrix while also maintaining a "current" factorization structure. To explain how this works let C1_ptr be set to the smart pointer to a basis matrix object created by this->basis_matrix_factory()->create().

If the client whats to recycle the factorization structure and nonzeros that are contained in *C1_ptr to analyze and factor A3, then the client would perform:

   direct_solver.analyze_and_factor(A3,&row_perm,&col_perm,&rank,C1_ptr.get());

If the factorization can not be performed for some reason then the exception FactorizationFailure will be thrown.

Scenarios / behavior:

basis_matrix->get_fact_struc().get()==NULL
- this->get_fact_struc().get()==NULL
  - New storage for the factorization structure will have to be allocated.
- this->get_fact_struc().get()!=NULL
  - this->get_fact_struc().count()==1
    - No other basis matrix object is using this factorization structure so the storage in this object can be recycled and used here.
  - this->get_fact_struc().count()>1
    - Some other basis matrix object or some other client is using this factorization structure so new storage for the factorization structure must be allocated.
- New storage for factorization nonzeros will have to be allocated.
basis_matrix->get_fact_struc().get()!=NULL
- this->get_fact_struc().get()==NULL
  - basis_matrix->get_fact_struc().count()==1
    - No other basis matrix ojbect is using this factorization structure so the storage in *basis_matrix->get_fact_struc() can be recycled and used here.
  - basis_matrix->get_fact_struc().count()>1
    - Some other basis matrix object or some other client is using this factorization structure so new storage for the factorization structure must be allocated.
- this->get_fact_struc().get()!=NULL
  - basis_matrix->get_fact_struc().get() == this->get_fact_struc().get()
    - basis_matrix->get_fact_struc().count()==2
      - These are the same factorization structure objects and no other basis matrix object is using this factorization structure so the storage in this object can be recycled and used here.
    - basis_matrix->get_fact_struc().count()>2
      - Some other basis matrix object or some other client is using this factorization structure so new storage for the factorization structure must be allocated.
  - basis_matrix->get_fact_struc().get() != this->get_fact_struc().get()
    - this->get_fact_struc().count()==1
      - basis_matrix->get_fact_struc().count()==1
        
        One of the factorization structure objects can be discarded and the other can be recycled here since no other basis system object or other client has a reference to these objects.
      - basis_matrix->get_fact_struc().count()>1
        
        The factorization structure associated with the basis matrix object basis_matrix is being shared by some other basis matrix object or other client but the current factorization structure in this is not. Therefore, we can recycle the storage in *this->get_fact_struc() here.
    - this->get_fact_struc().count()>1
      - basis_matrix->get_fact_struc().count()==1
        
        The factorization structure associated with the basis matrix object basis_matrix is not being shared by some other basis matrix object or other client but the current factorization structure in this is. Therefore, we can recycle the storage in *basis_matrix->get_fact_struc() here.
      - basis_matrix->get_fact_struc().count()>1
        
        The factorization structure objects associated with the basis matrix object basis_matrix and this are being shared by some other basis matrix objects or other clients. Therefore, we will have to allocate new factorization storage here.
- Current storage for the factorization nonzeros in *basis_matrix will be recycled here.

Implemented in AbstractLinAlgPack::DirectSparseSolverImp.

virtual void AbstractLinAlgPack::DirectSparseSolver::factor	(	const AbstractLinAlgPack::MatrixConvertToSparse &	A,
		BasisMatrix *	basis_matrix,
		const BasisMatrix::fact_struc_ptr_t &	fact_struc = `Teuchos::null`,
		std::ostream *	out = `NULL`
	)	`[pure virtual]`

Factor a matrix given its factorization structure.

Parameters:

	A	[in] Matrix to be factored given a previously determined factorization structure (see `fact_struc`).
	basis_matrix	[in/out] On input, this can be a previously initialized basis matrix object in which case its storage for the factorization nonzeros can be recycled if possible.
	fact_struc	[in] Smart pointer to the factorization structure to use. If `fact_struc.get()==NULL` then `this->get_fact_struc()` is used in its place.
	out	[in/out] If `out != NULL`, then some information about the operations performed internally may be printed to `*out`. The amount of this output should be very minimal and should not significantly increase with the size of the problem being solved.

Preconditions:

A has the same structure as used to form fact_struc (or this->get_fact_struc() if fact_struc.get()==NULL )
this->get_fact_struc().get() != NULL || fact_struc.get() != NULL (throw NoCurrentBasisException)
basis_matrix != NULL (throw std::invalid_argument)

Postconditions:

[fact_struc.get() != NULL] basis_matrix->get_fact_struc().get() == fact_struc.get() and fact_struc.count() is incremented by one on output.
[fact_struc.get() == NULL] basis_matrix->get_fact_struc().get() == this->get_fact_struc().get() and this->get_fact_struc().count() is incremented by one on output.
*this->get_fact_struc() is unchanged in any way (however, this->get_fact_struc().cout() will be incremented by one if fact_struct.get() == NULL).
If basis_matrix->get_fact_struct().get() != fact_struc.get() (or this->get_fact_struc().get() if fact_struc.get()==NULL) on input then if basis_matrix->get_fact_struct().count() == 1 on input then this factorization structure will be discarded on output.
If basis_matrix->get_fact_struc().get() == NULL on input, then new factorization nonzeros are allocated on output to hold the matrix factors. If basis_matrix->get_fact_struc().get() != NULL on input, then the storage for nonzeros in *basis_matrix is recycled.

This method allows clients to factor a matrix given it predetermined factorization structure. The client can use the factorization structure stored internally in this->fact_struc() (fact_struc.get()==NULL) or use the factorization structure from another precomputed basis matrix (fact_struc.get()!=NULL). If the passed in matrix A is obviously not compatible with the precomputed factorization structure being used then an IncompatibleMatrixStructureException exception will be thrown. If A1 was the matrix originally passed to this->analyze_and_factor(A1,...) and A2 is the matrix being passed to this->factor(A2,...) then if A2.rows()!=A1.rows() or A2.cols()!=A1.cols() or A2.nz()!=A1.nz() then the matrices are obviously not compatible and the exception will be thrown.

The argument basis_matrix can be used to recycle the nonzero values of the factorization of an existing basis matrix. For example, suppose that C1_ptr and C2_ptr point to two different basis matrices with different structure (i.e. C1_ptr->get_fact_struc().get() != C2_ptr->get_fact_struc().get()). Now suppose we would like to factor another matrix A3 that has the same structure of A2 by reusing the nonzero values referenced in C1_ptr. To do this the client would call:

   direct_solver.factor( A3, C1_ptr.get(), C2_ptr->get_fact_struc() );

Now consider the case where we would like to recycle the storage for the nonzeros in C1_ptr and also use the same structure as C1. To do this the client would call:

   direct_solver.factor(A3,C1_ptr.get(),C1_ptr->get_fact_struc());

Above, the client must explicitly pass C1_ptr->get_fact_struc() to use this structure or the current interally stored structure given by this->get_fact_struc() would be used instead.

If the factorization can not be performed for some reason then the exception FactorizationFailure will be thrown.

Implemented in AbstractLinAlgPack::DirectSparseSolverImp.

virtual const BasisMatrix::fact_struc_ptr_t& AbstractLinAlgPack::DirectSparseSolver::get_fact_struc ( ) const [pure virtual]

Get a smart pointer object to the current factorization structure.

This returns a smart reference counted pointer to the factorization structure computed from the last call to this->analyze_and_factor(). If no successful call to this->analyze_and_factor() has been made yet then this function will return return.get() == NULL.

Implemented in AbstractLinAlgPack::DirectSparseSolverImp.

virtual void AbstractLinAlgPack::DirectSparseSolver::set_uninitialized ( ) [pure virtual]

Release all allocated memory.

Postconditions:

this->get_fact_struc().get()==NULL.

Implemented in AbstractLinAlgPack::DirectSparseSolverImp.


Public Types
typedef Teuchos::RefCountPtr< const Teuchos::AbstractFactory< BasisMatrix > >	basis_matrix_factory_ptr_t

Public Member Functions
virtual	~DirectSparseSolver ()

virtual const basis_matrix_factory_ptr_t	basis_matrix_factory () const =0
	Return a factory object that can create the basis matrix.
virtual void	estimated_fillin_ratio (value_type estimated_fillin_ratio)=0
	Set the estimate of the fill-in ratio of the factorization.
virtual void	analyze_and_factor (const AbstractLinAlgPack::MatrixConvertToSparse &A, DenseLinAlgPack::IVector row_perm, DenseLinAlgPack::IVector col_perm, size_type rank, BasisMatrix basis_matrix, std::ostream *out=NULL)=0
	Analyze and factor a matrix.
virtual void	factor (const AbstractLinAlgPack::MatrixConvertToSparse &A, BasisMatrix basis_matrix, const BasisMatrix::fact_struc_ptr_t &fact_struc=Teuchos::null, std::ostream out=NULL)=0
	Factor a matrix given its factorization structure.
virtual const BasisMatrix::fact_struc_ptr_t &	get_fact_struc () const =0
	Get a smart pointer object to the current factorization structure.
virtual void	set_uninitialized ()=0
	Release all allocated memory.

AbstractLinAlgPack::DirectSparseSolver Class Reference

Public Types

Public Member Functions

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation