AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day

Abstract interface to serial direct sparse linear solvers. More...
#include <AbstractLinAlgPack_DirectSparseSolver.hpp>
Classes  
class  BasisMatrix 
Abstract class for objects that represent the factorized matrix and can be used to solve for different righthandsides. More...  
class  FactorizationFailure 
More...  
class  FactorizationStructure 
Abstract class for objects that represent the factorization structure of a particular matrix. More...  
class  IncompatibleMatrixStructureException 
More...  
class  InvalidObjectType 
More...  
class  NoCurrentBasisException 
More...  
class  UnsymmetricRankDeficientException 
More...  
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 fillin 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.  
Public Types  
typedef Teuchos::RCP< const Teuchos::AbstractFactory < BasisMatrix > >  basis_matrix_factory_ptr_t 

Abstract interface to serial direct sparse linear solvers.
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 fillin 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:
1) Maintain one factorization structure of an unsymmetric matrix.
// 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 noncompatible DSS objects then the InvalidObjectType
exception would be thrown. If the client follows the protocal defined in this interface, this will never happen.
Definition at line 239 of file AbstractLinAlgPack_DirectSparseSolver.hpp.
typedef Teuchos::RCP< const Teuchos::AbstractFactory<BasisMatrix> > AbstractLinAlgPack::DirectSparseSolver::basis_matrix_factory_ptr_t 
Definition at line 280 of file AbstractLinAlgPack_DirectSparseSolver.hpp.
virtual AbstractLinAlgPack::DirectSparseSolver::~DirectSparseSolver  (  )  [inline, virtual] 
Definition at line 310 of file AbstractLinAlgPack_DirectSparseSolver.hpp.
virtual const basis_matrix_factory_ptr_t AbstractLinAlgPack::DirectSparseSolver::basis_matrix_factory  (  )  const [pure virtual] 
Return a factory object that can create the basis matrix.
Implemented in AbstractLinAlgPack::DirectSparseSolverDense, and AbstractLinAlgPack::DirectSparseSolverMA28.
virtual void AbstractLinAlgPack::DirectSparseSolver::estimated_fillin_ratio  (  value_type  estimated_fillin_ratio  )  [pure virtual] 
Set the estimate of the fillin 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.
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'*A*Q 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'*A*Q 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 this>get_fact_struc()
.get()!=NULL this>get_fact_struc()
.count()==1 this>get_fact_struc()
.count()>1 basis_matrix>get_fact_struc()
.get()!=NULL this>get_fact_struc()
.get()==NULL basis_matrix>get_fact_struc()
.count()==1 *basis_matrix>get_fact_struc()
can be recycled and used here.basis_matrix>get_fact_struc()
.count()>1 this>get_fact_struc()
.get()!=NULL basis_matrix>get_fact_struc().get() == this>get_fact_struc().get()
basis_matrix>get_fact_struc()
.count()==2 basis_matrix>get_fact_struc()
.count()>2 basis_matrix>get_fact_struc().get() != this>get_fact_struc().get()
this>get_fact_struc()
.count()==1 basis_matrix>get_fact_struc()
.count()==1 basis_matrix>get_fact_struc()
.count()>1 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 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 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.*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.
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
). 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. 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.
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.