Reference documentation for deal.II version Git 035d1ad 2018-08-17 19:42:27 -0400
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
LAPACKFullMatrix< number > Class Template Reference

#include <deal.II/lac/lapack_full_matrix.h>

Inheritance diagram for LAPACKFullMatrix< number >:
[legend]

Public Types

using size_type = std::make_unsigned< types::blas_int >::type
 
- Public Types inherited from TransposeTable< number >
using size_type = typename TableBase< 2, number >::size_type
 
using value_type = typename AlignedVector< number >::value_type
 
using reference = typename AlignedVector< number >::reference
 
using const_reference = typename AlignedVector< number >::const_reference
 
using const_iterator = TransposeTableIterators::Iterator< number, true >
 
using iterator = TransposeTableIterators::Iterator< number, false >
 
- Public Types inherited from TableBase< 2, number >
using size_type = typename AlignedVector< number >::size_type
 

Public Member Functions

 LAPACKFullMatrix (const size_type size=0)
 
 LAPACKFullMatrix (const size_type rows, const size_type cols)
 
 LAPACKFullMatrix (const LAPACKFullMatrix &)
 
LAPACKFullMatrix< number > & operator= (const LAPACKFullMatrix< number > &)
 
template<typename number2 >
LAPACKFullMatrix< number > & operator= (const FullMatrix< number2 > &)
 
template<typename number2 >
LAPACKFullMatrix< number > & operator= (const SparseMatrix< number2 > &)
 
LAPACKFullMatrix< number > & operator= (const number d)
 
LAPACKFullMatrix< number > & operator*= (const number factor)
 
LAPACKFullMatrix< number > & operator/= (const number factor)
 
void set (const size_type i, const size_type j, const number value)
 
void add (const number a, const LAPACKFullMatrix< number > &B)
 
void rank1_update (const number a, const Vector< number > &v)
 
template<typename MatrixType >
void copy_from (const MatrixType &)
 
void reinit (const size_type size)
 
void grow_or_shrink (const size_type size)
 
void remove_row_and_column (const size_type row, const size_type col)
 
void reinit (const size_type rows, const size_type cols)
 
void set_property (const LAPACKSupport::Property property)
 
size_type m () const
 
size_type n () const
 
template<typename MatrixType >
void fill (const MatrixType &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0, const number factor=1., const bool transpose=false)
 
