`stdlib_specialmatrices`#

Introduction #

The stdlib_specialmatrices module provides derived types and specialized drivers for highly structured matrices often encountered in scientific computing as well as control and signal processing applications. These include:

Tridiagonal matrices
Symmetric tridiagonal matrices
Circulant matrices (not yet supported)
Toeplitz matrices (not yet supported)
Hankel matrices (not yet supported)

In addition, it also provides a Poisson2D matrix type (not yet supported) corresponding to the sparse block tridiagonal matrix obtained from discretizing the Laplace operator on a 2D grid with the standard second-order accurate central finite-difference scheme.

List of derived types for special matrices #

Below is a list of the currently supported derived types corresponding to different special matrices. Note that this module is under active development and this list will eventually grow.

Tridiagonal matrices {#Tridiagonal}#

Status#

Experimental

Description#

Tridiagonal matrices are ubiquitous in scientific computing and often appear when discretizing 1D differential operators. A generic tridiagonal matrix has the following structure:

\[\begin{split} A = \begin{bmatrix} a_1 & b_1 \\ c_1 & a_2 & b_2 \\ & \ddots & \ddots & \ddots \\ & & c_{n-2} & a_{n-1} & b_{n-1} \\ & & & c_{n-1} & a_n \end{bmatrix}. \end{split}\]

Hence, only one vector of size n and two of size n-1 need to be stored to fully represent the matrix. This particular structure also lends itself to specialized implementations for many linear algebra tasks. Interfaces to the most common ones will soon be provided by stdlib_specialmatrices. Tridiagonal matrices are available with all supported data types as tridiagonal_<kind>_type, for example:

tridiagonal_sp_type : Tridiagonal matrix of size n with real/single precision data.
tridiagonal_dp_type : Tridiagonal matrix of size n with real/double precision data.
tridiagonal_xdp_type : Tridiagonal matrix of size n with real/extended precision data.
tridiagonal_qp_type : Tridiagonal matrix of size n with real/quadruple precision data.
tridiagonal_csp_type : Tridiagonal matrix of size n with complex/single precision data.
tridiagonal_cdp_type : Tridiagonal matrix of size n with complex/double precision data.
tridiagonal_cxdp_type : Tridiagonal matrix of size n with complex/extended precision data.
tridiagonal_cqp_type : Tridiagonal matrix of size n with complex/quadruple precision data.

Syntax#

To construct a tridiagonal matrix from already allocated arrays dl (lower diagonal, size n-1), dv (main diagonal, size n) and du (upper diagonal, size n-1):

A = tridiagonal (dl, dv, du [, err])

To construct a tridiagonal matrix of size n x n with scalar diagonal elements dl (lower diagonal), dv (main diagonal), and du (upper diagonal):

A = tridiagonal (dl, dv, du, n [, err])

Where err is an optional intent(out) argument of type linalg_state_type. When err is not present, the corresponding tridiagonal call is pure and signals invalid sizes or array shape mismatches via an internal error stop. When err is present, the impure overload is used instead, and any such error is reported by setting the err variable rather than executing an error stop.

Arguments#

Constructing from arrays#

A = tridiagonal(dl, dv, du [, err])

dl: Shall be a rank-1 array of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
dv: Shall be a rank-1 array of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
du: Shall be a rank-1 array of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) variable. It is an intent(out) argument. If the input arrays are not of compatible sizes, the routine calls linalg_error_handling. If err is present, it is assigned LINALG_VALUE_ERROR. If err is not present, the program terminates via error stop. In either case, the elements of the resulting matrix A are unassigned.

Constructing from constant diagonal values#

A = tridiagonal(dl, dv, du, n [, err])

dl: Shall be a scalar of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
dv: Shall be a scalar of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
du: Shall be a scalar of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
n: Shall be a positive integer. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) variable. It is an intent(out) argument. If a non-positive matrix size is provided, the routine calls linalg_error_handling. If err is present, it is assigned LINALG_VALUE_ERROR. If err is not present, the program terminates via error stop. In either case, the elements of the resulting matrix A are unassigned.

Example#

program example_tridiagonal_matrix
   use stdlib_linalg_constants, only: dp
   use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal
   implicit none

   integer, parameter :: n = 5
   type(Tridiagonal_dp_type) :: A
   real(dp) :: dl(n - 1), dv(n), du(n - 1)

   ! Generate random tridiagonal elements.
   call random_number(dl)
   call random_number(dv)
   call random_number(du)

   ! Create the corresponding Tridiagonal matrix.
   A = Tridiagonal(dl, dv, du)

end program example_tridiagonal_matrix

