LU

Driver Routines

template<Number T, DataLocation L, template<typename, typename> typename S1, template<typename, typename> typename S2>
internal::LUSolverDriver<T, L, Matrix<T, L, S2>> lusolve(const Matrix<T, L, S1> &A, const Matrix<T, L, S2> &B)

Solves the linear system \( \mathbf{A} \; \mathbf{X} = \mathbf{B}\) using LU factorization. ( \( \mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}\)). This operation can be done in-place by calling X = lusolve(A, X).

Parameters:
  • A[in] coefficient matrix

  • B[in] right-hand side of the system

Returns:

An abstract object representing the solution of the linear system.

template<Number T, DataLocation L, template<typename, typename> typename S1, template<typename, typename> typename S2>
internal::LUSolverDriver<T, L, Matrix<T, L, S2>, kInPlace> lusolve(const Matrix<T, L, S1> &A, const Matrix<T, L, S2> &B, InPlace)

Solves the linear system \( \mathbf{A} \; \mathbf{X} = \mathbf{B}\) using LU factorization ( \( \mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}\)). This operation can be done in-place by calling X = lusolve(A, X).

Warning

The matrix A will be overwritten by the factor L and U of the factorization.

Parameters:
  • A[in] coefficient matrix

  • B[in] right-hand side of the system

Returns:

An abstract object representing the solution of the linear system.

template<Number T, DataLocation L, template<typename, typename> typename S1, template<typename, typename> typename S2>
internal::LUSolverDriver<T, L, Vector<T, L, S2>> lusolve(const Matrix<T, L, S1> &A, const Vector<T, L, S2> &b)

Solves the linear system \( \mathbf{A} \; \vec{x} = \vec{b}\) using LU factorization ( \( \mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}\)). This operation can be done in-place by calling x = lusolve(A, x).

Parameters:
  • A[in] coefficient matrix

  • b[in] right-hand side of the system

Returns:

An abstract object representing the solution of the linear system.

template<Number T, DataLocation L, template<typename, typename> typename S1, template<typename, typename> typename S2>
internal::LUSolverDriver<T, L, Vector<T, L, S2>> lusolve(const Matrix<T, L, S1> &A, const Vector<T, L, S2> &b, InPlace)

Solves the linear system \( \mathbf{A} \; \vec{x} = \vec{b}\) using LU factorization ( \( \mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}\)). This operation can be done in-place by calling x = lusolve(A, x).

Warning

The matrix A will be overwritten by the factor L and U of the factorization.

Parameters:
  • A[in] coefficient matrix

  • b[in] right-hand side of the system

Returns:

An abstract object representing the solution of the linear system.

template<Number T, DataLocation L, template<typename, typename> typename S>
LU<T, L> lu(const Matrix<T, L, S> &A)

Factorize the square matrix A as \( \mathbf{A} = \mathbf{PLU}\), where U is an upper triangular matrix, L is a lower triangular matrix with a unit diagonal and P is a permutation matrix.

Parameters:

A[in] input matrix

template<Number T, DataLocation L, template<typename, typename> typename S>
LU<T, L, kInPlace> lu(const Matrix<T, L, S> &A, InPlace)

Factorize the square matrix A as \( \mathbf{A} = \mathbf{PLU}\), where U is an upper triangular matrix, L is a lower triangular matrix with a unit diagonal and P is a permutation matrix.

Warning

The matrix A will be overwritten by the factor L and U of the factorization.

Parameters:

A[in] input matrix

Factorization

template<Number T, DataLocation Location, int Flags>
class LU

Object representing the LU factorization of a square matrix A: \( \mathbf{A} = \mathbf{PLU}\), where U is an upper triangular matrix, L is a lower triangular matrix with a unit diagonal and P is a permutation matrix.

Usage:

Matrix<double> A(n, n);
Matrix<double> invA(n, n);
Matrix<double> B(n, nrhs);
Matrix<double> X(n, nrhs);
LU<double> F = lu(A);           // Factorize AP = LU
X = F.solve(B);                 // Solve AX = B
double cond = F.cond();       // Calculates the reciprocal condition number using L1 norm
invA = F.inv();                 // Calculate inv(A)

Template Parameters:
  • T – Floating-point type

  • Location – Location of the data.

  • Flags – Flags for the QR decomposition. Available: kDefaultFlags (default flags) and kInPlace (the factorization is done in-place, overwriting A).

Public Functions

inline LU()

Default constructor.

template<template<typename, typename> typename S>
inline LU(const Matrix<T, Location, S> &A)

Creates a new solver and factorize the matrix A.

Parameters:

A[in] input matrix

virtual ~LU() = default

Destructor.

template<template<typename, typename> typename S>
inline LU &factorize(const Matrix<T, Location, S> &A)

Calculates the LU factorization of the matrix A.

Parameters:

A[in] input matrix

inline internal::scalar_type_t<T> rcond()
Returns:

The reciprocal of the condition number \( \kappa_1(\mathbf{A}) = \|\mathbf{A}\|_1 \|\mathbf{A}^{-1}\|_1\).

template<template<typename, typename> typename S>
inline internal::LUSolve<T, Location, Matrix<T, Location, S>, Flags> solve(Matrix<T, Location, S> &B, char trans = 'N')

Solves the linear system \( \text{op}(\mathbf{A}) \; \mathbf{X} = \mathbf{B}\). This operation can be done in-place by calling X = lu_fact.solve(X).

Parameters:
  • B[in] right hand side

  • trans[in] transpose operation: N (no transpose, \( \text{op}(\mathbf{A}) = \mathbf{A} \)), T (transpose, \( \text{op}(\mathbf{A}) = \mathbf{A}^T \)) or C (hermitian, complex only, \( \text{op}(\mathbf{A}) = \mathbf{A}^H \))

template<template<typename, typename> typename S>
inline internal::LUSolve<T, Location, Vector<T, Location, S>, Flags> solve(Vector<T, Location, S> &b, char trans = 'N')

Solves the linear system \( \text{op}(\mathbf{A}) \; \vec{x} = \vec{b}\). This operation can be done in-place by calling x = lu_fact.solve(x).

Parameters:
  • b[in] right hand side

  • trans[in] transpose operation: N (no transpose, \( \text{op}(\mathbf{A}) = \mathbf{A} \)), T (transpose, \( \text{op}(\mathbf{A}) = \mathbf{A}^T \)) or C (hermitian, complex only, \( \text{op}(\mathbf{A}) = \mathbf{A}^H \))

inline internal::LUInverse<T, Location, Flags> inv()

Calculate the inverse of A.

Warning

To solve the linear system \( \mathbf{A} \; \mathbf{X} = \mathbf{B}\), use call solve() instead of explicitly forming the matrix inverse \( \mathbf{A}^{-1}\) and then multiply \( \mathbf{X} = \mathbf{A}^{-1} \mathbf{B}\). This is more efficient and more accurate.

inline const mat_type &compressed_mat() const
Returns:

The underlying matrix. The matrix U is stored at and above the diagonal, while L is stored below the diagonal. The unit diagonal of L is not stored.

inline TriangularMap<T, kLowerUnit, Location> L() const
Returns:

The lower triangular matrix L from the LU factorization. The unit diagonal of L is not stored.

inline TriangularMap<T, kUpper, Location> U() const
Returns:

The upper triangular matrix U from the LU factorization.

inline const Vector<index_t, Location> &pivots() const
Returns:

The index of the pivots using during the factorization.

inline index_t rows() const
Returns:

The number of rows of the underlying matrix

inline index_t cols() const
Returns:

The number of cols of the underlying matrix