diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md new file mode 100644 index 000000000..d73866628 --- /dev/null +++ b/doc/specs/stdlib_specialmatrices.md @@ -0,0 +1,215 @@ +--- +title: specialmatrices +--- + +# The `stdlib_specialmatrices` module + +[TOC] + +## 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 (not yet supported) +- 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 ubiquituous in scientific computing and often appear when discretizing 1D differential operators. +A generic tridiagonal matrix has the following structure +$$ + 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}. +$$ +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`. +To date, `stdlib_specialmatrices` supports the following data types: + +- `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 = ` [[stdlib_specialmatrices(module):Tridiagonal(interface)]] `(dl, dv, du)` + +- To construct a tridiagonal matrix of size `n x n` with constant diagonal elements `dl`, `dv`, and `du`: + +`A = ` [[stdlib_specialmatrices(module):Tridiagonal(interface)]] `(dl, dv, du, n)` + +#### Example + +```fortran +{!example/specialmatrices/example_tridiagonal_dp_type.f90!} +``` + +## 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` matrices, the LAPACK `lagtm` backend is being used. + +#### Syntax + +`call ` [[stdlib_specialmatrices(module):spmv(interface)]] `(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. + +@warning +Due to some underlying `lapack`-related designs, `alpha` and `beta` can only take values in `[-1, 0, 1]` for `Tridiagonal` and `SymTridiagonal` matrices. See `lagtm` for more details. +@endwarning + +#### Examples + +```fortran +{!example/specialmatrices/example_specialmatrices_dp_spmv.f90!} +``` + +## 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 = ` [[stdlib_specialmatrices(module):dense(interface)]] `(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 = ` [[stdlib_specialmatrices(module):transpose(interface)]] `(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 = ` [[stdlib_specialmatrices(module):hermitian(interface)]] `(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 artihmetic 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` [[stdlib_specialmatrices(module):operator(+)(interface)]] `B` + +- Subtracting two matrices of the same type: + +`C = A` [[stdlib_specialmatrices(module):operator(-)(interface)]] `B` + +- Scalar multiplication + +`B = alpha` [[stdlib_specialmatrices(module):operator(*)(interface)]] `A` + +@note +For addition (`+`) and subtraction (`-`), the 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. +@endnote diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 3c61e8020..108b5b0d1 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -24,6 +24,7 @@ add_subdirectory(random) add_subdirectory(selection) add_subdirectory(sorting) add_subdirectory(specialfunctions_gamma) +add_subdirectory(specialmatrices) add_subdirectory(stats) add_subdirectory(stats_distribution_exponential) add_subdirectory(stats_distribution_normal) diff --git a/example/linalg/example_sparse_data_accessors.f90 b/example/linalg/example_sparse_data_accessors.f90 index e23164524..73ff440c2 100644 --- a/example/linalg/example_sparse_data_accessors.f90 +++ b/example/linalg/example_sparse_data_accessors.f90 @@ -1,49 +1,49 @@ program example_sparse_data_accessors - use stdlib_linalg_constants, only: dp - use stdlib_sparse - implicit none + use stdlib_linalg_constants, only: dp + use stdlib_sparse + implicit none - real(dp) :: mat(2,2) - real(dp), allocatable :: dense(:,:) - type(CSR_dp_type) :: CSR - type(COO_dp_type) :: COO - integer :: i, j, locdof(2) + real(dp) :: mat(2, 2) + real(dp), allocatable :: dense_matrix(:, :) + type(CSR_dp_type) :: CSR + type(COO_dp_type) :: COO + integer :: i, j, locdof(2) - ! Initial data - mat(:,1) = [1._dp,2._dp] - mat(:,2) = [2._dp,1._dp] - allocate(dense(5,5) , source = 0._dp) - do i = 0, 3 - dense(1+i:2+i,1+i:2+i) = dense(1+i:2+i,1+i:2+i) + mat - end do + ! Initial data + mat(:, 1) = [1._dp, 2._dp] + mat(:, 2) = [2._dp, 1._dp] + allocate (dense_matrix(5, 5), source=0._dp) + do i = 0, 3 + dense_matrix(1 + i:2 + i, 1 + i:2 + i) = dense_matrix(1 + i:2 + i, 1 + i:2 + i) + mat + end do - print *, 'Original Matrix' - do j = 1 , 5 - print '(5f8.1)',dense(j,:) - end do + print *, 'Original Matrix' + do j = 1, 5 + print '(5f8.1)', dense_matrix(j, :) + end do - ! Initialize CSR data and reset dense reference matrix - call dense2coo(dense,COO) - call coo2csr(COO,CSR) - CSR%data = 0._dp - dense = 0._dp + ! Initialize CSR data and reset dense reference matrix + call dense2coo(dense_matrix, COO) + call coo2csr(COO, CSR) + CSR%data = 0._dp + dense_matrix = 0._dp - ! Iteratively add blocks of data - do i = 0, 3 - locdof(1:2) = [1+i,2+i] - call CSR%add(locdof,locdof,mat) - ! lets print a dense view of every step - call csr2dense(CSR,dense) - print '(A,I2)', 'Add block :', i+1 - do j = 1 , 5 - print '(5f8.1)',dense(j,:) - end do - end do + ! Iteratively add blocks of data + do i = 0, 3 + locdof(1:2) = [1 + i, 2 + i] + call CSR%add(locdof, locdof, mat) + ! lets print a dense view of every step + call csr2dense(CSR, dense_matrix) + print '(A,I2)', 'Add block :', i + 1 + do j = 1, 5 + print '(5f8.1)', dense_matrix(j, :) + end do + end do - ! Request values from the matrix - print *, '' - print *, 'within sparse pattern :',CSR%at(2,1) - print *, 'outside sparse pattern :',CSR%at(5,2) - print *, 'outside matrix pattern :',CSR%at(7,7) - -end program example_sparse_data_accessors \ No newline at end of file + ! Request values from the matrix + print *, '' + print *, 'within sparse pattern :', CSR%at(2, 1) + print *, 'outside sparse pattern :', CSR%at(5, 2) + print *, 'outside matrix pattern :', CSR%at(7, 7) + +end program example_sparse_data_accessors diff --git a/example/specialmatrices/CMakeLists.txt b/example/specialmatrices/CMakeLists.txt new file mode 100644 index 000000000..1f691bbde --- /dev/null +++ b/example/specialmatrices/CMakeLists.txt @@ -0,0 +1,2 @@ +ADD_EXAMPLE(specialmatrices_dp_spmv) +ADD_EXAMPLE(tridiagonal_dp_type) diff --git a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 new file mode 100644 index 000000000..ebfcf34d0 --- /dev/null +++ b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 @@ -0,0 +1,26 @@ +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) + 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 diff --git a/example/specialmatrices/example_tridiagonal_dp_type.f90 b/example/specialmatrices/example_tridiagonal_dp_type.f90 new file mode 100644 index 000000000..7f221b5d6 --- /dev/null +++ b/example/specialmatrices/example_tridiagonal_dp_type.f90 @@ -0,0 +1,18 @@ +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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 933da34de..58e387104 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,14 +32,14 @@ set(fppFiles stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp stdlib_linalg_eigenvalues.fypp - stdlib_linalg_solve.fypp + stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_qr.fypp stdlib_linalg_inverse.fypp stdlib_linalg_pinv.fypp stdlib_linalg_norms.fypp stdlib_linalg_state.fypp - stdlib_linalg_svd.fypp + stdlib_linalg_svd.fypp stdlib_linalg_cholesky.fypp stdlib_linalg_schur.fypp stdlib_optval.fypp @@ -52,6 +52,8 @@ set(fppFiles stdlib_sparse_conversion.fypp stdlib_sparse_kinds.fypp stdlib_sparse_spmv.fypp + stdlib_specialmatrices.fypp + stdlib_specialmatrices_tridiagonal.fypp stdlib_specialfunctions_gamma.fypp stdlib_stats.fypp stdlib_stats_corr.fypp diff --git a/src/stdlib_sparse_kinds.fypp b/src/stdlib_sparse_kinds.fypp index ceba2a62d..bf03d2747 100644 --- a/src/stdlib_sparse_kinds.fypp +++ b/src/stdlib_sparse_kinds.fypp @@ -589,5 +589,5 @@ contains end subroutine #:endfor - -end module stdlib_sparse_kinds \ No newline at end of file + +end module stdlib_sparse_kinds diff --git a/src/stdlib_sparse_spmv.fypp b/src/stdlib_sparse_spmv.fypp index a92521099..76c78f9e1 100644 --- a/src/stdlib_sparse_spmv.fypp +++ b/src/stdlib_sparse_spmv.fypp @@ -589,5 +589,5 @@ contains end subroutine #:endfor - -end module \ No newline at end of file + +end module diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp new file mode 100644 index 000000000..d37b86a4a --- /dev/null +++ b/src/stdlib_specialmatrices.fypp @@ -0,0 +1,228 @@ +#:include "common.fypp" +#:set RANKS = range(1, 2+1) +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +module stdlib_specialmatrices + !! Provides derived-types and associated specialized linear algebra drivers + !! for highly-structured matrices commonly encountered in the discretization + !! of partial differential equations, as well as control and signal processing + !! applications. ([Specifications](../page/specs/stdlib_specialmatrices.html)) + use stdlib_linalg_constants + implicit none + private + public :: Tridiagonal + public :: spmv + public :: dense, transpose, hermitian + public :: operator(*), operator(+), operator(-) + + !-------------------------------------- + !----- ------ + !----- TYPE DEFINITIONS ------ + !----- ------ + !-------------------------------------- + + !--> Tridiagonal matrices + #:for k1, t1, s1 in (KINDS_TYPES) + type, public :: Tridiagonal_${s1}$_type + !! Base type to define a `Tridiagonal` matrix. + private + ${t1}$, allocatable :: dl(:), dv(:), du(:) + integer(ilp) :: n + end type + #:endfor + + !-------------------------------- + !----- ----- + !----- CONSTRUCTORS ----- + !----- ----- + !-------------------------------- + + interface Tridiagonal + !! ([Specifications](../page/specs/stdlib_specialmatrices.html#Tridiagonal)) This + !! interface provides different methods to construct a `Tridiagonal` matrix. Only + !! the non-zero elements of \( A \) are stored, i.e. + !! + !! \[ + !! 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}. + !! \] + !! + !! #### Syntax + !! + !! - Construct a real `Tridiagonal` matrix from rank-1 arrays: + !! + !! ```fortran + !! integer, parameter :: n + !! real(dp), allocatable :: dl(:), dv(:), du(:) + !! type(Tridiagonal_rdp_type) :: A + !! integer :: i + !! + !! dl = [(i, i=1, n-1)]; dv = [(2*i, i=1, n)]; du = [(3*i, i=1, n)] + !! A = Tridiagonal(dl, dv, du) + !! ``` + !! + !! - Construct a real `Tridiagonal` matrix with constant diagonals: + !! + !! ```fortran + !! integer, parameter :: n + !! real(dp), parameter :: a = 1.0_dp, b = 1.0_dp, c = 2.0_dp + !! type(Tridiagonal_rdp_type) :: A + !! + !! A = Tridiagonal(a, b, c, n) + !! ``` + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A) + !! Construct a `Tridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dl(:), dv(:), du(:) + !! Tridiagonal matrix elements. + type(Tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + end function + + pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A) + !! Construct a `Tridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dl, dv, du + !! Tridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(Tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + end function + #:endfor + end interface + + !---------------------------------- + !----- ----- + !----- LINEAR ALGEBRA ----- + !----- ----- + !---------------------------------- + + interface spmv + !! ([Specifications](../page/specs/stdlib_specialmatrices.html#spmv)) This + !! interface provides methods to compute the matrix-vector product + !! + !! $$ y = \alpha \mathrm{op}(A) x + \beta y$$ + !! + !! for the different matrix types defined by `stdlib_specialmatrices`. + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op) + type(Tridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$ + real(${k1}$), intent(in), optional :: alpha + real(${k1}$), intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + #:endfor + end interface + + !------------------------------------- + !----- ----- + !----- UTILITY FUNCTIONS ----- + !----- ----- + !------------------------------------- + + interface dense + !! This interface provides methods to convert a matrix of one of the + !! types defined by `stdlib_specialmatrices` to a standard rank-2 array. + !! ([Specifications](../page/specs/stdlib_specialmatrices.html#dense)) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function tridiagonal_to_dense_${s1}$(A) result(B) + !! Convert a `Tridiagonal` matrix to its dense representation. + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input Tridiagonal matrix. + ${t1}$, allocatable :: B(:, :) + !! Corresponding dense matrix. + end function + #:endfor + end interface + + interface transpose + !! This interface provides methods to compute the transpose operation for + !! the different matrix types defined by `stdlib_specialmatrices`. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#transpose) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function transpose_tridiagonal_${s1}$(A) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(Tridiagonal_${s1}$_type) :: B + end function + #:endfor + end interface + + interface hermitian + !! This interface provides methods to compute the hermitian operation for + !! the different matrix types defined by `stdlib_specialmatrices`. For + !! real-valued matrices, this is equivalent to the standard `transpose`. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#hermitian) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function hermitian_tridiagonal_${s1}$(A) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(Tridiagonal_${s1}$_type) :: B + end function + #:endfor + end interface + + !---------------------------------------- + !----- ----- + !----- ARITHMETIC OPERATORS ----- + !----- ----- + !---------------------------------------- + + interface operator(*) + !! Overload the `*` for scalar-matrix multiplications for the different matrix + !! types provided by `stdlib_specialmatrices`. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operators) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function scalar_multiplication_tridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type) :: B + end function + pure module function scalar_multiplication_bis_tridiagonal_${s1}$(A, alpha) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(Tridiagonal_${s1}$_type) :: B + end function + #:endfor + end interface + + interface operator(+) + !! Overload the `+` operator for matrix-matrix addition. The two matrices need to + !! be of the same type and kind. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operators) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C) + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type), intent(in) :: B + type(Tridiagonal_${s1}$_type) :: C + end function + #:endfor + end interface + + interface operator(-) + !! Overload the `-` operator for matrix-matrix subtraction. The two matrices need to + !! be of the same type and kind. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operators) + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C) + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type), intent(in) :: B + type(Tridiagonal_${s1}$_type) :: C + end function + #:endfor + end interface + +end module stdlib_specialmatrices diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp new file mode 100644 index 000000000..1bfe32f70 --- /dev/null +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -0,0 +1,190 @@ +#:include "common.fypp" +#:set RANKS = range(1, 2+1) +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +submodule (stdlib_specialmatrices) tridiagonal_matrices + use stdlib_linalg_lapack, only: lagtm + contains + + !-------------------------------- + !----- ----- + !----- CONSTRUCTORS ----- + !----- ----- + !-------------------------------- + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A) + !! Construct a `Tridiagonal` matrix from the rank-1 arrays + !! `dl`, `dv` and `du`. + ${t1}$, intent(in) :: dl(:), dv(:), du(:) + !! Tridiagonal matrix elements. + type(Tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + + ! Internal variables. + integer(ilp) :: n + + ! Sanity check. + n = size(dv) + if (size(dl) /= n-1) error stop "Vector dl does not have the correct length." + if (size(du) /= n-1) error stop "Vector du does not have the correct length." + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = dl ; A%dv = dv ; A%du = du + end function + + pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A) + !! Construct a `Tridiagonal` matrix with constant elements. + ${t1}$, intent(in) :: dl, dv, du + !! Tridiagonal matrix elements. + integer(ilp), intent(in) :: n + !! Matrix dimension. + type(Tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + + ! Internal variables. + integer(ilp) :: i + + ! Description of the matrix. + A%n = n + ! Matrix elements. + A%dl = [(dl, i = 1, n-1)] + A%dv = [(dv, i = 1, n)] + A%du = [(du, i = 1, n-1)] + end function + #:endfor + + !----------------------------------------- + !----- ----- + !----- MATRIX-VECTOR PRODUCT ----- + !----- ----- + !----------------------------------------- + + !! spmv_tridiag + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op) + type(Tridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$ + real(${k1}$), intent(in), optional :: alpha + real(${k1}$), intent(in), optional :: beta + character(1), intent(in), optional :: op + + ! Internal variables. + real(${k1}$) :: alpha_, beta_ + integer(ilp) :: n, nrhs, ldx, ldy + character(1) :: op_ + #:if rank == 1 + ${t1}$, pointer :: xmat(:, :), ymat(:, :) + #:endif + + ! Deal with optional arguments. + alpha_ = 1.0_${k1}$ ; if (present(alpha)) alpha_ = alpha + beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta + op_ = "N" ; if (present(op)) op_ = op + + ! Prepare Lapack arguments. + n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$ + nrhs = #{if rank==1}# 1 #{else}# size(x, 2) #{endif}# + + #:if rank == 1 + ! Pointer trick. + xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y + call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy) + #:else + call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy) + #:endif + end subroutine + #:endfor + #:endfor + + !------------------------------------- + !----- ----- + !----- UTILITY FUNCTIONS ----- + !----- ----- + !------------------------------------- + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function tridiagonal_to_dense_${s1}$(A) result(B) + !! Convert a `Tridiagonal` matrix to its dense representation. + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input Tridiagonal matrix. + ${t1}$, allocatable :: B(:, :) + !! Corresponding dense matrix. + + ! Internal variables. + integer(ilp) :: i + + associate (n => A%n) + allocate(B(n, n)) ; B = 0.0_${k1}$ + B(1, 1) = A%dv(1) ; B(1, 2) = A%du(1) + do concurrent (i=2:n-1) + B(i, i-1) = A%dl(i-1) + B(i, i) = A%dv(i) + B(i, i+1) = A%du(i) + enddo + B(n, n-1) = A%dl(n-1) ; B(n, n) = A%dv(n) + end associate + end function + #:endfor + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function transpose_tridiagonal_${s1}$(A) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(Tridiagonal_${s1}$_type) :: B + B%n = A%n ; B%dv = A%dv ; B%du = A%dl ; B%dl = A%du + end function + #:endfor + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function hermitian_tridiagonal_${s1}$(A) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + !! Input matrix. + type(Tridiagonal_${s1}$_type) :: B + B%n = A%n + #:if t1.startswith("complex") + B%dv = conjg(A%dv) ; B%du = conjg(A%dl) ; B%dl = conjg(A%du) + #:else + B%dv = A%dv ; B%du = A%dl ; B%dl = A%du + #:endif + end function + #:endfor + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function scalar_multiplication_tridiagonal_${s1}$(alpha, A) result(B) + ${t1}$, intent(in) :: alpha + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type) :: B + B = Tridiagonal(alpha*A%dl, alpha*A%dv, alpha*A%du) + end function + + pure module function scalar_multiplication_bis_tridiagonal_${s1}$(A, alpha) result(B) + type(Tridiagonal_${s1}$_type), intent(in) :: A + ${t1}$, intent(in) :: alpha + type(Tridiagonal_${s1}$_type) :: B + B = Tridiagonal(alpha*A%dl, alpha*A%dv, alpha*A%du) + end function + #:endfor + + #:for k1, t1, s1 in (KINDS_TYPES) + pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C) + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type), intent(in) :: B + type(Tridiagonal_${s1}$_type) :: C + C = Tridiagonal(A%dl+B%dl, A%dv+B%dv, A%du+B%du) + end function + + pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C) + type(Tridiagonal_${s1}$_type), intent(in) :: A + type(Tridiagonal_${s1}$_type), intent(in) :: B + type(Tridiagonal_${s1}$_type) :: C + C = Tridiagonal(A%dl-B%dl, A%dv-B%dv, A%du-B%du) + end function + #:endfor + +end submodule diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index c496bd2b3..93f340bc9 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -16,6 +16,7 @@ set( "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" "test_linalg_sparse.fypp" + "test_linalg_specialmatrices.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -35,3 +36,4 @@ ADDTEST(linalg_schur) ADDTEST(linalg_svd) ADDTEST(blas_lapack) ADDTEST(linalg_sparse) +ADDTEST(linalg_specialmatrices) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 21237f862..609054a33 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -25,7 +25,7 @@ contains new_unittest('sellc', test_sellc), & new_unittest('symmetries', test_symmetries), & new_unittest('diagonal', test_diagonal), & - new_unittest('add_get_values', test_add_get_values) & + new_unittest('add_get_values', test_add_get_values) & ] end subroutine diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp new file mode 100644 index 000000000..c4a9b45cd --- /dev/null +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -0,0 +1,96 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES +module test_specialmatrices + use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test + use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64 + use stdlib_linalg, only: hermitian + use stdlib_math, only: all_close + use stdlib_specialmatrices + + implicit none + +contains + + + !> Collect all exported unit tests + subroutine collect_suite(testsuite) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest('tridiagonal', test_tridiagonal) & + ] + end subroutine + + subroutine test_tridiagonal(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (KINDS_TYPES) + block + integer, parameter :: wp = ${k1}$ + integer, parameter :: n = 5 + type(Tridiagonal_${s1}$_type) :: A + ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:) + ${t1}$, allocatable :: x(:) + ${t1}$, allocatable :: y1(:), y2(:) + + ! Initialize matrix. + allocate(dl(n-1), dv(n), du(n-1)) + call random_number(dl) ; call random_number(dv) ; call random_number(du) + A = Tridiagonal(dl, dv, du) ; Amat = dense(A) + + ! Random vectors. + allocate(x(n), source = 0.0_wp) ; call random_number(x) + allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp) + + ! Test y = A @ x + y1 = matmul(Amat, x) ; call spmv(A, x, y2) + call check(error, all_close(y1, y2)) + if (allocated(error)) return + + ! Test y = A.T @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(transpose(Amat), x) ; call spmv(A, x, y2, op="T") + call check(error, all_close(y1, y2)) + if (allocated(error)) return + + #:if t1.startswith('complex') + ! Test y = A.H @ x + y1 = 0.0_wp ; y2 = 0.0_wp + y1 = matmul(hermitian(Amat), x) ; call spmv(A, x, y2, op="H") + call check(error, all_close(y1, y2)) + if (allocated(error)) return + #:endif + end block + #:endfor + end subroutine + +end module + + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_specialmatrices, only : collect_suite + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("sparse", collect_suite) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program