Symmetric tridiagonal matrices {#Sym_tridiagonal}#

Status#

Experimental

Description#

Symmetric tridiagonal matrices are a special case of tridiagonal matrices in which the subdiagonal and superdiagonal are identical. A generic symmetric tridiagonal matrix has the following structure:

\[\begin{split} A = \begin{bmatrix} a_1 & b_1 \\ b_1 & a_2 & b_2 \\ & \ddots & \ddots & \ddots \\ & & b_{n-2} & a_{n-1} & b_{n-1} \\ & & & b_{n-1} & a_n \end{bmatrix}. \end{split}\]

Hence, only one vector of size n and one vector of size n-1 need to be stored to fully represent the matrix. This particular structure also lends itself to specialized implementations for many linear algebra tasks. Interfaces to the most common ones will soon be provided by stdlib_specialmatrices. Symmetric tridiagonal matrices are available with all supported data types as sym_tridiagonal_<kind>_type, for example:

sym_tridiagonal_sp_type : Symmetric tridiagonal matrix of size n with real/single precision data.
sym_tridiagonal_dp_type : Symmetric tridiagonal matrix of size n with real/double precision data.
sym_tridiagonal_xdp_type : Symmetric tridiagonal matrix of size n with real/extended precision data.
sym_tridiagonal_qp_type : Symmetric tridiagonal matrix of size n with real/quadruple precision data.
sym_tridiagonal_csp_type : Symmetric tridiagonal matrix of size n with complex/single precision data.
sym_tridiagonal_cdp_type : Symmetric tridiagonal matrix of size n with complex/double precision data.
sym_tridiagonal_cxdp_type : Symmetric tridiagonal matrix of size n with complex/extended precision data.
sym_tridiagonal_cqp_type : Symmetric tridiagonal matrix of size n with complex/quadruple precision data.

Syntax#

To construct a symmetric tridiagonal matrix from already allocated arrays du (off-diagonal, size n-1) and dv (main diagonal, size n):

A = sym_tridiagonal (du, dv [, err])

To construct a symmetric tridiagonal matrix of size n x n with scalar diagonal elements du (off-diagonal) and dv (main diagonal):

A = sym_tridiagonal (du, dv, n [, err])

Arguments#

Constructing from arrays#

A = sym_tridiagonal(du, dv [, err])

du: Shall be a rank-1 array of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
dv: Shall be a rank-1 array of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) variable. It is an intent(out) argument. If the input arrays are not of compatible sizes, the routine calls linalg_error_handling. If err is present, it is assigned LINALG_VALUE_ERROR. If err is not present, the program terminates via error stop. In either case, the elements of the resulting matrix A are unassigned.

Constructing from constant diagonal values#

A = sym_tridiagonal(du, dv, n [, err])

du: Shall be a scalar of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
dv: Shall be a scalar of type real or complex with the same kind as the resulting matrix. It is an intent(in) argument.
n: Shall be a positive integer. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) variable. It is an intent(out) argument. If a non-positive matrix size is provided, the routine calls linalg_error_handling. If err is present, it is assigned LINALG_VALUE_ERROR. If err is not present, the program terminates via error stop. In either case, the elements of the resulting matrix A are unassigned.

Example#

program example_sym_tridiagonal_matrix
   use stdlib_linalg_constants, only: dp
   use stdlib_specialmatrices, only: sym_tridiagonal_dp_type, sym_tridiagonal
   implicit none

   integer, parameter :: n = 5
   type(sym_tridiagonal_dp_type) :: A
   real(dp) :: du(n - 1), dv(n)

   ! Generate random symmteric tridiagonal elements.
   call random_number(du)
   call random_number(dv)

   ! Create the corresponding symmetric tridiagonal matrix.
   A = sym_tridiagonal(du, dv)

end program example_sym_tridiagonal_matrix

Specialized drivers for linear algebra tasks #

Below is a list of all the specialized drivers for linear algebra tasks currently provided by the stdlib_specialmatrices module.

