From 3912f735904c2100a31607e540706eaacd2c2d0e Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Thu, 20 Mar 2025 09:04:13 +0100 Subject: [PATCH 01/12] Added the Tridiagonal type and associated spmv kernel. --- src/stdlib_sparse_kinds.fypp | 150 ++++++++++++++++++++++++++++++++++- src/stdlib_sparse_spmv.fypp | 54 ++++++++++++- 2 files changed, 201 insertions(+), 3 deletions(-) diff --git a/src/stdlib_sparse_kinds.fypp b/src/stdlib_sparse_kinds.fypp index ceba2a62d..1ebae83b5 100644 --- a/src/stdlib_sparse_kinds.fypp +++ b/src/stdlib_sparse_kinds.fypp @@ -13,6 +13,7 @@ module stdlib_sparse_kinds private public :: sparse_full, sparse_lower, sparse_upper public :: sparse_op_none, sparse_op_transpose, sparse_op_hermitian + public :: Tridiagonal, dense !! version: experimental !! !! Base sparse type holding the meta data related to the storage capacity of a matrix. @@ -128,6 +129,81 @@ module stdlib_sparse_kinds end type #:endfor + !------------------------------------------------------ + !----- ----- + !----- Tridiagonal matrix implementations ----- + !----- ----- + !------------------------------------------------------ + + !! version: experimental + !! + !! Tridiagonal matrix + #:for k1, t1, s1 in (KINDS_TYPES) + type, public, extends(sparse_type) :: Tridiagonal_${s1}$_type + ${t1}$, allocatable :: dl(:), dv(:), du(:) + end type + #:endfor + + interface 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) + module procedure initialize_tridiagonal_${s1}$ + module procedure initialize_constant_tridiagonal_${s1}$ + #:endfor + end interface + + interface dense + !! This interface provides methods to convert a `Tridiagonal` matrix + !! to a regular rank-2 array. + !! + !! #### Syntax + !! + !! ```fortran + !! B = dense(A) + !! ``` + #:for k1, t1, s1 in (KINDS_TYPES) + module procedure tridiagonal_to_dense_${s1}$ + #:endfor + end interface + contains !! (re)Allocate matrix memory for the COO type @@ -589,5 +665,77 @@ contains end subroutine #:endfor + + !------------------------------------------------------ + !----- ----- + !----- Tridiagonal matrix implementations ----- + !----- ----- + !------------------------------------------------------ + + #:for k1, t1, s1 in (KINDS_TYPES) + pure 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%nrows = n ; A%ncols = n ; A%nnz = n + 2*(n-1) + ! Matrix elements. + A%dl = dl ; A%dv = dv ; A%du = du + end function + + pure 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%nrows = n ; A%ncols = n ; A%nnz = n + 2*(n-1) + ! Matrix elements. + A%dl = [(dl, i = 1, n-1)] + A%dv = [(dv, i = 1, n)] + A%du = [(du, i = 1, n-1)] + end function + + pure 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%nrows) + 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 -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..c4ef3a3a7 100644 --- a/src/stdlib_sparse_spmv.fypp +++ b/src/stdlib_sparse_spmv.fypp @@ -13,6 +13,7 @@ module stdlib_sparse_spmv use stdlib_sparse_constants use stdlib_sparse_kinds + use stdlib_linalg_lapack, only: lagtm implicit none private @@ -27,6 +28,9 @@ module stdlib_sparse_spmv module procedure :: spmv_csr_${rank}$d_${s1}$ module procedure :: spmv_csc_${rank}$d_${s1}$ module procedure :: spmv_ell_${rank}$d_${s1}$ + #:if k1 != "qp" and k1 != "xdp" + module procedure :: spmv_tridiag_${rank}$d_${s1}$ + #:endif #:endfor module procedure :: spmv_sellc_${s1}$ #:endfor @@ -589,5 +593,51 @@ contains end subroutine #:endfor - -end module \ No newline at end of file + + !----------------------------------------------------------- + !----- ----- + !----- TRIDIAGONAL MATRIX SPMV IMPLEMENTATIONS ----- + !----- ----- + !----------------------------------------------------------- + !! spmv_tridiag + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + #:if k1 != "qp" and k1 != "xdp" + 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%nrows ; 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 + #:endif + #:endfor + #:endfor + +end module From 872777d18c6c058cb282d67beb2964e1592fba83 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Thu, 20 Mar 2025 09:04:34 +0100 Subject: [PATCH 02/12] Rename dense to dense_mat because of a naming conflict. --- .../linalg/example_sparse_data_accessors.f90 | 84 +++++++++---------- 1 file changed, 42 insertions(+), 42 deletions(-) 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 From abb69042c2ad6c378c238a32b60862c3749c1b4c Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Thu, 20 Mar 2025 09:04:46 +0100 Subject: [PATCH 03/12] Added test for Tridiagonal spmv. --- test/linalg/test_linalg_sparse.fypp | 40 ++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 21237f862..a63f7a110 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -25,7 +25,8 @@ 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), & + new_unittest('tridiagonal', test_tridiagonal) & ] end subroutine @@ -383,6 +384,43 @@ contains #:endfor end subroutine + subroutine test_tridiagonal(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k1, t1, s1 in (KINDS_TYPES) + #:if k1 != "qp" and k1 != "xdp" + 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(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(y1 == y2)) + if (allocated(error)) return + end block + #:endif + #:endfor + end subroutine + end module From 9d8710b4db6cfc88ec6af257814e0f6be87f2395 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Thu, 20 Mar 2025 14:56:57 +0100 Subject: [PATCH 04/12] Changed comparison to account for floating point errors. --- src/stdlib_sparse_spmv.fypp | 6 +++--- test/linalg/test_linalg_sparse.fypp | 14 ++++++++++++-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/stdlib_sparse_spmv.fypp b/src/stdlib_sparse_spmv.fypp index c4ef3a3a7..cb5772baf 100644 --- a/src/stdlib_sparse_spmv.fypp +++ b/src/stdlib_sparse_spmv.fypp @@ -620,9 +620,9 @@ contains #: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 + alpha_ = 1.0_${k1}$ ; if (present(alpha)) alpha_ = alpha + beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta + op_ = sparse_op_none ; if (present(op)) op_ = op ! Prepare Lapack arguments. n = A%nrows ; ldx = n ; ldy = n ; y = 0.0_${k1}$ diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index a63f7a110..fd9cb5240 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -4,6 +4,8 @@ module test_sparse_spmv 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_sparse implicit none @@ -408,14 +410,22 @@ contains ! Test y = A @ x y1 = matmul(Amat, x) ; call spmv(A, x, y2) - call check(error, all(y1 == 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(y1 == y2)) + 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 #:endif #:endfor From 9ae05544da6e709c3fde679ed16f052bf52a35a3 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 24 Mar 2025 08:29:55 +0100 Subject: [PATCH 05/12] Move implementations to a dedicated module. Moving all the implementations to a dedicated module (and submodule) will prevent cluter of the `sparse` module as well as make things clearer from an implementation and documentation point of view. --- src/CMakeLists.txt | 6 +- src/stdlib_sparse_kinds.fypp | 148 ---------------- src/stdlib_sparse_spmv.fypp | 50 ------ src/stdlib_specialmatrices.fypp | 173 +++++++++++++++++++ src/stdlib_specialmatrices_tridiagonal.fypp | 49 ++++++ test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_sparse.fypp | 50 +----- test/linalg/test_linalg_specialmatrices.fypp | 98 +++++++++++ 8 files changed, 327 insertions(+), 249 deletions(-) create mode 100644 src/stdlib_specialmatrices.fypp create mode 100644 src/stdlib_specialmatrices_tridiagonal.fypp create mode 100644 test/linalg/test_linalg_specialmatrices.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 933da34de..fe8f077fd 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_tridiagonal.fypp + stdlib_specialmatrices.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 1ebae83b5..bf03d2747 100644 --- a/src/stdlib_sparse_kinds.fypp +++ b/src/stdlib_sparse_kinds.fypp @@ -13,7 +13,6 @@ module stdlib_sparse_kinds private public :: sparse_full, sparse_lower, sparse_upper public :: sparse_op_none, sparse_op_transpose, sparse_op_hermitian - public :: Tridiagonal, dense !! version: experimental !! !! Base sparse type holding the meta data related to the storage capacity of a matrix. @@ -129,81 +128,6 @@ module stdlib_sparse_kinds end type #:endfor - !------------------------------------------------------ - !----- ----- - !----- Tridiagonal matrix implementations ----- - !----- ----- - !------------------------------------------------------ - - !! version: experimental - !! - !! Tridiagonal matrix - #:for k1, t1, s1 in (KINDS_TYPES) - type, public, extends(sparse_type) :: Tridiagonal_${s1}$_type - ${t1}$, allocatable :: dl(:), dv(:), du(:) - end type - #:endfor - - interface 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) - module procedure initialize_tridiagonal_${s1}$ - module procedure initialize_constant_tridiagonal_${s1}$ - #:endfor - end interface - - interface dense - !! This interface provides methods to convert a `Tridiagonal` matrix - !! to a regular rank-2 array. - !! - !! #### Syntax - !! - !! ```fortran - !! B = dense(A) - !! ``` - #:for k1, t1, s1 in (KINDS_TYPES) - module procedure tridiagonal_to_dense_${s1}$ - #:endfor - end interface - contains !! (re)Allocate matrix memory for the COO type @@ -666,76 +590,4 @@ contains #:endfor - !------------------------------------------------------ - !----- ----- - !----- Tridiagonal matrix implementations ----- - !----- ----- - !------------------------------------------------------ - - #:for k1, t1, s1 in (KINDS_TYPES) - pure 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%nrows = n ; A%ncols = n ; A%nnz = n + 2*(n-1) - ! Matrix elements. - A%dl = dl ; A%dv = dv ; A%du = du - end function - - pure 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%nrows = n ; A%ncols = n ; A%nnz = n + 2*(n-1) - ! Matrix elements. - A%dl = [(dl, i = 1, n-1)] - A%dv = [(dv, i = 1, n)] - A%du = [(du, i = 1, n-1)] - end function - - pure 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%nrows) - 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 - end module stdlib_sparse_kinds diff --git a/src/stdlib_sparse_spmv.fypp b/src/stdlib_sparse_spmv.fypp index cb5772baf..76c78f9e1 100644 --- a/src/stdlib_sparse_spmv.fypp +++ b/src/stdlib_sparse_spmv.fypp @@ -13,7 +13,6 @@ module stdlib_sparse_spmv use stdlib_sparse_constants use stdlib_sparse_kinds - use stdlib_linalg_lapack, only: lagtm implicit none private @@ -28,9 +27,6 @@ module stdlib_sparse_spmv module procedure :: spmv_csr_${rank}$d_${s1}$ module procedure :: spmv_csc_${rank}$d_${s1}$ module procedure :: spmv_ell_${rank}$d_${s1}$ - #:if k1 != "qp" and k1 != "xdp" - module procedure :: spmv_tridiag_${rank}$d_${s1}$ - #:endif #:endfor module procedure :: spmv_sellc_${s1}$ #:endfor @@ -593,51 +589,5 @@ contains end subroutine #:endfor - - !----------------------------------------------------------- - !----- ----- - !----- TRIDIAGONAL MATRIX SPMV IMPLEMENTATIONS ----- - !----- ----- - !----------------------------------------------------------- - !! spmv_tridiag - #:for k1, t1, s1 in (KINDS_TYPES) - #:for rank in RANKS - #:if k1 != "qp" and k1 != "xdp" - 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_ = sparse_op_none ; if (present(op)) op_ = op - - ! Prepare Lapack arguments. - n = A%nrows ; 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 - #:endif - #:endfor - #:endfor end module diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp new file mode 100644 index 000000000..a0111032e --- /dev/null +++ b/src/stdlib_specialmatrices.fypp @@ -0,0 +1,173 @@ +#: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 + use ieee_arithmetic + use stdlib_linalg_constants + implicit none + private + public :: Tridiagonal, spmv, dense + + !! version: experimental + !! + !! Tridiagonal matrix + #:for k1, t1, s1 in (KINDS_TYPES) + type, public :: Tridiagonal_${s1}$_type + ${t1}$, allocatable :: dl(:), dv(:), du(:) + integer(ilp) :: n, m + end type + #:endfor + + interface 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) + module procedure initialize_tridiagonal_${s1}$ + module procedure initialize_constant_tridiagonal_${s1}$ + #:endfor + end interface + + interface spmv + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + #:if k1 != "qp" and k1 != "xdp" + 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 + #:endif + #:endfor + #:endfor + end interface + + interface dense + !! This interface provides methods to convert a `Tridiagonal` matrix + !! to a regular rank-2 array. + !! + !! #### Syntax + !! + !! ```fortran + !! B = dense(A) + !! ``` + #:for k1, t1, s1 in (KINDS_TYPES) + module procedure tridiagonal_to_dense_${s1}$ + #:endfor + end interface + +contains + + !------------------------------------------------------ + !----- ----- + !----- Tridiagonal matrix implementations ----- + !----- ----- + !------------------------------------------------------ + + #: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 ; A%m = 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 ; A%m = 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 + + 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 +end module stdlib_specialmatrices diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp new file mode 100644 index 000000000..df206c8e9 --- /dev/null +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -0,0 +1,49 @@ +#: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 + !! spmv_tridiag + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + #:if k1 != "qp" and k1 != "xdp" + 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 + #:endif + #:endfor + #: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 fd9cb5240..609054a33 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -4,8 +4,6 @@ module test_sparse_spmv 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_sparse implicit none @@ -27,8 +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('tridiagonal', test_tridiagonal) & + new_unittest('add_get_values', test_add_get_values) & ] end subroutine @@ -386,51 +383,6 @@ contains #:endfor end subroutine - subroutine test_tridiagonal(error) - !> Error handling - type(error_type), allocatable, intent(out) :: error - #:for k1, t1, s1 in (KINDS_TYPES) - #:if k1 != "qp" and k1 != "xdp" - 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 - #:endif - #:endfor - end subroutine - end module diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp new file mode 100644 index 000000000..bb05a3a03 --- /dev/null +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -0,0 +1,98 @@ +#: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) + #:if k1 != "qp" and k1 != "xdp" + 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 + #:endif + #: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 From fc61941104dbaa85f825fa618772c0c2c179eb40 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 24 Mar 2025 11:13:48 +0100 Subject: [PATCH 06/12] Implementation of basic operations for Tridiagonal matrices is done. --- src/CMakeLists.txt | 2 +- src/stdlib_specialmatrices.fypp | 220 +++++++++++++------- src/stdlib_specialmatrices_tridiagonal.fypp | 143 +++++++++++++ 3 files changed, 287 insertions(+), 78 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fe8f077fd..58e387104 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -52,8 +52,8 @@ set(fppFiles stdlib_sparse_conversion.fypp stdlib_sparse_kinds.fypp stdlib_sparse_spmv.fypp - stdlib_specialmatrices_tridiagonal.fypp stdlib_specialmatrices.fypp + stdlib_specialmatrices_tridiagonal.fypp stdlib_specialfunctions_gamma.fypp stdlib_stats.fypp stdlib_stats_corr.fypp diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index a0111032e..1275b1460 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -8,18 +8,34 @@ module stdlib_specialmatrices use stdlib_linalg_constants implicit none private - public :: Tridiagonal, spmv, dense - - !! version: experimental + public :: Tridiagonal + public :: spmv + public :: dense, transpose, hermitian + public :: operator(*), operator(+), operator(-) + + !-------------------------------------- + !----- ------ + !----- TYPE DEFINITIONS ------ + !----- ------ + !-------------------------------------- + + !! Version: experimental !! !! Tridiagonal matrix #:for k1, t1, s1 in (KINDS_TYPES) type, public :: Tridiagonal_${s1}$_type + private ${t1}$, allocatable :: dl(:), dv(:), du(:) - integer(ilp) :: n, m + integer(ilp) :: n end type #:endfor + !-------------------------------- + !----- ----- + !----- CONSTRUCTORS ----- + !----- ----- + !-------------------------------- + interface Tridiagonal !! This interface provides different methods to construct a !! `Tridiagonal` matrix. Only the non-zero elements of \( A \) are @@ -61,11 +77,37 @@ module stdlib_specialmatrices !! A = Tridiagonal(a, b, c, n) !! ``` #:for k1, t1, s1 in (KINDS_TYPES) - module procedure initialize_tridiagonal_${s1}$ - module procedure initialize_constant_tridiagonal_${s1}$ + 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 ----- + !----- ----- + !---------------------------------- + + !! Version: experimental + !! + !! Apply the matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_specialmatrices.html#spmv) interface spmv #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS @@ -83,6 +125,16 @@ module stdlib_specialmatrices #:endfor end interface + !------------------------------------- + !----- ----- + !----- UTILITY FUNCTIONS ----- + !----- ----- + !------------------------------------- + + !! Version: experimental + !! + !! Convert a matrix of type `Tridiagonal` to its dense representation. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#dense) interface dense !! This interface provides methods to convert a `Tridiagonal` matrix !! to a regular rank-2 array. @@ -93,81 +145,95 @@ module stdlib_specialmatrices !! B = dense(A) !! ``` #:for k1, t1, s1 in (KINDS_TYPES) - module procedure tridiagonal_to_dense_${s1}$ + 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 -contains + !! Version: experimental + !! + !! Returns the transpose of a `Tridiagonal` matrix. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#transpose) + interface 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 - !------------------------------------------------------ - !----- ----- - !----- Tridiagonal matrix implementations ----- - !----- ----- - !------------------------------------------------------ + !! Version: experimental + !! + !! Returns the Hermitian of a `Tridiagonal` matrix. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#hermitian) + interface 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 - #: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 + !---------------------------------------- + !----- ----- + !----- ARITHMETIC OPERATORS ----- + !----- ----- + !---------------------------------------- + + !! Version: experimental + !! + !! Overloads the scalar multiplication `*` for `Tridiagonal` matrices. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(*)) + interface operator(*) + #: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 + + !! Version: experimental + !! + !! Overloads the addition `+` for `Tridiagonal` matrices. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(+)) + interface operator(+) + #: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 + + !! Version: experimental + !! + !! Overloads the subtraction `-` for `Tridiagonal` matrices. + !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(-)) + interface operator(-) + #: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 - ! 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 ; A%m = 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 ; A%m = 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 - - 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 end module stdlib_specialmatrices diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index df206c8e9..6013d9406 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -6,6 +6,63 @@ 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 @@ -46,4 +103,90 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #:endif #: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 From 9b8f6500e406d1409afe9846e769d7e5184881b1 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Tue, 25 Mar 2025 15:09:42 +0100 Subject: [PATCH 07/12] In-code documentation. --- src/stdlib_specialmatrices.fypp | 77 +++++++++++++++------------------ 1 file changed, 34 insertions(+), 43 deletions(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 1275b1460..e4f41bcd0 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -4,7 +4,10 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES module stdlib_specialmatrices - use ieee_arithmetic + !! 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 @@ -19,11 +22,10 @@ module stdlib_specialmatrices !----- ------ !-------------------------------------- - !! Version: experimental - !! - !! Tridiagonal matrix + !--> 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 @@ -37,9 +39,9 @@ module stdlib_specialmatrices !-------------------------------- interface Tridiagonal - !! This interface provides different methods to construct a - !! `Tridiagonal` matrix. Only the non-zero elements of \( A \) are - !! stored, i.e. + !! ([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 @@ -104,11 +106,13 @@ module stdlib_specialmatrices !----- ----- !---------------------------------- - !! Version: experimental - !! - !! Apply the matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ - !! [Specifications](../page/specs/stdlib_specialmatrices.html#spmv) 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 #:if k1 != "qp" and k1 != "xdp" @@ -131,19 +135,10 @@ module stdlib_specialmatrices !----- ----- !------------------------------------- - !! Version: experimental - !! - !! Convert a matrix of type `Tridiagonal` to its dense representation. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#dense) interface dense - !! This interface provides methods to convert a `Tridiagonal` matrix - !! to a regular rank-2 array. - !! - !! #### Syntax - !! - !! ```fortran - !! B = dense(A) - !! ``` + !! 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. @@ -155,11 +150,10 @@ module stdlib_specialmatrices #:endfor end interface - !! Version: experimental - !! - !! Returns the transpose of a `Tridiagonal` matrix. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#transpose) 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 @@ -169,11 +163,11 @@ module stdlib_specialmatrices #:endfor end interface - !! Version: experimental - !! - !! Returns the Hermitian of a `Tridiagonal` matrix. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#hermitian) 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 @@ -189,11 +183,10 @@ module stdlib_specialmatrices !----- ----- !---------------------------------------- - !! Version: experimental - !! - !! Overloads the scalar multiplication `*` for `Tridiagonal` matrices. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(*)) 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 @@ -208,11 +201,10 @@ module stdlib_specialmatrices #:endfor end interface - !! Version: experimental - !! - !! Overloads the addition `+` for `Tridiagonal` matrices. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(+)) 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 @@ -222,11 +214,10 @@ module stdlib_specialmatrices #:endfor end interface - !! Version: experimental - !! - !! Overloads the subtraction `-` for `Tridiagonal` matrices. - !! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(-)) 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 From de3098b0913a1e070c0f8596a62cbc8d80bcff7e Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Tue, 25 Mar 2025 15:09:49 +0100 Subject: [PATCH 08/12] Added examples. --- .../example_specialmatrices_dp_spmv.f90 | 25 +++++++++++++++++++ .../example_tridiagonal_dp_type.f90 | 18 +++++++++++++ 2 files changed, 43 insertions(+) create mode 100644 example/specialmatrices/example_specialmatrices_dp_spmv.f90 create mode 100644 example/specialmatrices/example_tridiagonal_dp_type.f90 diff --git a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 new file mode 100644 index 000000000..474581edf --- /dev/null +++ b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 @@ -0,0 +1,25 @@ +program example_tridiagonal_matrix + use stdlib_constants, only: dp + use stdlib_specialmatrices + implicit none + + integer, parameter :: n = 5 + type(Tridiagonal_dp_type) :: A + real(dp) :: dl(n), dv(n), du(n) + real(dp) :: x(n), y(n), y_dense(n) + + ! Create an arbitrary tridiagonal matrix. + dl = [(i, i=1, n - 1)]; dv = [(2*i, i=1, n)]; du = [(3*i, i=1, n)] + 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..e1edee1c3 --- /dev/null +++ b/example/specialmatrices/example_tridiagonal_dp_type.f90 @@ -0,0 +1,18 @@ +program example_tridiagonal_matrix + use stdlib_constants, only: dp + use stdlib_specialmatrices + implicit none + + integer, parameter :: n = 5 + type(Tridiagonal_dp_type) :: A + real(dp) :: dl(n), dv(n), du(n) + + ! 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 From 6ead4bad96e5e243f52a69937741d34d5634d148 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Tue, 25 Mar 2025 15:10:06 +0100 Subject: [PATCH 09/12] Specifications page. --- doc/specs/stdlib_specialmatrices.md | 215 ++++++++++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 doc/specs/stdlib_specialmatrices.md 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 From 46986625125530cad6e6a94189998c5ac74fba2e Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Tue, 25 Mar 2025 15:29:19 +0100 Subject: [PATCH 10/12] Fix typos in module header. --- src/stdlib_specialmatrices.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index e4f41bcd0..98d9395a2 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -7,7 +7,7 @@ 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)) + !! applications. ([Specifications](../page/specs/stdlib_specialmatrices.html)) use stdlib_linalg_constants implicit none private From c9d25319f72d317eaa3aa2f66cf1badf61995b7b Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Tue, 25 Mar 2025 15:54:46 +0100 Subject: [PATCH 11/12] Fix errors in examples. --- example/CMakeLists.txt | 1 + example/specialmatrices/CMakeLists.txt | 2 ++ .../specialmatrices/example_specialmatrices_dp_spmv.f90 | 7 ++++--- example/specialmatrices/example_tridiagonal_dp_type.f90 | 4 ++-- 4 files changed, 9 insertions(+), 5 deletions(-) create mode 100644 example/specialmatrices/CMakeLists.txt 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/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 index 474581edf..ebfcf34d0 100644 --- a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 +++ b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 @@ -1,15 +1,16 @@ program example_tridiagonal_matrix - use stdlib_constants, only: dp + use stdlib_linalg_constants, only: dp use stdlib_specialmatrices implicit none integer, parameter :: n = 5 type(Tridiagonal_dp_type) :: A - real(dp) :: dl(n), dv(n), du(n) + 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)] + 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. diff --git a/example/specialmatrices/example_tridiagonal_dp_type.f90 b/example/specialmatrices/example_tridiagonal_dp_type.f90 index e1edee1c3..7f221b5d6 100644 --- a/example/specialmatrices/example_tridiagonal_dp_type.f90 +++ b/example/specialmatrices/example_tridiagonal_dp_type.f90 @@ -1,11 +1,11 @@ program example_tridiagonal_matrix - use stdlib_constants, only: dp + use stdlib_linalg_constants, only: dp use stdlib_specialmatrices implicit none integer, parameter :: n = 5 type(Tridiagonal_dp_type) :: A - real(dp) :: dl(n), dv(n), du(n) + real(dp) :: dl(n - 1), dv(n), du(n - 1) ! Generate random tridiagonal elements. call random_number(dl) From 9e91af3026efa452bfedf08810e2891c9e91f356 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 26 Mar 2025 12:27:36 +0100 Subject: [PATCH 12/12] Enabled `xdp` and `qp` for `spmv` kernels. --- src/stdlib_specialmatrices.fypp | 2 -- src/stdlib_specialmatrices_tridiagonal.fypp | 2 -- test/linalg/test_linalg_specialmatrices.fypp | 2 -- 3 files changed, 6 deletions(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 98d9395a2..d37b86a4a 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -115,7 +115,6 @@ module stdlib_specialmatrices !! for the different matrix types defined by `stdlib_specialmatrices`. #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - #:if k1 != "qp" and k1 != "xdp" 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)}$ @@ -124,7 +123,6 @@ module stdlib_specialmatrices real(${k1}$), intent(in), optional :: beta character(1), intent(in), optional :: op end subroutine - #:endif #:endfor #:endfor end interface diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 6013d9406..1bfe32f70 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -66,7 +66,6 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices !! spmv_tridiag #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - #:if k1 != "qp" and k1 != "xdp" 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)}$ @@ -100,7 +99,6 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy) #:endif end subroutine - #:endif #:endfor #:endfor diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index bb05a3a03..c4a9b45cd 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -27,7 +27,6 @@ contains !> Error handling type(error_type), allocatable, intent(out) :: error #:for k1, t1, s1 in (KINDS_TYPES) - #:if k1 != "qp" and k1 != "xdp" block integer, parameter :: wp = ${k1}$ integer, parameter :: n = 5 @@ -64,7 +63,6 @@ contains if (allocated(error)) return #:endif end block - #:endif #:endfor end subroutine