stdlib_specialmatrices moduleThe 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:
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.
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.
Experimental
Tridiagonal matrices are ubiquituous in scientific computing and often appear when discretizing 1D differential operators.
A generic tridiagonal matrix has the following structure:
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.dl (lower diagonal, size n-1), dv (main diagonal, size n) and du (upper diagonal, size n-1):A = tridiagonal (dl, dv, du)
n x n with constant diagonal elements dl, dv, and du:A = tridiagonal (dl, dv, du, n)
program example_tridiagonal_matrix
use stdlib_linalg_constants, only: dp
use stdlib_specialmatrices
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
Below is a list of all the specialized drivers for linear algebra tasks currently provided by the stdlib_specialmatrices module.
spmvExperimental
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.
tridiagonal matrices, the LAPACK lagtm backend is being used.call spmv (A, x, y [, alpha, beta, op])
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.
Warning
Due to limitations of the underlying lapack driver, currently alpha and beta can only take one of the values [-1, 0, 1] for tridiagonal and symtridiagonal matrices. See lagtm for more details.
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
dense : converting a special matrix to a standard Fortran arrayExperimental
Utility function to convert all the matrix types provided by stdlib_specialmatrices to a standard rank-2 array of the appropriate kind.
B = dense (A)
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 matrixExperimental
Utility function returning the transpose of a special matrix. The returned matrix is of the same type and kind as the input one.
B = transpose (A)
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 matrixExperimental
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.
B = hermitian (A)
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.
+, -, *)Experimental
The definition of all standard artihmetic operators have been overloaded to be applicable for the matrix types defined by stdlib_specialmatrices:
+ operator for adding two matrices of the same type and kind.- operator for subtracting two matrices of the same type and kind.* for scalar-matrix multiplication.C = A [[stdlib_specialmatrices(module):operator(+)(interface)]] B
C = A [[stdlib_specialmatrices(module):operator(-)(interface)]] B
B = alpha [[stdlib_specialmatrices(module):operator(*)(interface)]] A
Note
For addition (+) and subtraction (-), matrices A, B and C all need to be of the same type and kind. 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.