Matrix-vector products with `spmv` {#spmv}#

Status#

Experimental

Description#

With the exception of extended precision and quadruple precision, all the types provided by stdlib_specialmatrices benefit from specialized kernels for matrix-vector products accessible via the common spmv interface.

For tridiagonal and symmetric tridiagonal matrices, the backend is either LAPACK lagtm or the generalized routine glagtm, depending on the values and types of alpha and beta.

Syntax#

call spmv (A, x, y [, alpha, beta, op])

Arguments#

A : Shall be a matrix of one of the types provided by stdlib_specialmatrices. It is an intent(in) argument.
x : Shall be a rank-1 or rank-2 array with the same kind as A. It is an intent(in) argument.
y : Shall be a rank-1 or rank-2 array with the same kind as A. It is an intent(inout) argument.
alpha (optional) : Scalar value of the same type as x. It is an intent(in) argument. By default, alpha = 1.
beta (optional) : Scalar value of the same type as y. It is an intent(in) argument. By default beta = 0.
op (optional) : In-place operator identifier. Shall be a character(1) argument. It can have any of the following values: N: no transpose, T: transpose, H: hermitian or complex transpose.

Examples#

program example_tridiagonal_matrix
   use stdlib_linalg_constants, only: dp
   use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal, dense, spmv
   implicit none

   integer, parameter :: n = 5
   type(tridiagonal_dp_type) :: A
   real(dp) :: dl(n - 1), dv(n), du(n - 1)
   real(dp) :: x(n), y(n), y_dense(n)
   integer :: i

   ! Create an arbitrary tridiagonal matrix.
   dl = [(i, i=1, n - 1)]; dv = [(2*i, i=1, n)]; du = [(3*i, i=1, n - 1)]
   A = tridiagonal(dl, dv, du)

   ! Initialize vectors.
   x = 1.0_dp; y = 0.0_dp; y_dense = 0.0_dp

   ! Perform matrix-vector products.
   call spmv(A, x, y)
   y_dense = matmul(dense(A), x)

   print *, 'dense       :', y_dense
   print *, 'Tridiagonal :', y

end program example_tridiagonal_matrix

program example_tridiagonal_matrix_cdp
    use stdlib_linalg_constants, only: dp
    use stdlib_specialmatrices, only: tridiagonal_cdp_type, tridiagonal, dense, spmv
    implicit none

    integer, parameter :: n = 5
    type(tridiagonal_cdp_type) :: A
    complex(dp) :: dl(n-1), dv(n), du(n-1)
    complex(dp) :: x(n), y(n), y_dense(n)
    integer :: i
    complex(dp) :: alpha, beta

    dl = [(cmplx(i,i, dp), i=1, n - 1)]
    dv = [(cmplx(2*i,2*i, dp), i=1, n)]
    du = [(cmplx(3*i,3*i, dp), i=1, n - 1)]

    A = tridiagonal(dl, dv, du)

    x = (1.0_dp, 0.0_dp)
    y = (3.0_dp, -7.0_dp)
    y_dense = (0.0_dp, 0.0_dp)
    alpha = cmplx(2.0_dp, 3.0_dp)
    beta = cmplx(-1.0_dp, 5.0_dp)

    y_dense = alpha * matmul(dense(A), x) + beta * y
    call spmv(A, x, y, alpha, beta)

    print *, 'dense       :', y_dense
    print *, 'Tridiagonal :', y
end program example_tridiagonal_matrix_cdp

Utility functions #

`dense` : converting a special matrix to a standard Fortran array {#dense}#

Status#

Experimental

Description#

Utility function to convert all the matrix types provided by stdlib_specialmatrices to a standard rank-2 array of the appropriate kind.

Syntax#

B = dense (A)

Arguments#

A : Shall be a matrix of one of the types provided by stdlib_specialmatrices. It is an intent(in) argument.
B : Shall be a rank-2 allocatable array of the appropriate real or complex kind.

`transpose` : Transposition of a special matrix {#transpose}#

Status#

Experimental

Description#

Utility function returning the transpose of a special matrix. The returned matrix is of the same type and kind as the input one.

Syntax#

B = transpose (A)

Arguments#

A : Shall be a matrix of one of the types provided by stdlib_specialmatrices. It is an intent(in) argument.
B : Shall be a matrix of one of the same type and kind as A.

`hermitian` : Complex-conjugate transpose of a special matrix {#hermitian}#

Status#

Experimental

Description#

Utility function returning the complex-conjugate transpose of a special matrix. The returned matrix is of the same type and kind as the input one. For real-valued matrices, hermitian is equivalent to transpose.

Syntax#

B = hermitian (A)

Arguments#

A : Shall be a matrix of one of the types provided by stdlib_specialmatrices. It is an intent(in) argument.
B : Shall be a matrix of one of the same type and kind as A.

Operator overloading (`+`, `-`, `*`) {#operators}#

Status#

Experimental

Description#

The definition of all standard arithmetic operators have been overloaded to be applicable for the matrix types defined by stdlib_specialmatrices:

Overloading the + operator for adding two matrices of the same type and kind.
Overloading the - operator for subtracting two matrices of the same type and kind.
Overloading the * for scalar-matrix multiplication.

Syntax#

Adding two matrices of the same type:

C = A operator(+) B

Subtracting two matrices of the same type:

C = A operator(-) B

Scalar multiplication

B = alpha operator(*) A

Note

For addition (+) and subtraction (-), the operand matrices A and B must be of the same type and kind and must have matching dimensions (n). The result C of the expression will then have the same dimensions as A (and B) when the operation is valid. For scalar multiplication (*), A and B need to be of the same type and kind, while alpha is either real or complex (with the same kind again) depending on the type being used.

stdlib_specialmatrices#

Status#

Description#

Syntax#

Arguments#

Constructing from arrays#

Constructing from constant diagonal values#

Example#

Status#

Description#

Syntax#

Arguments#

Constructing from arrays#

Constructing from constant diagonal values#

Example#

Status#

Description#

Syntax#

Arguments#

Examples#

Status#

Description#

Syntax#

Arguments#

Status#

Description#

Syntax#

Arguments#

Status#

Description#

Syntax#

Arguments#

Status#

Description#

Syntax#

`stdlib_specialmatrices`#