AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects
Version of the Day
|
Mix-in Interface for updating a serial symmetric matrix by adding and deleting rows and columns. More...
#include <AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp>
Classes | |
struct | Inertia |
Struct for the inertia of the matrix. More... | |
class | MaxSizeExceededException |
Thrown if the maximum size is exceeded in augment_update(...). More... | |
struct | PivotTolerances |
Struct for pivot tolerances to be used when initializing, and augmenting and deleting rows and columns. More... | |
class | SingularUpdateException |
Thrown if the matrix is singular and should not have been. More... | |
class | WarnNearSingularUpdateException |
Thrown if the matrix is near singular as a warning. More... | |
class | WrongInertiaUpdateException |
Thrown if matrix has the wrong inertia from what was expected. More... | |
Public types | |
enum | EEigenValType |
Public members to be overridden | |
virtual | ~MatrixSymAddDelUpdateable () |
virtual void | initialize (value_type alpha, size_type max_size)=0 |
Initialize to a 1x1 matrix. More... | |
virtual void | initialize (const DMatrixSliceSym &A, size_type max_size, bool force_factorization, Inertia inertia, PivotTolerances pivot_tols=PivotTolerances())=0 |
Initialize given a symmetric matrix. More... | |
virtual size_type | max_size () const =0 |
Return the maximum size the matrix is allowed to become. More... | |
virtual Inertia | inertia () const =0 |
Return the inertia of the matrix (if it is known). If any of the members of the inertia is not known then they may be set to Inertia::UNKNOWN . If the matrix is nonsingular then return.zero_eigens == 0 will be true. More... | |
virtual void | set_uninitialized ()=0 |
Set the matrix to uninitialized. More... | |
virtual void | augment_update (const DVectorSlice *t, value_type alpha, bool force_refactorization=true, EEigenValType add_eigen_val=EIGEN_VAL_UNKNOWN, PivotTolerances pivot_tols=PivotTolerances())=0 |
Update by adding a symmetric row and column. More... | |
virtual void | delete_update (size_type jd, bool force_refactorization=true, EEigenValType drop_eigen_val=EIGEN_VAL_UNKNOWN, PivotTolerances pivot_tols=PivotTolerances())=0 |
Update by deleteing a symmetric row and column. More... | |
Additional Inherited Members | |
Public Member Functions inherited from AbstractLinAlgPack::MatrixBase | |
virtual | ~MatrixBase () |
Virtual destructor. More... | |
virtual const VectorSpace & | space_cols () const =0 |
Vector space for vectors that are compatible with the columns of the matrix. More... | |
virtual const VectorSpace & | space_rows () const =0 |
Vector space for vectors that are compatible with the rows of the matrix. More... | |
virtual size_type | rows () const |
Return the number of rows in the matrix. More... | |
virtual size_type | cols () const |
Return the number of columns in the matrix. More... | |
virtual size_type | nz () const |
Return the number of nonzero elements in the matrix. More... | |
Mix-in Interface for updating a serial symmetric matrix by adding and deleting rows and columns.
This interface is designed to allow objects of this type to be used in several different situations. Generally, the matrix being updated would be nonsingular but it does not have to be. An important property of symmetric matrix is its inertia. The inertia of a matrix is the number of negative, zero, and positive eigenvalues respectively. While the eigenvalues themselves can be very difficult to compute, in many cases the inertia is very easy to determine. For instance, if a A = L*D*L'
factorization is being used, the inertia of the diagonal matrix D
is the same as for A
. Likewise, if a cholesky factorization A = (+-)C*C'
is used then it is easy to prove if the matrix is positive definite (all positive eigen values) or negative definite (all negative eigen values). With other factorizations (such as QR for instance) it is more difficult to determine the inertia and therefore it may not be available.
In inexact floating point arithmetic, it can be difficult to distingish if a matrix is singular or nonsingular or if a matrix really has the wrong inertia or is just singular. In order to do this, tolerances have to be identified. In many contexts it is important for the client to be able to specify these tolerances. In order to establish a frame of reference that is independent of the actual implementation of the factorizations for the subclasses of this interface, we will use the LU factorization:
A = L*U
where L
is lower unit triangular and U
is upper nonunit triangular. We will define the quantity:
gamma = min{|U(i,i)|,i=1..n}/max{|U(i,i)|,i=1..n}
as a measure of singularity. Of course subclasses may not actually use a LU factorization but by establishing this simple frame of reference the tolerances can be properly interpreted by the subclasses.
The classification of an factorized matrix will be as follows:
if (correct inertia) and (gamma > warning_tol) then The matrix is nonsingular and has the correct inertia, the initialization or update will succeed and all is good :-) elseif (correct inertia) and (singular_tol < gamma <= warning_tol) then The matrix will be considered nonsingular and the initialization or the update will succeed but a WarnNearSingularUpdateException will be thrown containing gamma and a warning message. elseif (correct inertia) and (0.0 < gamma <= singular_tol) then The matrix is considered singular, the initialization or update will not succeeed and a SingularUpdateException will be thrown containing gamma and an error message. elseif (gamma == 0.0) then The matrix is exactly singular, the initialization or update will not succeed and a SingularUpdateException will be thrown containing gamma and an error message. elseif (incorrect inertia) and (0.0 < gamma < wrong_inertia_tol) then The matrix will be considered singular, the initialization or update will not succeed and a SingularUpdateException will be thrown containing gamma and an error message. elseif (incorrect inertia) and (gamma >= wrong_inertia_tol) then The matrix is considered to be nonsingular but to have the wrong inertia, the initialization or update will not succeed and a WrongInertiaException will be thrown containing gamma and an error message. endif
The tolerances warning_tol
, singular_tol
and wrong_inertia_tol
are passed in as part of the struct PivotTolerances
to each of the methods that need them. The default initialization for these tolerances is PivotTolerances::UNKNOWN
which means that the matrix object can do what ever it wants to do. It may use its own tolerances or none at all. In other words, the behavior is completely up to the subclass. The idea of warning_tol
and throwing an exception except
of type WarnNearSingularUpdateException
is to allow the updates to succeed but return a warning message and the value of gamma
as information to the user.
Definition at line 126 of file AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp.
Definition at line 135 of file AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp.
|
inlinevirtual |
Definition at line 216 of file AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp.
|
pure virtual |
Initialize given a symmetric matrix.
The behavior of this function will vary based on the subclass that implements it. Some subclasses may require that A
be nonsingular and therefore inertia.zero_eigens
should be zero.
A | [in] Symetric matrix that this is initialized with. |
max_size | [in] The maximum size rows() and cols() can become. |
force_factorization | [in] If true, the factorization of the matrix will be forced and any possible exceptions will be thrown. If false then the factorization may not be forced, in which case the client may not know immediatly that the matrix is singular or has the wrong inertia. |
inertia | [in] The estimated inertia of the matrix. If the user knows any of the members of inertia then they should be set. Some subclasses may rely on this estimate of the inertia to determine what should be done. |
pivot_tols | [in] Tolerances to use to determine singularity, nonsingularity etc. See the intro. Default is no tolerances. |
|
pure virtual |
Return the maximum size the matrix is allowed to become.
|
pure virtual |
Return the inertia of the matrix (if it is known). If any of the members of the inertia is not known then they may be set to Inertia::UNKNOWN
. If the matrix is nonsingular then return.zero_eigens == 0
will be true.
|
pure virtual |
Set the matrix to uninitialized.
|
pure virtual |
Update by adding a symmetric row and column.
The update performed is:
[ A t ] [ t' alpha ] ==> A_new
Preconditions:
{itemize} [t != NULL] t->size() == this->rows()
(throw std::length_error
) this->rows() < this->max_size()
(throw MaxSizeExceededException
) {itemize}
Postcondiditons:
The update gives a legal update depending on the context of the subclass (nonsigular, positive definite etc.). If the subclass requires the matrix to be nonsingular but inertia.zero_eigens == 0
or the matrix is determined to be singular then the exception SingularUpdateException
will be thrown. If the matrix is found to not have the propper inertia then the exception WrongInertiaUpdateException
will be thrown. This subclass may not be able to determine the inertia in which case this exception will never be thrown. If no exceptions are thrown then this->rows()
and this->cols()
will increase by one and this->inertia()
will return the new inertia if it is known.
t | [in] DVectorSlice (size == rows() ) where t may be NULL in which case t is considered zero. |
alpha | [in] Scalar added. |
force_refactorization | [in] If true, then the factorization of the matrix will be performed before the function returns. If something goes wrong then an exeception will be thrown here. |
add_eigen_val | [in] Gives the estimate of the new eigen value added to the matrix. If the matrix does not agree with this then an exception will be thrown. |
pivot_tols | [in] Tolerances to use to determine singularity, nonsingularity etc. See the intro. Default is no tolerances. |
|
pure virtual |
Update by deleteing a symmetric row and column.
jd [ A11 a12 A13 ] A = [ a12' a22 a23' ] jd ==> A_new = [ A11 A13 ] [ A13' a23 A33 ] [ A13' A33 ]
Preconditions:
{itemize} 1 <= jd && jd <= this->rows()
(throw std::out_of_range
) {itemize}
Postcondiditons:
The update give a legal update depending on the context of the subclass (nonsigular, positive definite etc.). Also rows()
and cols()
will decrease by one so this_after->rows()
== this_before->rows()
- 1.
jd | [in] The jth row and column to be removed from the matrix. |
force_refactorization | [in] If true, then the factorization of the matrix will be performed before the function returns. If something goes wrong then an exeception will be thrown here. |
drop_eigen_val | [in] Gives the estimate of the eigen value dropped from the matrix. If the matrix does not agree with this then an exception will be thrown. |
pivot_tols | [in] Tolerances to use to determine singularity, nonsingularity etc. See the intro. Default is no tolerances. |