#include <AbstractLinAlgPack_DirectSparseSolver.hpp>
Inheritance diagram for AbstractLinAlgPack::DirectSparseSolver:
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 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. 
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 226 of file AbstractLinAlgPack_DirectSparseSolver.hpp.

Definition at line 267 of file AbstractLinAlgPack_DirectSparseSolver.hpp. 

Definition at line 297 of file AbstractLinAlgPack_DirectSparseSolver.hpp. 

Return a factory object that can create the basis matrix.
Implemented in AbstractLinAlgPack::DirectSparseSolverDense, and AbstractLinAlgPack::DirectSparseSolverMA28. 

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. 

Analyze and factor a matrix.
Postconditions:
This function extracts the nonzero elements from the matrix
Given the
If a symmetric system is being solved and it is feared that
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
If the client whats to recycle the factorization structure and nonzeros that are contained in
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 Scenarios / behavior:
Implemented in AbstractLinAlgPack::DirectSparseSolverImp. 

Factor a matrix given its factorization structure.
Postconditions:
This method allows clients to factor a matrix given it predetermined factorization structure. The client can use the factorization structure stored internally in
The argument
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
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
If the factorization can not be performed for some reason then the exception Implemented in AbstractLinAlgPack::DirectSparseSolverImp. 

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 Implemented in AbstractLinAlgPack::DirectSparseSolverImp. 

Release all allocated memory. Postconditions:
Implemented in AbstractLinAlgPack::DirectSparseSolverImp. 