template<typename number2 >
void vmult (Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
 
void vmult (Vector< number > &w, const Vector< number > &v, const bool adding=false) const
 
template<typename number2 >
void vmult_add (Vector< number2 > &w, const Vector< number2 > &v) const
 
void vmult_add (Vector< number > &w, const Vector< number > &v) const
 
template<typename number2 >
void Tvmult (Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
 
void Tvmult (Vector< number > &w, const Vector< number > &v, const bool adding=false) const
 
template<typename number2 >
void Tvmult_add (Vector< number2 > &w, const Vector< number2 > &v) const
 
void Tvmult_add (Vector< number > &w, const Vector< number > &v) const
 
void mmult (LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void mmult (FullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void Tmmult (LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void Tmmult (FullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void Tmmult (LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const Vector< number > &V, const bool adding=false) const
 
void mTmult (LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void mTmult (FullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void TmTmult (LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void TmTmult (FullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
 
void scale_rows (const Vector< number > &V)
 
void compute_lu_factorization ()
 
void compute_cholesky_factorization ()
 
number reciprocal_condition_number (const number l1_norm) const
 
number reciprocal_condition_number () const
 
number determinant () const
 
number l1_norm () const
 
number linfty_norm () const
 
number frobenius_norm () const
 
number trace () const
 
void invert ()
 
void solve (Vector< number > &v, const bool transposed=false) const
 
void solve (LAPACKFullMatrix< number > &B, const bool transposed=false) const
 
void apply_lu_factorization (Vector< number > &v, const bool transposed) const
 
void apply_lu_factorization (LAPACKFullMatrix< number > &B, const bool transposed) const
 
void compute_eigenvalues (const bool right_eigenvectors=false, const bool left_eigenvectors=false)
 
void compute_eigenvalues_symmetric (const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
 
void compute_generalized_eigenvalues_symmetric (LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number >> &eigenvectors, const types::blas_int itype=1)
 
void compute_generalized_eigenvalues_symmetric (LAPACKFullMatrix< number > &B, std::vector< Vector< number >> &eigenvectors, const types::blas_int itype=1)
 
void compute_svd ()
 
void compute_inverse_svd (const double threshold=0.)
 
void compute_inverse_svd_with_kernel (const unsigned int kernel_size)
 
std::complex< number > eigenvalue (const size_type i) const
 
number singular_value (const size_type i) const
 
void print_formatted (std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1., const double threshold=0.) const
 
- Public Member Functions inherited from TransposeTable< number >
 TransposeTable ()=default
 
 TransposeTable (const size_type size1, const size_type size2)
 
void reinit (const size_type size1, const size_type size2, const bool omit_default_initialization=false)
 
const_reference operator() (const size_type i, const size_type j) const
 
reference operator() (const size_type i, const size_type j)
 
size_type n_rows () const
 
size_type n_cols () const
 
iterator begin ()
 
const_iterator begin () const
 
iterator end ()
 
const_iterator end () const
 
- Public Member Functions inherited from TableBase< 2, number >
 TableBase ()=default
 
 TableBase (const TableIndices< N > &sizes)
 
 TableBase (const TableIndices< N > &sizes, InputIterator entries, const bool C_style_indexing=true)
 
 TableBase (const TableBase< N, number > &src)
 
 TableBase (const TableBase< N, T2 > &src)
 
 TableBase (TableBase< N, number > &&src) noexcept
 
 ~TableBase () override=default
 
TableBase< N, number > & operator= (const TableBase< N, number > &src)
 
TableBase< N, number > & operator= (const TableBase< N, T2 > &src)
 
TableBase< N, number > & operator= (TableBase< N, number > &&src) noexcept
 
bool operator== (const TableBase< N, number > &T2) const
 
void reset_values ()
 
void reinit (const TableIndices< N > &new_size, const bool omit_default_initialization=false)
 
size_type size (const unsigned int i) const
 
const TableIndices< N > & size () const
 
size_type n_elements () const
 
bool empty () const
 
void fill (InputIterator entries, const bool C_style_indexing=true)
 
void fill (const number &value)
 
AlignedVector< number >::reference operator() (const TableIndices< N > &indices)
 
AlignedVector< number >::const_reference operator() (const TableIndices< N > &indices) const
 
void swap (TableBase< N, number > &v)
 
std::size_t memory_consumption () const
 
void serialize (Archive &ar, const unsigned int version)
 
- Public Member Functions inherited from Subscriptor
 Subscriptor ()
 
 Subscriptor (const Subscriptor &)
 
 Subscriptor (Subscriptor &&) noexcept
 
virtual ~Subscriptor ()
 
Subscriptoroperator= (const Subscriptor &)
 
Subscriptoroperator= (Subscriptor &&) noexcept
 
void subscribe (const char *identifier=nullptr) const
 
void unsubscribe (const char *identifier=nullptr) const
 
unsigned int n_subscriptions () const
 
template<typename StreamType >
void list_subscribers (StreamType &stream) const
 
void list_subscribers () const
 
template<class Archive >
void serialize (Archive &ar, const unsigned int version)
 

Private Member Functions

number norm (const char type) const
 

Private Attributes

LAPACKSupport::State state
 
LAPACKSupport::Property property
 
std::vector< number > work
 
std::vector< types::blas_intiwork
 
std::vector< types::blas_intipiv
 
std::vector< number > inv_work
 
std::vector< typename numbers::NumberTraits< number >::real_type > wr
 
std::vector< number > wi
 
std::vector< number > vl
 
std::vector< number > vr
 
std::unique_ptr< LAPACKFullMatrix< number > > svd_u
 
std::unique_ptr< LAPACKFullMatrix< number > > svd_vt
 
Threads::Mutex mutex
 

Additional Inherited Members

- Static Public Member Functions inherited from Subscriptor
static::ExceptionBase & ExcInUse (int arg1, std::string arg2, std::string arg3)
 
static::ExceptionBase & ExcNoSubscriber (std::string arg1, std::string arg2)
 
- Protected Member Functions inherited from TransposeTable< number >
reference el (const size_type i, const size_type j)
 
const_reference el (const size_type i, const size_type j) const
 
- Protected Member Functions inherited from TableBase< 2, number >
size_type position (const TableIndices< N > &indices) const
 
AlignedVector< number >::reference el (const TableIndices< N > &indices)
 
AlignedVector< number >::const_reference el (const TableIndices< N > &indices) const
 
- Protected Attributes inherited from TableBase< 2, number >
AlignedVector< number > values
 
TableIndices< N > table_size
 

Detailed Description

template<typename number>
class LAPACKFullMatrix< number >

A variant of FullMatrix using LAPACK functions wherever possible. In order to do this, the matrix is stored in transposed order. The element access functions hide this fact by reverting the transposition.

Note
In order to perform LAPACK functions, the class contains a lot of auxiliary data in the private section. The names of these data vectors are usually the names chosen for the arguments in the LAPACK documentation.
Author
Guido Kanschat, 2005, Denis Davydov, 2017, 2018

Definition at line 42 of file full_matrix.h.

Member Typedef Documentation

template<typename number>
using LAPACKFullMatrix< number >::size_type = std::make_unsigned<types::blas_int>::type

Declare type for container size.

Definition at line 65 of file lapack_full_matrix.h.

Constructor & Destructor Documentation

template<typename number >
LAPACKFullMatrix< number >::LAPACKFullMatrix ( const size_type  size = 0)
explicit

Constructor. Initialize the matrix as a square matrix with dimension size.

In order to avoid the implicit conversion of integers and other types to a matrix, this constructor is declared explicit.

By default, no memory is allocated.

Definition at line 243 of file lapack_full_matrix.cc.

template<typename number >
LAPACKFullMatrix< number >::LAPACKFullMatrix ( const size_type  rows,
const size_type  cols 
)

Constructor. Initialize the matrix as a rectangular matrix \(\rm{rows} \times \rm{cols}\).

Definition at line 251 of file lapack_full_matrix.cc.

template<typename number >
LAPACKFullMatrix< number >::LAPACKFullMatrix ( const LAPACKFullMatrix< number > &  M)

Copy constructor. This constructor does a deep copy of the matrix. Therefore, it poses a possible efficiency problem, if for example, function arguments are passed by value rather than by reference. Unfortunately, we can't mark this copy constructor explicit, since that prevents the use of this class in containers, such as std::vector. The responsibility to check performance of programs must therefore remain with the user of this class.

Definition at line 259 of file lapack_full_matrix.cc.

Member Function Documentation

template<typename number>
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator= ( const LAPACKFullMatrix< number > &  M)

Assignment operator.

Definition at line 268 of file lapack_full_matrix.cc.

template<typename number >
template<typename number2 >
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator= ( const FullMatrix< number2 > &  M)

Assignment operator from a regular FullMatrix.

Note
Since LAPACK expects matrices in transposed order, this transposition is included here.

Definition at line 341 of file lapack_full_matrix.cc.

template<typename number >
template<typename number2 >
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator= ( const SparseMatrix< number2 > &  M)

Assignment operator from a regular SparseMatrix.

Note
Since LAPACK expects matrices in transposed order, this transposition is included here.

Definition at line 358 of file lapack_full_matrix.cc.

template<typename number>
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator= ( const number  d)

This operator assigns a scalar to a matrix. To avoid confusion with constructors, zero (when cast to the number type) is the only value allowed for d.

Definition at line 374 of file lapack_full_matrix.cc.

template<typename number>
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator*= ( const number  factor)

This operator multiplies all entries by a fixed factor.

Definition at line 389 of file lapack_full_matrix.cc.

template<typename number>
LAPACKFullMatrix< number > & LAPACKFullMatrix< number >::operator/= ( const number  factor)

This operator divides all entries by a fixed factor.

Definition at line 417 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::set ( const size_type  i,
const size_type  j,
const number  value 
)
inline

Set a particular entry of the matrix to a value. Thus, calling A.set(1,2,3.141); is entirely equivalent to the operation A(1,2) = 3.141;. This function exists for compatibility with the various sparse matrix objects.

Parameters
iThe row index of the element to be set.
jThe column index of the element to be set.
valueThe value to be written into the element.

Definition at line 983 of file lapack_full_matrix.h.

template<typename number>
void LAPACKFullMatrix< number >::add ( const number  a,
const LAPACKFullMatrix< number > &  B 
)

Simple addition of a scaled matrix, i.e. \(\mathbf A \mathrel{+}= a \, \mathbf B\).

Definition at line 448 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::rank1_update ( const number  a,
const Vector< number > &  v 
)

Perform a rank-1 update of a symmetric matrix \( \mathbf A \leftarrow \mathbf A + a \, \mathbf v \mathbf v^T \).

This function also works for Cholesky factorization. In that case, updating ( \(a>0\)) is performed via Givens rotations, whereas downdating ( \(a<0\)) via hyperbolic rotations. Note that the latter case might lead to a negative definite matrix in which case the error will be thrown (because Cholesky factorizations are only valid for symmetric and positive definite matrices).

Definition at line 575 of file lapack_full_matrix.cc.

template<typename number >
template<typename MatrixType >
void LAPACKFullMatrix< number >::copy_from ( const MatrixType &  M)
inline

Assignment from different matrix classes, performing the usual conversion to the transposed format expected by LAPACK. This assignment operator uses iterators of the typename MatrixType. Therefore, sparse matrices are possible sources.

Definition at line 1008 of file lapack_full_matrix.h.

template<typename number >
void LAPACKFullMatrix< number >::reinit ( const size_type  size)

Regenerate the current matrix by one that has the same properties as if it were created by the constructor of this class with the same argument list as this present function.

Definition at line 280 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::grow_or_shrink ( const size_type  size)

Same as above but will preserve the values of matrix upon resizing. The original values of the matrix are kept on increasing the size

\[ \mathbf A \rightarrow \left( \begin{array}{cc} \mathbf A & \mathbf 0 \\ \mathbf 0 & \mathbf 0 \end{array} \right) \]

Whereas if the new size is smaller, the matrix will contain the upper left block of the original one

\[ \left( \begin{array}{cc} \mathbf A_{11} & \mathbf A_{12} \\ \mathbf A_{21} & \mathbf A_{22} \end{array} \right) \rightarrow \mathbf A_{11} \]

Definition at line 290 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::remove_row_and_column ( const size_type  row,
const size_type  col 
)

Remove row row and column col from the matrix.

\[ \left( \begin{array}{ccc} \mathbf A_{11} & \mathbf a_{12} & \mathbf A_{13} \\ \mathbf a_{21}^T & a_{22} & \mathbf a_{23}^T \\ \mathbf A_{31} & \mathbf a_{32} & \mathbf A_{33} \end{array} \right) \rightarrow \left( \begin{array}{cc} \mathbf A_{11} & \mathbf A_{13} \\ \mathbf A_{31} & \mathbf A_{33} \end{array} \right) \]

Definition at line 304 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::reinit ( const size_type  rows,
const size_type  cols 
)

Regenerate the current matrix by one that has the same properties as if it were created by the constructor of this class with the same argument list as this present function.

Definition at line 331 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::set_property ( const LAPACKSupport::Property  property)

Assign property to this matrix.

Definition at line 1331 of file lapack_full_matrix.cc.

template<typename number >
LAPACKFullMatrix< number >::size_type LAPACKFullMatrix< number >::m ( ) const
inline

Return the dimension of the codomain (or range) space.

Note
The matrix is of dimension \(m \times n\).

Definition at line 993 of file lapack_full_matrix.h.

template<typename number >
LAPACKFullMatrix< number >::size_type LAPACKFullMatrix< number >::n ( ) const
inline

Return the dimension of the domain space.

Note
The matrix is of dimension \(m \times n\).

Definition at line 1000 of file lapack_full_matrix.h.

template<typename number>
template<typename MatrixType >
void LAPACKFullMatrix< number >::fill ( const MatrixType &  src,
const size_type  dst_offset_i = 0,
const size_type  dst_offset_j = 0,
const size_type  src_offset_i = 0,
const size_type  src_offset_j = 0,
const number  factor = 1.,
const bool  transpose = false 
)
inline

Fill rectangular block.

A rectangular block of the matrix src is copied into this. The upper left corner of the block being copied is (src_offset_i,src_offset_j). The upper left corner of the copied block is (dst_offset_i,dst_offset_j). The size of the rectangular block being copied is the maximum size possible, determined either by the size of this or src.

The final two arguments allow to enter a multiple of the source or its transpose.

Definition at line 1032 of file lapack_full_matrix.h.

template<typename number >
template<typename number2 >
void LAPACKFullMatrix< number >::vmult ( Vector< number2 > &  w,
const Vector< number2 > &  v,
const bool  adding = false 
) const

Matrix-vector-multiplication.

Depending on previous transformations recorded in state, the result of this function is one of

The optional parameter adding determines, whether the result is stored in the vector \(\mathbf w = \mathbf A \cdot \mathbf v\) or added to it \(\mathbf w \mathrel{+}= \mathbf A \cdot \mathbf v\).

Note
Source and destination must not be the same vector.
The template with number2 only exists for compile-time compatibility with FullMatrix. Only the case number2 = number is implemented due to limitations in the underlying LAPACK interface. All other variants throw an error upon invocation.

Definition at line 1066 of file lapack_full_matrix.h.

template<typename number>
void LAPACKFullMatrix< number >::vmult ( Vector< number > &  w,
const Vector< number > &  v,
const bool  adding = false 
) const

Specialization of above function for compatible Vector::value_type.

Definition at line 615 of file lapack_full_matrix.cc.

template<typename number >
template<typename number2 >
void LAPACKFullMatrix< number >::vmult_add ( Vector< number2 > &  w,
const Vector< number2 > &  v 
) const

Adding Matrix-vector-multiplication \(\mathbf w \mathrel{+}= \mathbf A \cdot \mathbf v\).

See the documentation of vmult() for details on the implementation.

Definition at line 1079 of file lapack_full_matrix.h.

template<typename number>
void LAPACKFullMatrix< number >::vmult_add ( Vector< number > &  w,
const Vector< number > &  v 
) const

Specialization of above function for compatible Vector::value_type.

Definition at line 886 of file lapack_full_matrix.cc.

template<typename number >
template<typename number2 >
void LAPACKFullMatrix< number >::Tvmult ( Vector< number2 > &  w,
const Vector< number2 > &  v,
const bool  adding = false 
) const

Transpose matrix-vector-multiplication.

The optional parameter adding determines, whether the result is stored in the vector \(\mathbf w = \mathbf A^T \cdot \mathbf v\) or added to it \(\mathbf w \mathrel{+}= \mathbf A^T \cdot \mathbf v\).

See the documentation of vmult() for details on the implementation.

Definition at line 1091 of file lapack_full_matrix.h.

template<typename number>
void LAPACKFullMatrix< number >::Tvmult ( Vector< number > &  w,
const Vector< number > &  v,
const bool  adding = false 
) const

Specialization of above function for compatible Vector::value_type.

Definition at line 749 of file lapack_full_matrix.cc.

template<typename number >
template<typename number2 >
void LAPACKFullMatrix< number >::Tvmult_add ( Vector< number2 > &  w,
const Vector< number2 > &  v 
) const

Adding transpose matrix-vector-multiplication \(\mathbf w \mathrel{+}= \mathbf A^T \cdot \mathbf v\).

See the documentation of vmult() for details on the implementation.

Definition at line 1104 of file lapack_full_matrix.h.

template<typename number>
void LAPACKFullMatrix< number >::Tvmult_add ( Vector< number > &  w,
const Vector< number > &  v 
) const

Specialization of above function for compatible Vector::value_type.

Definition at line 895 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::mmult ( LAPACKFullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Matrix-matrix-multiplication.

The optional parameter adding determines, whether the result is stored in the matrix \(\mathbf C = \mathbf A \cdot \mathbf B\) or added to it \(\mathbf C \mathrel{+}= \mathbf A \cdot \mathbf B\).

Note
It is assumed that A and B have compatible sizes and that C already has the right size.

This function uses the BLAS function Xgemm.

Definition at line 904 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::mmult ( FullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Same as before, but stores the result in a FullMatrix, not in a LAPACKFullMatrix.

Definition at line 938 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::Tmmult ( LAPACKFullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Matrix-matrix-multiplication using transpose of this.

The optional parameter adding determines, whether the result is stored in the matrix \(\mathbf C = \mathbf A^T \cdot \mathbf B\) or added to it \(\mathbf C \mathrel{+}= \mathbf A^T \cdot \mathbf B\).

Note
It is assumed that A and B have compatible sizes and that C already has the right size.
This function uses the BLAS function Xgemm.

Definition at line 1053 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::Tmmult ( FullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Same as before, but stores the result in a FullMatrix, not in a LAPACKFullMatrix.

Definition at line 1110 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::Tmmult ( LAPACKFullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const Vector< number > &  V,
const bool  adding = false 
) const

Matrix-matrix-multiplication using transpose of this and a diagonal vector V.

If the adding=false then the result is stored in the matrix \(\mathbf C = \mathbf A^T \cdot \rm{diag}(\mathbf V) \cdot \mathbf B\) otherwise it is added \(\mathbf C \mathrel{+}= \mathbf A^T \cdot \rm{diag}(\mathbf V) \cdot \mathbf B\).

Note
It is assumed that A, B and V have compatible sizes and that C already has the right size.
This function is not provided by LAPACK. The function first forms \(\rm{diag}(\mathbf V) \cdot \mathbf B\) product and then uses Xgemm function.

Definition at line 974 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::mTmult ( LAPACKFullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Matrix-matrix-multiplication using transpose of B.

The optional parameter adding determines, whether the result is stored in the matrix \(\mathbf C = \mathbf A \cdot \mathbf B^T\) or added to it \(\mathbf C \mathrel{+}= \mathbf A \cdot \mathbf B^T\).

Note
It is assumed that A and B have compatible sizes and that C already has the right size.
This function uses the BLAS function Xgemm.

Definition at line 1145 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::mTmult ( FullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Same as before, but stores the result in a FullMatrix, not in a LAPACKFullMatrix.

Definition at line 1203 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::TmTmult ( LAPACKFullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Matrix-matrix-multiplication using transpose of this and B.

The optional parameter adding determines, whether the result is stored in the matrix \(\mathbf C = \mathbf A^T \cdot \mathbf B^T\) or added to it \(\mathbf C \mathrel{+}= \mathbf A^T \cdot \mathbf B^T\).

Note
It is assumed that A and B have compatible sizes and that C already has the right size.
This function uses the BLAS function Xgemm.

Definition at line 1238 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::TmTmult ( FullMatrix< number > &  C,
const LAPACKFullMatrix< number > &  B,
const bool  adding = false 
) const

Same as before, but stores the result in a FullMatrix, not in a LAPACKFullMatrix.

Definition at line 1272 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::scale_rows ( const Vector< number > &  V)

Scale rows of this matrix by V . This is equivalent to premultiplication with a diagonal matrix \(\mathbf A\leftarrow {\rm diag}(\mathbf V)\mathbf A\).

Definition at line 1036 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_lu_factorization ( )

Compute the LU factorization of the matrix using LAPACK function Xgetrf.

Definition at line 1307 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_cholesky_factorization ( )

Compute the Cholesky factorization of the matrix using LAPACK function Xpotrf.

Note
The factorization is stored in the lower-triangular part of the matrix.

Definition at line 1421 of file lapack_full_matrix.cc.

template<typename number>
number LAPACKFullMatrix< number >::reciprocal_condition_number ( const number  l1_norm) const

Estimate the reciprocal of the condition number \(1/k(\mathbf A)\) in \(L_1\) norm ( \(1/(||\mathbf A||_1 \, ||\mathbf A^{-1}||_1)\)) of a symmetric positive definite matrix using Cholesky factorization. This function can only be called if the matrix is already factorized.

Note
The condition number \(k(\mathbf A)\) can be used to estimate the numerical error related to the matrix inversion or the solution of the system of linear algebraic equations as error = std::numeric_limits<Number>::epsilon * k. Alternatively one can get the number of accurate digits std::floor(std::log10(k)).
The function computes reciprocal of the condition number to avoid possible overflow if the matrix is nearly singular.
Parameters
[in]l1_normIs the \(L_1\) norm of the matrix before calling Cholesky factorization. It can be obtained by calling l1_norm().

Definition at line 1448 of file lapack_full_matrix.cc.

template<typename number>
number LAPACKFullMatrix< number >::reciprocal_condition_number ( ) const

Estimate the reciprocal of the condition number \(1/k(\mathbf A)\) in \(L_1\) norm for triangular matrices. The matrix has to have the LAPACKSupport::Property set to either LAPACKSupport::Property::upper_triangular or LAPACKSupport::Property::lower_triangular, see set_property().

Definition at line 1481 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::determinant ( ) const

Compute the determinant of a matrix. As it requires the LU factorization of the matrix, this function can only be called after compute_lu_factorization() has been called.

Definition at line 1800 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::l1_norm ( ) const

Compute \(L_1\) norm.

Definition at line 1340 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::linfty_norm ( ) const

Compute \(L_\infty\) norm.

Definition at line 1350 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::frobenius_norm ( ) const

Compute Frobenius norm

Definition at line 1360 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::trace ( ) const

Compute trace of the matrix, i.e. the sum of the diagonal values. Obviously, the matrix needs to be quadratic for this function.

Definition at line 1403 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::invert ( )

Invert the matrix by first computing an LU/Cholesky factorization with the LAPACK function Xgetrf/Xpotrf and then building the actual inverse using Xgetri/Xpotri.

Definition at line 1639 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::solve ( Vector< number > &  v,
const bool  transposed = false 
) const

Solve the linear system with right hand side v and put the solution back to v. The matrix should be either triangular or LU/Cholesky factorization should be previously computed.

The flag transposed indicates whether the solution of the transposed system is to be performed.

Definition at line 1682 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::solve ( LAPACKFullMatrix< number > &  B,
const bool  transposed = false 
) const

Same as above but for multiple right hand sides (as many as there are columns in the matrix B).

Definition at line 1725 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::apply_lu_factorization ( Vector< number > &  v,
const bool  transposed 
) const

Solve the linear system with right hand side given by applying forward/backward substitution to the previously computed LU factorization. Uses LAPACK function Xgetrs.

The flag transposed indicates whether the solution of the transposed system is to be performed.

Deprecated:
use solve() instead.

Definition at line 1780 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::apply_lu_factorization ( LAPACKFullMatrix< number > &  B,
const bool  transposed 
) const

Solve the linear system with multiple right hand sides (as many as there are columns in the matrix b) given by applying forward/backward substitution to the previously computed LU factorization. Uses LAPACK function Xgetrs.

The flag transposed indicates whether the solution of the transposed system is to be performed.

Deprecated:
use solve() instead.

Definition at line 1790 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_eigenvalues ( const bool  right_eigenvectors = false,
const bool  left_eigenvectors = false 
)

Compute eigenvalues of the matrix. After this routine has been called, eigenvalues can be retrieved using the eigenvalue() function. The matrix itself will be LAPACKSupport::unusable after this operation.

The optional arguments allow to compute left and right eigenvectors as well.

Note that the function does not return the computed eigenvalues right away since that involves copying data around between the output arrays of the LAPACK functions and any return array. This is often unnecessary since one may not be interested in all eigenvalues at once, but for example only the extreme ones. In that case, it is cheaper to just have this function compute the eigenvalues and have a separate function that returns whatever eigenvalue is requested.

Note
Calls the LAPACK function Xgeev.

Definition at line 1831 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::compute_eigenvalues_symmetric ( const number  lower_bound,
const number  upper_bound,
const number  abs_accuracy,
Vector< number > &  eigenvalues,
FullMatrix< number > &  eigenvectors 
)

Compute eigenvalues and eigenvectors of a real symmetric matrix. Only eigenvalues in the interval \((\rm{lower\_bound}, \rm{upper\_bound}]\) are computed with the absolute tolerance \(\rm abs\_accuracy\). An approximate eigenvalue is accepted as converged when it is determined to lie in an interval \([a,b]\) of width less than or equal to \(\rm{abs\_accuracy} + eps * \rm{max}(|a|,|b|)\), where \(eps\) is the machine precision. If \(\rm{abs\_accuracy}\) is less than or equal to zero, then \(eps\,|\mathbf{T}|_1\) will be used in its place, where \(|\mathbf{T}|_1\) is the 1-norm of the tridiagonal matrix obtained by reducing \(\mathbf A\) to tridiagonal form. Eigenvalues will be computed most accurately when \(\rm{abs\_accuracy}\) is set to twice the underflow threshold, not zero. After this routine has been called, all eigenvalues in \((\rm{lower\_bound}, \rm{upper\_bound}]\) will be stored in eigenvalues and the corresponding eigenvectors will be stored in the columns of eigenvectors, whose dimension is set accordingly.

Note
Calls the LAPACK function Xsyevx.

Definition at line 1911 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::compute_generalized_eigenvalues_symmetric ( LAPACKFullMatrix< number > &  B,
const number  lower_bound,
const number  upper_bound,
const number  abs_accuracy,
Vector< number > &  eigenvalues,
std::vector< Vector< number >> &  eigenvectors,
const types::blas_int  itype = 1 
)

Compute generalized eigenvalues and eigenvectors of a real generalized symmetric eigenproblem of the form

  • itype = 1: \(\mathbf A \cdot \mathbf x=\lambda \mathbf B \cdot \mathbf x\)
  • itype = 2: \(\mathbf A \cdot \mathbf B \cdot \mathbf x=\lambda \mathbf x\)
  • itype = 3: \(\mathbf B \cdot \mathbf A \cdot \mathbf x=\lambda \mathbf x\)

where \(\mathbf A\) is this matrix. \(\mathbf A\) and \(\mathbf B\) are assumed to be symmetric, and \(\mathbf B\) has to be positive definite. Only eigenvalues in the interval \((\rm{lower\_bound}, \rm{upper\_bound}]\) are computed with the absolute tolerance \(\rm{abs\_accuracy}\). An approximate eigenvalue is accepted as converged when it is determined to lie in an interval \([a,b]\) of width less than or equal to \(\rm{abs\_accuracy} + eps * \rm{max}( |a|,|b| )\), where \(eps\) is the machine precision. If \(\rm{abs\_accuracy}\) is less than or equal to zero, then \(eps \, |\mathbf{T}|_1\) will be used in its place, where \(|\mathbf{T}|_1\) is the 1-norm of the tridiagonal matrix obtained by reducing \(\mathbf A\) to tridiagonal form. Eigenvalues will be computed most accurately when \(\rm{abs\_accuracy}\) is set to twice the underflow threshold, not zero. After this routine has been called, all eigenvalues in \((\rm{lower\_bound}, \rm{upper\_bound}]\) will be stored in eigenvalues and the corresponding eigenvectors will be stored in eigenvectors, whose dimension is set accordingly.

Note
Calls the LAPACK function Xsygvx.

Definition at line 2022 of file lapack_full_matrix.cc.

template<typename number>
void LAPACKFullMatrix< number >::compute_generalized_eigenvalues_symmetric ( LAPACKFullMatrix< number > &  B,
std::vector< Vector< number >> &  eigenvectors,
const types::blas_int  itype = 1 
)

Same as the other compute_generalized_eigenvalues_symmetric function except that all eigenvalues are computed and the tolerance is set automatically. Note that this function does not return the computed eigenvalues right away since that involves copying data around between the output arrays of the LAPACK functions and any return array. This is often unnecessary since one may not be interested in all eigenvalues at once, but for example only the extreme ones. In that case, it is cheaper to just have this function compute the eigenvalues and have a separate function that returns whatever eigenvalue is requested. Eigenvalues can be retrieved using the eigenvalue() function. The number of computed eigenvectors is equal to eigenvectors.size()

Note
Calls the LAPACK function Xsygv.

Definition at line 2147 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_svd ( )

Compute the singular value decomposition of the matrix using LAPACK function Xgesdd.

Requires that the state is LAPACKSupport::matrix, fills the data members wr, svd_u, and svd_vt, and leaves the object in the state LAPACKSupport::svd.

Definition at line 1519 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_inverse_svd ( const double  threshold = 0.)

Compute the inverse of the matrix by singular value decomposition.

Requires that state is either LAPACKSupport::matrix or LAPACKSupport::svd. In the first case, this function calls compute_svd(). After this function, the object will have the state LAPACKSupport::inverse_svd.

For a singular value decomposition, the inverse is simply computed by replacing all singular values by their reciprocal values. If the matrix does not have maximal rank, singular values 0 are not touched, thus computing the minimal norm right inverse of the matrix.

The parameter threshold determines, when a singular value should be considered zero. It is the ratio of the smallest to the largest nonzero singular value \(s_{max}\). Thus, the inverses of all singular values less than \(s_{max}/\rm{threshold}\) will be set to zero.

Definition at line 1595 of file lapack_full_matrix.cc.

template<typename number >
void LAPACKFullMatrix< number >::compute_inverse_svd_with_kernel ( const unsigned int  kernel_size)

Same as above but provide the size of the kernel instead of a threshold, i.e. the kernel_size smallest eigenvalues.

Definition at line 1618 of file lapack_full_matrix.cc.

template<typename number >
std::complex< number > LAPACKFullMatrix< number >::eigenvalue ( const size_type  i) const
inline

Retrieve eigenvalue after compute_eigenvalues() was called.

Definition at line 1116 of file lapack_full_matrix.h.

template<typename number >
number LAPACKFullMatrix< number >::singular_value ( const size_type  i) const
inline

Retrieve singular values after compute_svd() or compute_inverse_svd() was called.

Definition at line 1132 of file lapack_full_matrix.h.

template<typename number >
void LAPACKFullMatrix< number >::print_formatted ( std::ostream &  out,
const unsigned int  precision = 3,
const bool  scientific = true,
const unsigned int  width = 0,
const char *  zero_string = " ",
const double  denominator = 1.,
const double  threshold = 0. 
) const

Print the matrix and allow formatting of entries.

The parameters allow for a flexible setting of the output format:

Parameters
precisiondenotes the number of trailing digits.
scientificis used to determine the number format, where scientific = false means fixed point notation.
widthdenotes the with of each column. A zero entry for width makes the function compute a width, but it may be changed to a positive value, if output is crude.
zero_stringspecifies a string printed for zero entries.
denominatorMultiply the whole matrix by this common denominator to get nicer numbers.
thresholdall entries with absolute value smaller than this are considered zero.
Note
The entries stored resemble a matrix only if the state is either LAPACKSupport::matrix or LAPACK::inverse_matrix. Otherwise, calling this function is not allowed.

Definition at line 2239 of file lapack_full_matrix.cc.

template<typename number >
number LAPACKFullMatrix< number >::norm ( const char  type) const
private

Internal function to compute various norms.

Definition at line 1370 of file lapack_full_matrix.cc.

Member Data Documentation

template<typename number>
LAPACKSupport::State LAPACKFullMatrix< number >::state
private

Since LAPACK operations notoriously change the meaning of the matrix entries, we record the current state after the last operation here.

Definition at line 878 of file lapack_full_matrix.h.

template<typename number>
LAPACKSupport::Property LAPACKFullMatrix< number >::property
private

Additional property of the matrix which may help to select more efficient LAPACK functions.

Definition at line 884 of file lapack_full_matrix.h.

template<typename number>
std::vector<number> LAPACKFullMatrix< number >::work
mutableprivate

The working array used for some LAPACK functions.

Definition at line 889 of file lapack_full_matrix.h.

template<typename number>
std::vector<types::blas_int> LAPACKFullMatrix< number >::iwork
mutableprivate

Integer working array used for some LAPACK functions.

Definition at line 894 of file lapack_full_matrix.h.

template<typename number>
std::vector<types::blas_int> LAPACKFullMatrix< number >::ipiv
private

The vector storing the permutations applied for pivoting in the LU- factorization.

Also used as the scratch array IWORK for LAPACK functions needing it.

Definition at line 902 of file lapack_full_matrix.h.

template<typename number>
std::vector<number> LAPACKFullMatrix< number >::inv_work
private

Workspace for calculating the inverse matrix from an LU factorization.

Definition at line 907 of file lapack_full_matrix.h.

template<typename number>
std::vector<typename numbers::NumberTraits<number>::real_type> LAPACKFullMatrix< number >::wr
private

Real parts of eigenvalues or the singular values. Filled by compute_eigenvalues() or compute_svd().

Definition at line 913 of file lapack_full_matrix.h.

template<typename number>
std::vector<number> LAPACKFullMatrix< number >::wi
private

Imaginary parts of eigenvalues, or, in the complex scalar case, the eigenvalues themselves. Filled by compute_eigenvalues.

Definition at line 919 of file lapack_full_matrix.h.

template<typename number>
std::vector<number> LAPACKFullMatrix< number >::vl
private

Space where left eigenvectors can be stored.

Definition at line 924 of file lapack_full_matrix.h.

template<typename number>
std::vector<number> LAPACKFullMatrix< number >::vr
private

Space where right eigenvectors can be stored.

Definition at line 929 of file lapack_full_matrix.h.

template<typename number>
std::unique_ptr<LAPACKFullMatrix<number> > LAPACKFullMatrix< number >::svd_u
private

The matrix \(\mathbf U\) in the singular value decomposition \(\mathbf U \cdot \mathbf S \cdot \mathbf V^T\).

Definition at line 935 of file lapack_full_matrix.h.

template<typename number>
std::unique_ptr<LAPACKFullMatrix<number> > LAPACKFullMatrix< number >::svd_vt
private

The matrix \(\mathbf V^T\) in the singular value decomposition \(\mathbf U \cdot \mathbf S \cdot \mathbf V^T\).

Definition at line 941 of file lapack_full_matrix.h.

template<typename number>
Threads::Mutex LAPACKFullMatrix< number >::mutex
mutableprivate

Thread mutex.

Definition at line 946 of file lapack_full_matrix.h.


The documentation for this class was generated from the following files: