Classes | Functions
o2scl_linalg Namespace Reference

The namespace for linear algebra classes and functions. More...

Classes

class  lanczos
 Lanczos diagonalization. More...
 
class  linear_solver
 A generic solver for the linear system $ A x = b $ [abstract base]. More...
 
class  linear_solver_arma
 Armadillo linear solver. More...
 
class  linear_solver_eigen_colQR
 Eigen linear solver using QR decomposition with column pivoting. More...
 
class  linear_solver_eigen_fullLU
 Eigen linear solver using LU decomposition with full pivoting. More...
 
class  linear_solver_eigen_fullQR
 Eigen linear solver using QR decomposition with full pivoting. More...
 
class  linear_solver_eigen_houseQR
 Eigen linear solver using QR decomposition with column pivoting. More...
 
class  linear_solver_eigen_LDLT
 Eigen linear solver using LDLT decomposition with full pivoting. More...
 
class  linear_solver_eigen_LLT
 Eigen linear solver using LLT decomposition with full pivoting. More...
 
class  linear_solver_eigen_partLU
 Eigen linear solver using LU decomposition with partial pivoting. More...
 
class  linear_solver_HH
 Generic Householder linear solver. More...
 
class  linear_solver_LU
 Generic linear solver using LU decomposition. More...
 
class  linear_solver_QR
 Generic linear solver using QR decomposition. More...
 
class  ubvector_2_mem
 Allocation object for 2 arrays of equal size. More...
 
class  ubvector_4_mem
 Allocation object for 4 arrays of equal size. More...
 
class  ubvector_5_mem
 Allocation object for 5 arrays of equal size. More...
 

Functions

template<class mat_t , class vec_t , class vec2_t >
int bidiag_decomp (size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V)
 Factor a matrix into bidiagonal form. More...
 
template<class mat_t , class vec_t , class mat2_t , class vec2_t , class mat3_t , class vec3_t , class vec4_t >
int bidiag_unpack (size_t M, size_t N, const mat_t &A, const vec_t &tau_U, mat2_t &U, const vec2_t &tau_V, mat3_t &V, vec3_t &diag, vec4_t &superdiag)
 Unpack a matrix A with the bidiagonal decomposition and create matrices U, V, diagonal diag and superdiagonal superdiag. More...
 
template<class mat_t , class vec_t , class vec2_t , class mat2_t >
int bidiag_unpack2 (size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V, mat2_t &V)
 Unpack a matrix A with the bidiagonal decomposition and create matrix V.
 
template<class mat_t , class vec_t , class vec2_t >
int bidiag_unpack_B (size_t M, size_t N, const mat_t &A, vec_t &diag, vec2_t &superdiag)
 Unpack the diagonal and superdiagonal of the bidiagonal decomposition of A into diag and superdiag.
 
template<>
void cholesky_decomp< Eigen::MatrixXd > (const size_t M, Eigen::MatrixXd &A)
 Eigen specialization of cholesky_decomp()
 
template<class mat_t >
void cholesky_decomp (const size_t M, mat_t &A)
 Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix. More...
 
template<class mat_t , class vec_t , class vec2_t >
void cholesky_solve (const size_t N, const mat_t &LLT, const vec_t &b, vec2_t &x)
 Solve a symmetric positive-definite linear system after a Cholesky decomposition. More...
 
template<class mat_t , class vec_t >
void cholesky_svx (const size_t N, const mat_t &LLT, vec_t &x)
 Solve a linear system in place using a Cholesky decomposition.
 
template<class mat_t >
void cholesky_invert (const size_t N, mat_t &LLT)
 Compute the inverse of a symmetric positive definite matrix given the Cholesky decomposition. More...
 
template<class mat_t , class vec_t >
int cholesky_decomp_unit (const size_t N, mat_t &A, vec_t &D)
 Cholesky decomposition with unit-diagonal triangular parts. More...
 
void create_givens (const double a, const double b, double &c, double &s)
 Create a Givens rotation matrix. More...
 
template<class mat1_t , class mat2_t >
void apply_givens_qr (size_t M, size_t N, mat1_t &Q, mat2_t &R, size_t i, size_t j, double c, double s)
 Apply a rotation to matrices from the QR decomposition. More...
 
template<class mat1_t , class mat2_t >
void apply_givens_lq (size_t M, size_t N, mat1_t &Q, mat2_t &L, size_t i, size_t j, double c, double s)
 Apply a rotation to matrices from the LQ decomposition. More...
 
template<class vec_t >
void apply_givens_vec (vec_t &v, size_t i, size_t j, double c, double s)
 Apply a rotation to a vector, $ v \rightarrow G^{T} v $.
 
template<class mat_t , class vec_t >
int HH_svx (size_t N, size_t M, mat_t &A, vec_t &x)
 Solve a linear system after Householder decomposition in place. More...
 
template<class mat_t , class vec_t , class vec2_t >
int HH_solve (size_t n, mat_t &A, const vec_t &b, vec2_t &x)
 Solve linear system after Householder decomposition.
 
template<class vec_t >
double householder_transform (const size_t n, vec_t &v)
 Replace the vector v with a Householder vector and a coefficient tau that annihilates v[1] through v[n-1] (inclusive) More...
 
template<class mat_t >
double householder_transform_subcol (mat_t &A, const size_t ir, const size_t ic, const size_t M)
 Compute the Householder transform of a vector formed with n rows of a column of a matrix. More...
 
template<class mat_t >
double householder_transform_subrow (mat_t &A, const size_t ir, const size_t ic, const size_t N)
 Compute the Householder transform of a vector formed with the last n columns of a row of a matrix. More...
 
template<class vec_t , class mat_t >
void householder_hm (const size_t M, const size_t N, double tau, const vec_t &v, mat_t &A)
 Apply a Householder transformation $ (v,\tau) $ to matrix $ A $ of size (M,N) More...
 
template<class mat_t >
void householder_hm_subcol (mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat_t &M2, const size_t ir2, const size_t ic2, double tau)
 Apply a Householder transformation to the lower-right part of M when the transformation is stored in a column of M2. More...
 
template<class mat_t >
void householder_hm_subrow (mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat_t &M2, const size_t ir2, const size_t ic2, double tau)
 Apply a Householder transformation to the lower-right part of M when the transformation is stored in a row of M2. More...
 
template<class vec_t , class vec2_t >
void householder_hv (const size_t N, double tau, const vec_t &v, vec2_t &w)
 Apply a Householder transformation v to vector w.
 
template<class mat_t , class vec_t >
void householder_hv_subcol (const mat_t &A, vec_t &w, double tau, const size_t ie, const size_t N)
 Apply a Householder transformation v to vector w where v is stored as a column in a matrix A. More...
 
template<class mat_t >
void householder_hm1 (const size_t M, const size_t N, double tau, mat_t &A)
 Apply a Householder transformation $ (v,\tau) $ to a matrix being build up from the identity matrix, using the first column of A as a Householder vector.
 
template<class mat_t >
void householder_hm1_sub (const size_t M, const size_t N, double tau, mat_t &A, size_t ir, size_t ic)
 Apply a Householder transformation $ (v,\tau) $ to a matrix being build up from the identity matrix, using the first column of A as a Householder vector. More...
 
template<class vec_t , class mat_t >
void householder_mh (const size_t M, const size_t N, double tau, const vec_t &v, mat_t &A)
 Apply the Householder transformation (v,tau) to the right-hand side of the matrix A.
 
template<class mat_t , class mat2_t >
void householder_mh_subrow (mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat2_t &M2, const size_t ir2, const size_t ic2, double tau)
 Apply the Householder transformation (v,tau) to the right-hand side of the matrix A. More...
 
template<class mat_t >
int diagonal_has_zero (const size_t N, mat_t &A)
 Return 1 if at least one of the elements in the diagonal is zero.
 
template<class mat_t >
int LU_decomp (const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
 Compute the LU decomposition of the matrix A. More...
 
template<class mat_t , class vec_t >
int LU_svx (const size_t N, const mat_t &LU, const o2scl::permutation &p, vec_t &x)
 Solve a linear system after LU decomposition in place. More...
 
template<class mat_t , class mat_row_t >
int LU_decomp_alt (const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
 An alternate form of LU decomposition with matrix row objects.
 
template<class mat_t , class vec_t , class vec2_t >
int LU_solve (const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
 Solve a linear system after LU decomposition. More...
 
template<class mat_t , class mat2_t , class vec_t , class vec2_t , class vec3_t >
int LU_refine (const size_t N, const mat_t &A, const mat2_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x, vec3_t &residual)
 Refine the solution of a linear system. More...
 
template<class mat_t , class mat2_t , class mat_col_t >
int LU_invert (const size_t N, const mat_t &LU, const o2scl::permutation &p, mat2_t &inverse)
 Compute the inverse of a matrix from its LU decomposition. More...
 
template<class mat_t >
double LU_det (const size_t N, const mat_t &LU, int signum)
 Compute the determinant of a matrix from its LU decomposition. More...
 
template<class mat_t >
double LU_lndet (const size_t N, const mat_t &LU)
 Compute the logarithm of the absolute value of the determinant of a matrix from its LU decomposition. More...
 
template<class mat_t >
int LU_sgndet (const size_t N, const mat_t &LU, int signum)
 Compute the sign of the determinant of a matrix from its LU decomposition. More...
 
template<>
void QR_decomp_unpack< arma::mat, arma::mat, arma::mat > (const size_t M, const size_t N, arma::mat &A, arma::mat &Q, arma::mat &R)
 Armadillo specialization of QR_decomp_unpack().
 
template<>
void QR_decomp_unpack< Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd > (const size_t M, const size_t N, Eigen::MatrixXd &A, Eigen::MatrixXd &Q, Eigen::MatrixXd &R)
 Eigen specialization of QR_decomp_unpack().
 
template<class mat_t , class vec_t >
void QR_decomp (size_t M, size_t N, mat_t &A, vec_t &tau)
 Compute the QR decomposition of matrix A.
 
template<class mat_t , class vec_t , class vec2_t >
void QR_QTvec (const size_t M, const size_t N, const mat_t &QR, const vec_t &tau, vec2_t &v)
 Form the product Q^T v from a QR factorized matrix.
 
template<class mat1_t , class mat2_t , class mat3_t , class vec_t >
void QR_unpack (const size_t M, const size_t N, const mat1_t &QR, const vec_t &tau, mat2_t &Q, mat3_t &R)
 Unpack the QR matrix to the individual Q and R components.
 
template<class mat_t , class vec_t , class vec2_t >
void QR_svx (size_t M, size_t N, const mat_t &QR, const vec_t &tau, vec2_t &x)
 Solve the system A x = b in place using the QR factorization.
 
template<class mat_t , class vec_t , class vec2_t , class vec3_t >
void QR_solve (size_t N, const mat_t &QR, const vec_t &tau, const vec2_t &b, vec3_t &x)
 Solve the system A x = b using the QR factorization.
 
template<class mat1_t , class mat2_t , class vec1_t , class vec2_t >
void QR_update (size_t M, size_t N, mat1_t &Q, mat2_t &R, vec1_t &w, vec2_t &v)
 Update a QR factorisation for A= Q R, A' = A + u v^T,. More...
 
template<class mat_t , class mat2_t , class mat3_t >
void QR_decomp_unpack (const size_t M, const size_t N, mat_t &A, mat2_t &Q, mat3_t &R)
 Compute the unpacked QR decomposition of matrix A. More...
 
template<class mat_t , class vec_t , class vec2_t >
void QRPT_decomp (size_t M, size_t N, mat_t &A, vec_t &tau, o2scl::permutation &p, int &signum, vec2_t &norm)
 Compute the QR decomposition of matrix A.
 
template<class mat_t , class mat2_t , class vec_t , class vec2_t >
void SV_decomp (size_t M, size_t N, mat_t &A, mat2_t &V, vec_t &S, vec2_t &work)
 Factorise a general matrix into its SV decomposition using the Golub-Reinsch algorithm. More...
 
template<class mat_t , class mat2_t , class mat3_t , class vec_t , class vec2_t >
void SV_decomp_mod (size_t M, size_t N, mat_t &A, mat2_t &X, mat3_t &V, vec_t &S, vec2_t &work)
 SV decomposition by the modified Golub-Reinsch algorithm which is better for $ M \gg N $. More...
 
template<class mat_t , class mat2_t , class vec_t , class vec2_t , class vec3_t >
void SV_solve (size_t M, size_t N, mat_t &U, mat2_t &V, vec_t &S, vec2_t &b, vec3_t &x)
 Solve the system A x = b using the SV decomposition. More...
 
template<class mat_t , class mat2_t , class vec_t >
void SV_decomp_jacobi (size_t M, size_t N, mat_t &A, mat2_t &Q, vec_t &S)
 SV decomposition using one-sided Jacobi orthogonalization. More...
 
template<class mat_t , class vec_t >
void balance_columns (size_t M, size_t N, mat_t &A, vec_t &D)
 Balance a general matrix A by scaling the columns by the diagonal matrix D. More...
 
template<class vec_t , class vec2_t >
void chop_small_elements (size_t N, vec_t &d, vec2_t &f)
 Zero out small elements in f according to the scales set in d. More...
 
template<class vec_t , class vec2_t >
double trailing_eigenvalue (size_t n, const vec_t &d, const vec2_t &f)
 Desc. More...
 
void create_schur (double d0, double f0, double d1, double &c, double &s)
 Desc. More...
 
template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void svd2 (size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
 2-variable SVD More...
 
template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void svd2_sub (size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a)
 Shifted 2-variable SVD. More...
 
template<class vec_t , class vec2_t , class mat_t >
void chase_out_intermediate_zero (size_t M, size_t n, vec_t &d, vec2_t &f, mat_t &U, size_t k0)
 Desc. More...
 
template<class vec_t , class vec2_t , class mat_t >
void chase_out_trailing_zero (size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &V)
 Desc. More...
 
template<class vec_t , class vec2_t , class mat_t >
void chase_out_trailing_zero_sub (size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &V, size_t a)
 Desc. More...
 
template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void qrstep (size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
 Desc. More...
 
template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void qrstep_sub (size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a)
 A special form of qrstep() for SV_decomp() More...
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class mem_t , class mem_vec_t >
void solve_tridiag_sym (const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
 Solve a symmetric tridiagonal linear system with user-specified memory. More...
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t , class mem_t , class mem_vec_t >
void solve_tridiag_nonsym (const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
 Solve an asymmetric tridiagonal linear system with user-specified memory. More...
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class mem_t , class mem_vec_t >
void solve_cyc_tridiag_sym (const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
 Solve a symmetric cyclic tridiagonal linear system with user specified memory. More...
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t , class mem_t , class mem_vec_t >
void solve_cyc_tridiag_nonsym (const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
 Solve an asymmetric cyclic tridiagonal linear system with user-specified memory. More...
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t >
void solve_tridiag_sym (const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N)
 Solve a symmetric tridiagonal linear system.
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t >
void solve_tridiag_nonsym (const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N)
 Solve an asymmetric tridiagonal linear system.
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t >
void solve_cyc_tridiag_sym (const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N)
 Solve a symmetric cyclic tridiagonal linear system.
 
template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t >
void solve_cyc_tridiag_nonsym (const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N)
 Solve an asymmetric cyclic tridiagonal linear system.
 

Detailed Description

See Linear Algebra for more complete information about linear algebra in O2scl .

This namespace documentation is in the file src/base/cblas.h

Function Documentation

◆ apply_givens_lq()

template<class mat1_t , class mat2_t >
void o2scl_linalg::apply_givens_lq ( size_t  M,
size_t  N,
mat1_t &  Q,
mat2_t &  L,
size_t  i,
size_t  j,
double  c,
double  s 
)

This performs $ Q \rightarrow Q G $ and $ L \rightarrow L G^{T} $.

Definition at line 86 of file givens_base.h.

◆ apply_givens_qr()

template<class mat1_t , class mat2_t >
void o2scl_linalg::apply_givens_qr ( size_t  M,
size_t  N,
mat1_t &  Q,
mat2_t &  R,
size_t  i,
size_t  j,
double  c,
double  s 
)

This performs $ Q \rightarrow Q G $ and $ R \rightarrow G^{T} R $.

Definition at line 58 of file givens_base.h.

◆ balance_columns()

template<class mat_t , class vec_t >
void o2scl_linalg::balance_columns ( size_t  M,
size_t  N,
mat_t &  A,
vec_t &  D 
)

The function computes $ A_{\mathrm{new}} = A D^{-1} $ where $ D $ is a diagonal matrix.

Definition at line 615 of file svd_base.h.

◆ bidiag_decomp()

template<class mat_t , class vec_t , class vec2_t >
int o2scl_linalg::bidiag_decomp ( size_t  M,
size_t  N,
mat_t &  A,
vec_t &  tau_U,
vec2_t &  tau_V 
)

Factor matrix A of size (M,N) with $ M\geq N $ into $ A = U B V^T $ where U and V are orthogonal and B is upper bidiagonal.

After the function call, the matrix $ B $ is stored the diagonal and first superdiagonal of A. The matrices $ U $ and $ V $ are stored as packed sets of Householder transformations in the lower and upper triangular parts of A, respectively.

Adapted from the GSL version which was based on algorithm 5.4.2 in Golub96.

Definition at line 65 of file bidiag_base.h.

◆ bidiag_unpack()

template<class mat_t , class vec_t , class mat2_t , class vec2_t , class mat3_t , class vec3_t , class vec4_t >
int o2scl_linalg::bidiag_unpack ( size_t  M,
size_t  N,
const mat_t &  A,
const vec_t &  tau_U,
mat2_t &  U,
const vec2_t &  tau_V,
mat3_t &  V,
vec3_t &  diag,
vec4_t &  superdiag 
)

Given a matrix A of size (M,N) with $ M \geq N $ created by bidiag_decomp(), this function creates the matrix U of size (M,N), the matrix V of size (N,N), the diagonal diag of size N and the super-diagonal superdiag of size N-1.

Definition at line 113 of file bidiag_base.h.

◆ chase_out_intermediate_zero()

template<class vec_t , class vec2_t , class mat_t >
void o2scl_linalg::chase_out_intermediate_zero ( size_t  M,
size_t  n,
vec_t &  d,
vec2_t &  f,
mat_t &  U,
size_t  k0 
)

The vector d should be of size n, the vector f should be of size n, and the matrix U should be of size (M,n)

Used in qrstep() and qrstep_sub().

Definition at line 457 of file svdstep_base.h.

◆ chase_out_trailing_zero()

template<class vec_t , class vec2_t , class mat_t >
void o2scl_linalg::chase_out_trailing_zero ( size_t  N,
size_t  n,
vec_t &  d,
vec2_t &  f,
mat_t &  V 
)

The vector d should be of size n, the vector f should be of size n, and the matrix V should be of size (N,n)

Used in qrstep().

Definition at line 507 of file svdstep_base.h.

◆ chase_out_trailing_zero_sub()

template<class vec_t , class vec2_t , class mat_t >
void o2scl_linalg::chase_out_trailing_zero_sub ( size_t  N,
size_t  n,
vec_t &  d,
vec2_t &  f,
mat_t &  V,
size_t  a 
)

The vector d should be of size n, the vector f should be of size n, and the matrix V should be of size (N,n)

Used in qrstep_sub().

Definition at line 556 of file svdstep_base.h.

◆ cholesky_decomp()

template<class mat_t >
void o2scl_linalg::cholesky_decomp ( const size_t  M,
mat_t &  A 
)

On input, the upper triangular part of A is ignored (only the lower triangular part and diagonal are used). On output, the diagonal and lower triangular part contain the matrix L and the upper triangular part contains L^T.

If the matrix is not positive-definite, the error handler will be called.

Definition at line 61 of file cholesky_base.h.

◆ cholesky_decomp_unit()

template<class mat_t , class vec_t >
int o2scl_linalg::cholesky_decomp_unit ( const size_t  N,
mat_t &  A,
vec_t &  D 
)

A = L D L^T, where diag(L) = (1,1,...,1). Upon exit, A contains L and L^T as for Cholesky, and the diagonal of A is (1,1,...,1). The vector Dis set to the diagonal elements of the diagonal matrix D.

Definition at line 299 of file cholesky_base.h.

◆ cholesky_invert()

template<class mat_t >
void o2scl_linalg::cholesky_invert ( const size_t  N,
mat_t &  LLT 
)

Given a Cholesky decomposition produced by cholesky_decomp(), this function returns the inverse of that matrix in LLT.

Definition at line 204 of file cholesky_base.h.

◆ cholesky_solve()

template<class mat_t , class vec_t , class vec2_t >
void o2scl_linalg::cholesky_solve ( const size_t  N,
const mat_t &  LLT,
const vec_t &  b,
vec2_t &  x 
)

Given the Cholesky decomposition of a matrix A in LLT, this function solves the system A*x=b.

Definition at line 157 of file cholesky_base.h.

◆ chop_small_elements()

template<class vec_t , class vec2_t >
void o2scl_linalg::chop_small_elements ( size_t  N,
vec_t &  d,
vec2_t &  f 
)

The parameter N is the size of d.

Used in SV_decomp().

Definition at line 58 of file svdstep_base.h.

◆ create_givens()

void o2scl_linalg::create_givens ( const double  a,
const double  b,
double &  c,
double &  s 
)

Given values a and b, create entries c and s of a matrix for which

\[ \left[ \begin{array}{cc} c & -s \\ s & c \end{array} \right] \left[ \begin{array}{c} a \\ b \end{array} \right] = \left[ \begin{array}{c} r \\ 0 \end{array} \right] \]

◆ create_schur()

void o2scl_linalg::create_schur ( double  d0,
double  f0,
double  d1,
double &  c,
double &  s 
)

Used in svd2() and svd2_sub().

Definition at line 120 of file svdstep_base.h.

◆ HH_svx()

template<class mat_t , class vec_t >
int o2scl_linalg::HH_svx ( size_t  N,
size_t  M,
mat_t &  A,
vec_t &  x 
)
Idea for Future:
Handle memory allocation like the tridiagonal functions.

Definition at line 59 of file hh_base.h.

◆ householder_hm()

template<class vec_t , class mat_t >
void o2scl_linalg::householder_hm ( const size_t  M,
const size_t  N,
double  tau,
const vec_t &  v,
mat_t &  A 
)

The vector v must have at least N entries, with the exception that the vector element v[0] is never referenced by this function.

Definition at line 210 of file householder_base.h.

◆ householder_hm1_sub()

template<class mat_t >
void o2scl_linalg::householder_hm1_sub ( const size_t  M,
const size_t  N,
double  tau,
mat_t &  A,
size_t  ir,
size_t  ic 
)

Used in SV_decomp_mod() and bidiag_unpack2().

Definition at line 436 of file householder_base.h.

◆ householder_hm_subcol()

template<class mat_t >
void o2scl_linalg::householder_hm_subcol ( mat_t &  M,
const size_t  ir,
const size_t  ic,
const size_t  nr,
const size_t  nc,
const mat_t &  M2,
const size_t  ir2,
const size_t  ic2,
double  tau 
)

This applies a householder transformation (v,tau) to a lower-right submatrix of M. The submatrix has nr-ir rows and nc-ic columns and starts at row ir of column ic of the original matrix M. The vector containing the transformation is taken from a column of M2 starting at row ir2 and column ic2. The matrix M2 must have at least ic2+1 columns and at least nr-ir+ir2 rows.

This function is used in QR_decomp() and QR_unpack() .

Definition at line 255 of file householder_base.h.

◆ householder_hm_subrow()

template<class mat_t >
void o2scl_linalg::householder_hm_subrow ( mat_t &  M,
const size_t  ir,
const size_t  ic,
const size_t  nr,
const size_t  nc,
const mat_t &  M2,
const size_t  ir2,
const size_t  ic2,
double  tau 
)

This applies a householder transformation (v,tau) to a lower-right submatrix of M. The submatrix has nr-ir rows and nc-ic columns and starts at row ir of column ic of the original matrix M. The vector containing the transformation is taken from a row of M2 starting at row ir2 and column ic2. The matrix M2 must have at least ir2+1 rows and nr-ir+ic2 columns.

Used in bidiag_unpack().

Definition at line 301 of file householder_base.h.

◆ householder_hv_subcol()

template<class mat_t , class vec_t >
void o2scl_linalg::householder_hv_subcol ( const mat_t &  A,
vec_t &  w,
double  tau,
const size_t  ie,
const size_t  N 
)

Used in QR_QTvec().

Definition at line 364 of file householder_base.h.

◆ householder_mh_subrow()

template<class mat_t , class mat2_t >
void o2scl_linalg::householder_mh_subrow ( mat_t &  M,
const size_t  ir,
const size_t  ic,
const size_t  nr,
const size_t  nc,
const mat2_t &  M2,
const size_t  ir2,
const size_t  ic2,
double  tau 
)

Used in bidiag_decomp().

Definition at line 519 of file householder_base.h.

◆ householder_transform()

template<class vec_t >
double o2scl_linalg::householder_transform ( const size_t  n,
vec_t &  v 
)

On exit, this function returns the value of $ \tau = 2/ (v^{T} v) $. If n is less than or equal to 1 then this function returns zero without calling the error handler.

Definition at line 68 of file householder_base.h.

◆ householder_transform_subcol()

template<class mat_t >
double o2scl_linalg::householder_transform_subcol ( mat_t &  A,
const size_t  ir,
const size_t  ic,
const size_t  M 
)

This performs a Householder transform of a vector defined by a column of a matrix A which starts at element A(ir,ic) and ends at element A(M-1,ic). If M-1 is equal to ir+1, this function quietly does nothing.

Used in QR_decomp() and SV_decomp_mod().

Definition at line 118 of file householder_base.h.

◆ householder_transform_subrow()

template<class mat_t >
double o2scl_linalg::householder_transform_subrow ( mat_t &  A,
const size_t  ir,
const size_t  ic,
const size_t  N 
)

This performs a Householder transform of a vector defined by a row of a matrix A which starts at element A(ir,ic) and ends at element A(ir,N-1) If N-1 is equal to ic, this function quietly does nothing.

Definition at line 165 of file householder_base.h.

◆ LU_decomp()

template<class mat_t >
int o2scl_linalg::LU_decomp ( const size_t  N,
mat_t &  A,
o2scl::permutation p,
int &  signum 
)

On output the diagonal and upper triangular part of the input matrix A contain the matrix U. The lower triangular part of the input matrix (excluding the diagonal) contains L. The diagonal elements of L are unity, and are not stored.

The permutation matrix P is encoded in the permutation p. The j-th column of the matrix P is given by the k-th column of the identity matrix, where k = p_j the j-th element of the permutation vector. The sign of the permutation is given by signum. It has the value (-1)^n, where n is the number of interchanges in the permutation.

The algorithm used in the decomposition is Gaussian Elimination with partial pivoting (Golub & Van Loan, Matrix Computations, Algorithm 3.4.1).

Idea for Future:
The "swap rows j and i_pivot" section could probably be made more efficient using a "matrix_row"-like object as done in GSL. (7/16/09 - I've tried this, and it doesn't seem to improve the speed significantly.)

Definition at line 86 of file lu_base.h.

◆ LU_det()

template<class mat_t >
double o2scl_linalg::LU_det ( const size_t  N,
const mat_t &  LU,
int  signum 
)

Compute the determinant of a matrix A from its LU decomposition, LU. The determinant is computed as the product of the diagonal elements of U and the sign of the row permutation signum.

Definition at line 348 of file lu_base.h.

◆ LU_invert()

template<class mat_t , class mat2_t , class mat_col_t >
int o2scl_linalg::LU_invert ( const size_t  N,
const mat_t &  LU,
const o2scl::permutation p,
mat2_t &  inverse 
)

These functions compute the inverse of a matrix A from its LU decomposition (LU,p), storing the result in the matrix inverse. The inverse is computed by solving the system A x = b for each column of the identity matrix. It is preferable to avoid direct use of the inverse whenever possible, as the linear solver functions can obtain the same result more efficiently and reliably.

Idea for Future:
Could rewrite to avoid mat_col_t, (9/16/09 - However, the function may be faster if mat_col_t is left in, so it's unclear what's best.)

Definition at line 309 of file lu_base.h.

◆ LU_lndet()

template<class mat_t >
double o2scl_linalg::LU_lndet ( const size_t  N,
const mat_t &  LU 
)

Compute the logarithm of the absolute value of the determinant of a matrix A, $ \ln|\det(A)| $, from its LU decomposition, LU. This function may be useful if the direct computation of the determinant would overflow or underflow.

Definition at line 370 of file lu_base.h.

◆ LU_refine()

template<class mat_t , class mat2_t , class vec_t , class vec2_t , class vec3_t >
int o2scl_linalg::LU_refine ( const size_t  N,
const mat_t &  A,
const mat2_t &  LU,
const o2scl::permutation p,
const vec_t &  b,
vec2_t &  x,
vec3_t &  residual 
)

These functions apply an iterative improvement to x, the solution of A x = b, using the LU decomposition of A into (LU,p). The initial residual r = A x - b is also computed and stored in residual.

Definition at line 271 of file lu_base.h.

◆ LU_sgndet()

template<class mat_t >
int o2scl_linalg::LU_sgndet ( const size_t  N,
const mat_t &  LU,
int  signum 
)

Compute the sign or phase factor of the determinant of a matrix A, $ \det(A)/|\det(A)| $, from its LU decomposition, LU.

Definition at line 389 of file lu_base.h.

◆ LU_solve()

template<class mat_t , class vec_t , class vec2_t >
int o2scl_linalg::LU_solve ( const size_t  N,
const mat_t &  LU,
const o2scl::permutation p,
const vec_t &  b,
vec2_t &  x 
)

This function solve the square system A x = b using the LU decomposition of A into (LU, p) given by gsl_linalg_LU_decomp or gsl_linalg_complex_LU_decomp.

Definition at line 248 of file lu_base.h.

◆ LU_svx()

template<class mat_t , class vec_t >
int o2scl_linalg::LU_svx ( const size_t  N,
const mat_t &  LU,
const o2scl::permutation p,
vec_t &  x 
)

These functions solve the square system A x = b in-place using the LU decomposition of A into (LU,p). On input x should contain the right-hand side b, which is replaced by the solution on output.

Definition at line 148 of file lu_base.h.

◆ QR_decomp_unpack()

template<class mat_t , class mat2_t , class mat3_t >
void o2scl_linalg::QR_decomp_unpack ( const size_t  M,
const size_t  N,
mat_t &  A,
mat2_t &  Q,
mat3_t &  R 
)

If O2scl is compiled with Armadillo support, this is specialized for arma::mat to use arma::qr_econ. If O2scl is compiled with Eigen support, this is specialized for Eigen::MatrixXd.

Definition at line 231 of file qr_base.h.

◆ QR_update()

template<class mat1_t , class mat2_t , class vec1_t , class vec2_t >
void o2scl_linalg::QR_update ( size_t  M,
size_t  N,
mat1_t &  Q,
mat2_t &  R,
vec1_t &  w,
vec2_t &  v 
)

The parameters M and N are the number of rows and columns of the matrix R.

* Q' R' = QR + u v^T
*       = Q (R + Q^T u v^T)
*       = Q (R + w v^T)
*
* where w = Q^T u.
*
* Algorithm from Golub and Van Loan, "Matrix Computations", Section
* 12.5 (Updating Matrix Factorizations, Rank-One Changes)

Definition at line 165 of file qr_base.h.

◆ qrstep()

template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void o2scl_linalg::qrstep ( size_t  M,
size_t  N,
size_t  n,
vec_t &  d,
vec2_t &  f,
mat_t &  U,
mat2_t &  V 
)

The vector d should be of size n, the vector f should be of size n, the matrix U should be of size (M,N), and the matrix V should be of size (N,N).

Definition at line 603 of file svdstep_base.h.

◆ qrstep_sub()

template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void o2scl_linalg::qrstep_sub ( size_t  M,
size_t  N,
size_t  n,
vec_t &  d,
vec2_t &  f,
mat_t &  U,
mat2_t &  V,
size_t  a 
)

The vector d should be of size n, the vector f should be of size n, the matrix U should be of size (M,n), and the matrix V should be of size (N,n).

A version of qrstep(), but starting at the a'th column of U, the a'th column of V, and the a'th entries of d and f.

This function is used by SV_decomp().

Definition at line 762 of file svdstep_base.h.

◆ solve_cyc_tridiag_nonsym()

template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t , class mem_t , class mem_vec_t >
void o2scl_linalg::solve_cyc_tridiag_nonsym ( const vec_t &  diag,
const vec2_t &  abovediag,
const vec3_t &  belowdiag,
const vec4_t &  rhs,
vec5_t &  x,
size_t  N,
mem_t &  m 
)

This function solves the system $ A x = b $ where $ A $ is a matrix of the form

*
*        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
*   belowdiag[0]       diag[1]  abovediag[1]   .....
*              0  belowdiag[1]       diag[2]
*              0             0  belowdiag[2]   .....
*            ...           ...
* abovediag[N-1]           ...

This function solves the following system without the corner elements and then use Sherman-Morrison formula to compensate for them.

Idea for Future:
Offer an option to avoid throwing on divide by zero?

Definition at line 338 of file tridiag_base.h.

◆ solve_cyc_tridiag_sym()

template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class mem_t , class mem_vec_t >
void o2scl_linalg::solve_cyc_tridiag_sym ( const vec_t &  diag,
const vec2_t &  offdiag,
const vec3_t &  b,
vec4_t &  x,
size_t  N,
mem_t &  m 
)

This function solves the system $ A x = b $ where $ A $ is a matrix of the form

*
*      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
*   offdiag[0]     diag[1]    offdiag[1]   .....
*            0  offdiag[1]       diag[2]
*            0           0    offdiag[2]   .....
*          ...         ...
* offdiag[N-1]         ...

See EngelnMullges96 .

Definition at line 241 of file tridiag_base.h.

◆ solve_tridiag_nonsym()

template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class vec5_t , class mem_t , class mem_vec_t >
void o2scl_linalg::solve_tridiag_nonsym ( const vec_t &  diag,
const vec2_t &  abovediag,
const vec3_t &  belowdiag,
const vec4_t &  rhs,
vec5_t &  x,
size_t  N,
mem_t &  m 
)

This function solves the system $ A x = b $ where $ A $ is a matrix of the form

* 
*       diag[0]  abovediag[0]             0   .....
*  belowdiag[0]       diag[1]  abovediag[1]   .....
*             0  belowdiag[1]       diag[2]
*             0             0  belowdiag[2]   .....

This function uses plain Gauss elimination, not bothering with the zeroes.

Idea for Future:
Offer an option to avoid throwing on divide by zero?

Definition at line 183 of file tridiag_base.h.

◆ solve_tridiag_sym()

template<class vec_t , class vec2_t , class vec3_t , class vec4_t , class mem_t , class mem_vec_t >
void o2scl_linalg::solve_tridiag_sym ( const vec_t &  diag,
const vec2_t &  offdiag,
const vec3_t &  b,
vec4_t &  x,
size_t  N,
mem_t &  m 
)

This function solves the system $ A x = b $ where $ A $ is a matrix of the form

*
*     diag[0]  offdiag[0]             0   .....
*  offdiag[0]     diag[1]    offdiag[1]   .....
*           0  offdiag[1]       diag[2]
*           0           0    offdiag[2]   .....

given the N diagonal elements in diag, N-1 diagonal elements in offdiag, and the N elements b from the RHS.

See EngelnMullges96 .

Definition at line 118 of file tridiag_base.h.

◆ SV_decomp()

template<class mat_t , class mat2_t , class vec_t , class vec2_t >
void o2scl_linalg::SV_decomp ( size_t  M,
size_t  N,
mat_t &  A,
mat2_t &  V,
vec_t &  S,
vec2_t &  work 
)

This factors matrix A of size (M,N) into

\[ A = U~D~V^T \]

where U is a column-orthogonal matrix of size (M,N) (stored in A), D is a diagonal matrix of size (N,N) (stored in the vector S of size N), and V is a orthogonal matrix of size (N,N). The vector work is a workspace vector of size N. The matrices U and V are constructed so that

\[ U^T~U = I \qquad \mathrm{and} \qquad V^T~V = V~V^T = I \]

This algorithm requres $ M \geq N $.

Todo:
Test N=1 case, N=2 case, and non-square matrices.

Definition at line 73 of file svd_base.h.

◆ SV_decomp_jacobi()

template<class mat_t , class mat2_t , class vec_t >
void o2scl_linalg::SV_decomp_jacobi ( size_t  M,
size_t  N,
mat_t &  A,
mat2_t &  Q,
vec_t &  S 
)

This factors matrix A of size (M,N) into

\[ A = U~D~V^T \]

where U is a column-orthogonal matrix of size (M,N) (stored in A), D is a diagonal matrix of size (N,N) (stored in the vector S of size N), and V is a orthogonal matrix of size (N,N).

This function computes singular values to higher relative accuracy than SV_decomp() and SV_decomp_mod().

Idea for Future:
There were originally a few GSL_COERCE_DBL calls which have been temporarily removed and could be restored.

Definition at line 438 of file svd_base.h.

◆ SV_decomp_mod()

template<class mat_t , class mat2_t , class mat3_t , class vec_t , class vec2_t >
void o2scl_linalg::SV_decomp_mod ( size_t  M,
size_t  N,
mat_t &  A,
mat2_t &  X,
mat3_t &  V,
vec_t &  S,
vec2_t &  work 
)

This factors matrix A of size (M,N) into

\[ A = U~D~V^T \]

where U is a column-orthogonal matrix of size (M,N) (stored in A), D is a diagonal matrix of size (N,N) (stored in the vector S of size N), and V is a orthogonal matrix of size (N,N). The vector work is a workspace vector of size N and the matrix X is a workspace of size (N,N).

Definition at line 290 of file svd_base.h.

◆ SV_solve()

template<class mat_t , class mat2_t , class vec_t , class vec2_t , class vec3_t >
void o2scl_linalg::SV_solve ( size_t  M,
size_t  N,
mat_t &  U,
mat2_t &  V,
vec_t &  S,
vec2_t &  b,
vec3_t &  x 
)

Solves a linear system using the output of SV_decomp(). Only non-zero singular values are used in computing solution. In the over-determined case, $ M>N $, the system is solved in the least-squares sense.

Definition at line 373 of file svd_base.h.

◆ svd2()

template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void o2scl_linalg::svd2 ( size_t  M,
size_t  N,
vec_t &  d,
vec2_t &  f,
mat_t &  U,
mat2_t &  V 
)

The parameter M is the number of rows in U and N is the number of rows in V. Both U and V assumed to have two columns.

Used in qrstep().

Definition at line 176 of file svdstep_base.h.

◆ svd2_sub()

template<class vec_t , class vec2_t , class mat_t , class mat2_t >
void o2scl_linalg::svd2_sub ( size_t  M,
size_t  N,
vec_t &  d,
vec2_t &  f,
mat_t &  U,
mat2_t &  V,
size_t  a 
)

The parameter M is the number of rows in U and N is the number of rows in V. Both U and V assumed to have two columns.

Used in qrstep_sub().

Definition at line 316 of file svdstep_base.h.

◆ trailing_eigenvalue()

template<class vec_t , class vec2_t >
double o2scl_linalg::trailing_eigenvalue ( size_t  n,
const vec_t &  d,
const vec2_t &  f 
)

The parameter n is the size of the vector d.

Used in qrstep() and qrstep_sub().

Definition at line 85 of file svdstep_base.h.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).