From 3912f735904c2100a31607e540706eaacd2c2d0e Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Thu, 20 Mar 2025 09:04:13 +0100 Subject: [PATCH 01/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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/27] 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 From 78133fbcf8466849744b7c701771e996f19cb5a7 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 13 Jun 2025 14:13:03 +0200 Subject: [PATCH 13/27] Update test/linalg/test_linalg_sparse.fypp Remove trailing white space. Co-authored-by: Jeremie Vandenplas --- test/linalg/test_linalg_sparse.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 609054a33..21237f862 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 From 78b50859bd9dbcc952744405d3b93372e9e0a840 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 13 Jun 2025 14:13:23 +0200 Subject: [PATCH 14/27] Update example/specialmatrices/example_specialmatrices_dp_spmv.f90 Import only required elements from the module. Co-authored-by: Jeremie Vandenplas --- example/specialmatrices/example_specialmatrices_dp_spmv.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 index ebfcf34d0..7313672b6 100644 --- a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 +++ b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 @@ -1,6 +1,6 @@ program example_tridiagonal_matrix use stdlib_linalg_constants, only: dp - use stdlib_specialmatrices + use stdlib_specialmatrices, only: tridiaonal_dp_type, tridiagonal, dense implicit none integer, parameter :: n = 5 From 6777aa0ea272ff1e65350ec24ce42534b0fa70b4 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 13 Jun 2025 14:53:35 +0200 Subject: [PATCH 15/27] Consistent lower-casing and imports. --- doc/specs/stdlib_specialmatrices.md | 26 ++++---- .../example_specialmatrices_dp_spmv.f90 | 6 +- src/stdlib_specialmatrices.fypp | 60 +++++++++---------- test/linalg/test_linalg_specialmatrices.fypp | 11 +++- 4 files changed, 54 insertions(+), 49 deletions(-) diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md index d73866628..d688c7ce6 100644 --- a/doc/specs/stdlib_specialmatrices.md +++ b/doc/specs/stdlib_specialmatrices.md @@ -33,7 +33,7 @@ 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 generic tridiagonal matrix has the following structure: $$ A = @@ -50,25 +50,25 @@ This particular structure also lends itself to specialized implementations for m 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. +- `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)` +`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)` +`A = ` [[stdlib_specialmatrices(module):tridiagonal(interface)]] `(dl, dv, du, n)` #### Example @@ -90,7 +90,7 @@ Experimental 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. +- For `tridiagonal` matrices, the LAPACK `lagtm` backend is being used. #### Syntax @@ -111,7 +111,7 @@ With the exception of `extended precision` and `quadruple precision`, all the ty - `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. +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 diff --git a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 index 7313672b6..6a3c4e0ab 100644 --- a/example/specialmatrices/example_specialmatrices_dp_spmv.f90 +++ b/example/specialmatrices/example_specialmatrices_dp_spmv.f90 @@ -1,17 +1,17 @@ program example_tridiagonal_matrix use stdlib_linalg_constants, only: dp - use stdlib_specialmatrices, only: tridiaonal_dp_type, tridiagonal, dense + use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal, dense, spmv implicit none integer, parameter :: n = 5 - type(Tridiagonal_dp_type) :: A + 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) + A = tridiagonal(dl, dv, du) ! Initialize vectors. x = 1.0_dp; y = 0.0_dp; y_dense = 0.0_dp diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index d37b86a4a..df821bdc7 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -11,7 +11,7 @@ module stdlib_specialmatrices use stdlib_linalg_constants implicit none private - public :: Tridiagonal + public :: tridiagonal public :: spmv public :: dense, transpose, hermitian public :: operator(*), operator(+), operator(-) @@ -24,8 +24,8 @@ module stdlib_specialmatrices !--> Tridiagonal matrices #:for k1, t1, s1 in (KINDS_TYPES) - type, public :: Tridiagonal_${s1}$_type - !! Base type to define a `Tridiagonal` matrix. + type, public :: tridiagonal_${s1}$_type + !! Base type to define a `tridiagonal` matrix. private ${t1}$, allocatable :: dl(:), dv(:), du(:) integer(ilp) :: n @@ -38,9 +38,9 @@ module stdlib_specialmatrices !----- ----- !-------------------------------- - interface Tridiagonal + interface tridiagonal !! ([Specifications](../page/specs/stdlib_specialmatrices.html#Tridiagonal)) This - !! interface provides different methods to construct a `Tridiagonal` matrix. Only + !! interface provides different methods to construct a `tridiagonal` matrix. Only !! the non-zero elements of \( A \) are stored, i.e. !! !! \[ @@ -57,44 +57,44 @@ module stdlib_specialmatrices !! !! #### Syntax !! - !! - Construct a real `Tridiagonal` matrix from rank-1 arrays: + !! - Construct a real `tridiagonal` matrix from rank-1 arrays: !! !! ```fortran !! integer, parameter :: n !! real(dp), allocatable :: dl(:), dv(:), du(:) - !! type(Tridiagonal_rdp_type) :: A + !! 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: + !! - 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 + !! 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 + !! 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 + 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. + !! 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 + type(tridiagonal_${s1}$_type) :: A !! Corresponding Tridiagonal matrix. end function #:endfor @@ -116,7 +116,7 @@ module 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 + 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 @@ -139,8 +139,8 @@ module stdlib_specialmatrices !! ([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 + !! Convert a `tridiagonal` matrix to its dense representation. + type(tridiagonal_${s1}$_type), intent(in) :: A !! Input Tridiagonal matrix. ${t1}$, allocatable :: B(:, :) !! Corresponding dense matrix. @@ -154,9 +154,9 @@ module 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 + type(tridiagonal_${s1}$_type), intent(in) :: A !! Input matrix. - type(Tridiagonal_${s1}$_type) :: B + type(tridiagonal_${s1}$_type) :: B end function #:endfor end interface @@ -168,9 +168,9 @@ module stdlib_specialmatrices !! [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 + type(tridiagonal_${s1}$_type), intent(in) :: A !! Input matrix. - type(Tridiagonal_${s1}$_type) :: B + type(tridiagonal_${s1}$_type) :: B end function #:endfor end interface @@ -188,13 +188,13 @@ module stdlib_specialmatrices #: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 + 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 + type(tridiagonal_${s1}$_type), intent(in) :: A ${t1}$, intent(in) :: alpha - type(Tridiagonal_${s1}$_type) :: B + type(tridiagonal_${s1}$_type) :: B end function #:endfor end interface @@ -205,9 +205,9 @@ module stdlib_specialmatrices !! [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 + type(tridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C end function #:endfor end interface @@ -218,9 +218,9 @@ module stdlib_specialmatrices !! [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 + type(tridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C end function #:endfor end interface diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index c4a9b45cd..205e4c63b 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -6,7 +6,12 @@ module test_specialmatrices 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 + use stdlib_specialmatrices, only: tridiagonal, & + tridiagonal_sp_type, & + tridiagonal_dp_type, & + tridiagonal_xdp_type, & + tridiagonal_qp_type, & + dense, spmv implicit none @@ -30,7 +35,7 @@ contains block integer, parameter :: wp = ${k1}$ integer, parameter :: n = 5 - type(Tridiagonal_${s1}$_type) :: A + type(tridiagonal_${s1}$_type) :: A ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:) ${t1}$, allocatable :: x(:) ${t1}$, allocatable :: y1(:), y2(:) @@ -38,7 +43,7 @@ contains ! 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) + A = tridiagonal(dl, dv, du) ; Amat = dense(A) ! Random vectors. allocate(x(n), source = 0.0_wp) ; call random_number(x) From 48c734ef8384227d16b09aef4ee6c8b43d2a72fe Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Sun, 22 Jun 2025 14:24:07 +0200 Subject: [PATCH 16/27] Fixed allocation from source for complex types. --- src/stdlib_specialmatrices.fypp | 17 +++++++++++++++++ src/stdlib_specialmatrices_tridiagonal.fypp | 6 +++++- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index df821bdc7..6a4e63f97 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -16,6 +16,23 @@ module stdlib_specialmatrices public :: dense, transpose, hermitian public :: operator(*), operator(+), operator(-) + real(sp), parameter :: zero_sp = 0._sp + real(sp), parameter :: one_sp = 1._sp + real(dp), parameter :: zero_dp = 0._dp + real(dp), parameter :: one_dp = 1._dp + real(xdp), parameter :: zero_xdp = 0._xdp + real(xdp), parameter :: one_xdp = 1._xdp + real(qp), parameter :: zero_qp = 0._qp + real(qp), parameter :: one_qp = 1._qp + complex(sp), parameter :: zero_csp = (0._sp,0._sp) + complex(sp), parameter :: one_csp = (1._sp,1._sp) + complex(dp), parameter :: zero_cdp = (0._dp,0._dp) + complex(dp), parameter :: one_cdp = (1._dp,1._dp) + complex(xdp), parameter :: zero_cxdp = (0._xdp,0._xdp) + complex(xdp), parameter :: one_cxdp = (1._xdp,1._xdp) + complex(qp), parameter :: zero_cqp = (0._qp,0._qp) + complex(qp), parameter :: one_cqp = (1._qp,1._qp) + !-------------------------------------- !----- ------ !----- TYPE DEFINITIONS ------ diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 1bfe32f70..55f9e11a1 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -120,7 +120,11 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices integer(ilp) :: i associate (n => A%n) - allocate(B(n, n)) ; B = 0.0_${k1}$ + #:if t1.startswith('complex') + allocate(B(n, n), source=zero_c${k1}$) + #:else + allocate(B(n, n), source=zero_${k1}$) + #:endif 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) From f102167f745851fd653373504f0950be88d16291 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 27 Jun 2025 09:53:57 +0200 Subject: [PATCH 17/27] Make use of stdlib_constants to define the various zeros and ones. --- src/stdlib_specialmatrices.fypp | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 6a4e63f97..939eed294 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -9,6 +9,7 @@ module stdlib_specialmatrices !! of partial differential equations, as well as control and signal processing !! applications. ([Specifications](../page/specs/stdlib_specialmatrices.html)) use stdlib_linalg_constants + use stdlib_constants implicit none private public :: tridiagonal @@ -16,22 +17,22 @@ module stdlib_specialmatrices public :: dense, transpose, hermitian public :: operator(*), operator(+), operator(-) - real(sp), parameter :: zero_sp = 0._sp - real(sp), parameter :: one_sp = 1._sp - real(dp), parameter :: zero_dp = 0._dp - real(dp), parameter :: one_dp = 1._dp - real(xdp), parameter :: zero_xdp = 0._xdp - real(xdp), parameter :: one_xdp = 1._xdp - real(qp), parameter :: zero_qp = 0._qp - real(qp), parameter :: one_qp = 1._qp - complex(sp), parameter :: zero_csp = (0._sp,0._sp) - complex(sp), parameter :: one_csp = (1._sp,1._sp) - complex(dp), parameter :: zero_cdp = (0._dp,0._dp) - complex(dp), parameter :: one_cdp = (1._dp,1._dp) - complex(xdp), parameter :: zero_cxdp = (0._xdp,0._xdp) - complex(xdp), parameter :: one_cxdp = (1._xdp,1._xdp) - complex(qp), parameter :: zero_cqp = (0._qp,0._qp) - complex(qp), parameter :: one_cqp = (1._qp,1._qp) + ! real(sp), parameter :: zero_sp = 0._sp + ! real(sp), parameter :: one_sp = 1._sp + ! real(dp), parameter :: zero_dp = 0._dp + ! real(dp), parameter :: one_dp = 1._dp + ! real(xdp), parameter :: zero_xdp = 0._xdp + ! real(xdp), parameter :: one_xdp = 1._xdp + ! real(qp), parameter :: zero_qp = 0._qp + ! real(qp), parameter :: one_qp = 1._qp + ! complex(sp), parameter :: zero_csp = (0._sp,0._sp) + ! complex(sp), parameter :: one_csp = (1._sp,1._sp) + ! complex(dp), parameter :: zero_cdp = (0._dp,0._dp) + ! complex(dp), parameter :: one_cdp = (1._dp,1._dp) + ! complex(xdp), parameter :: zero_cxdp = (0._xdp,0._xdp) + ! complex(xdp), parameter :: one_cxdp = (1._xdp,1._xdp) + ! complex(qp), parameter :: zero_cqp = (0._qp,0._qp) + ! complex(qp), parameter :: one_cqp = (1._qp,1._qp) !-------------------------------------- !----- ------ From 05be1d613704229b619ad8b61ba39719bb651909 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 27 Jun 2025 09:54:29 +0200 Subject: [PATCH 18/27] Prevent temporary arrays + lower casing every names. --- src/stdlib_specialmatrices_tridiagonal.fypp | 71 +++++++++++---------- 1 file changed, 37 insertions(+), 34 deletions(-) diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 55f9e11a1..22a42d4df 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -15,12 +15,12 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #: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 + !! 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. + !! tridiagonal matrix elements. + type(tridiagonal_${s1}$_type) :: A + !! Corresponding tridiagonal matrix. ! Internal variables. integer(ilp) :: n @@ -37,13 +37,13 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices end function pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A) - !! Construct a `Tridiagonal` matrix with constant elements. + !! Construct a `tridiagonal` matrix with constant elements. ${t1}$, intent(in) :: dl, dv, du - !! Tridiagonal matrix elements. + !! tridiagonal matrix elements. integer(ilp), intent(in) :: n !! Matrix dimension. - type(Tridiagonal_${s1}$_type) :: A - !! Corresponding Tridiagonal matrix. + type(tridiagonal_${s1}$_type) :: A + !! Corresponding tridiagonal matrix. ! Internal variables. integer(ilp) :: i @@ -67,7 +67,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #: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 + 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 @@ -110,9 +110,9 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #: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. + !! Convert a `tridiagonal` matrix to its dense representation. + type(tridiagonal_${s1}$_type), intent(in) :: A + !! Input tridiagonal matrix. ${t1}$, allocatable :: B(:, :) !! Corresponding dense matrix. @@ -138,23 +138,22 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #:for k1, t1, s1 in (KINDS_TYPES) pure module function transpose_tridiagonal_${s1}$(A) result(B) - type(Tridiagonal_${s1}$_type), intent(in) :: A + 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 + type(tridiagonal_${s1}$_type) :: B + B = tridiagonal(A%du, A%dv, A%dl) 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 + type(tridiagonal_${s1}$_type), intent(in) :: A !! Input matrix. - type(Tridiagonal_${s1}$_type) :: B - B%n = A%n + type(tridiagonal_${s1}$_type) :: B #:if t1.startswith("complex") - B%dv = conjg(A%dv) ; B%du = conjg(A%dl) ; B%dl = conjg(A%du) + B = tridiagonal(conjg(A%du), conjg(A%dv), conjg(A%dl)) #:else - B%dv = A%dv ; B%du = A%dl ; B%dl = A%du + B = tridiagonal(A%du, A%dv, A%dl) #:endif end function #:endfor @@ -162,32 +161,36 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices #: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) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type) :: B + B = tridiagonal(A%dl, A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = alpha*B%du end function pure module function scalar_multiplication_bis_tridiagonal_${s1}$(A, alpha) result(B) - type(Tridiagonal_${s1}$_type), intent(in) :: A + 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) + type(tridiagonal_${s1}$_type) :: B + B = tridiagonal(A%dl, A%dv, A%du) + B%dl = alpha*B%dl; B%dv = alpha*B%dv; B%du = alpha*B%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) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl + B%dl; C%dv = C%dv + B%dv; C%du = C%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) + type(tridiagonal_${s1}$_type), intent(in) :: A + type(tridiagonal_${s1}$_type), intent(in) :: B + type(tridiagonal_${s1}$_type) :: C + C = tridiagonal(A%dl, A%dv, A%du) + C%dl = C%dl - B%dl; C%dv = C%dv - B%dv; C%du = C%du - B%du end function #:endfor From 94220f4b17181c7c7613daf3095e45be4e50acf0 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 27 Jun 2025 09:58:49 +0200 Subject: [PATCH 19/27] Import everything from stdlib_kinds --- test/linalg/test_linalg_specialmatrices.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index 205e4c63b..cf7af63f3 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -3,7 +3,7 @@ #: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_kinds use stdlib_linalg, only: hermitian use stdlib_math, only: all_close use stdlib_specialmatrices, only: tridiagonal, & From b4b4f8756252ea663d86e48e374a56f06f45046e Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 27 Jun 2025 10:39:46 +0200 Subject: [PATCH 20/27] Fix imports. --- test/linalg/test_linalg_specialmatrices.fypp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index cf7af63f3..3fa2d2038 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -6,13 +6,7 @@ module test_specialmatrices use stdlib_kinds use stdlib_linalg, only: hermitian use stdlib_math, only: all_close - use stdlib_specialmatrices, only: tridiagonal, & - tridiagonal_sp_type, & - tridiagonal_dp_type, & - tridiagonal_xdp_type, & - tridiagonal_qp_type, & - dense, spmv - + use stdlib_specialmatrices implicit none contains From b8366a1d0ce34d96fe08f0de5bb836414d0a7fb1 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Fri, 27 Jun 2025 13:15:35 +0200 Subject: [PATCH 21/27] Fix missing argument in check. --- test/linalg/test_linalg_specialmatrices.fypp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index 3fa2d2038..2748fcda1 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/test/linalg/test_linalg_specialmatrices.fypp @@ -45,20 +45,20 @@ contains ! Test y = A @ x y1 = matmul(Amat, x) ; call spmv(A, x, y2) - call check(error, all_close(y1, y2)) + call check(error, all_close(y1, y2), .true.) 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)) + call check(error, all_close(y1, y2), .true.) 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)) + call check(error, all_close(y1, y2), .true.) if (allocated(error)) return #:endif end block From 2fcc330f806464a13fb154a12985eae569062e00 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 16:47:16 +0200 Subject: [PATCH 22/27] Update doc/specs/stdlib_specialmatrices.md Co-authored-by: Federico Perini --- doc/specs/stdlib_specialmatrices.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md index d688c7ce6..59d1218cd 100644 --- a/doc/specs/stdlib_specialmatrices.md +++ b/doc/specs/stdlib_specialmatrices.md @@ -48,7 +48,7 @@ $$ 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 matrices are available with all supported data types as `tridiagonal__type`, for example: - `tridiagonal_sp_type` : Tridiagonal matrix of size `n` with `real`/`single precision` data. - `tridiagonal_dp_type` : Tridiagonal matrix of size `n` with `real`/`double precision` data. From 94ad82e52aa868992e1464eab00edee989ab3cc0 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 16:47:43 +0200 Subject: [PATCH 23/27] Update doc/specs/stdlib_specialmatrices.md Co-authored-by: Federico Perini --- doc/specs/stdlib_specialmatrices.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md index 59d1218cd..fbae1f99f 100644 --- a/doc/specs/stdlib_specialmatrices.md +++ b/doc/specs/stdlib_specialmatrices.md @@ -111,7 +111,7 @@ With the exception of `extended precision` and `quadruple precision`, all the ty - `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. +Due to limitations of the underlying `lapack` driver, currently `alpha` and `beta` can only take one of the values `[-1, 0, 1]` for `tridiagonal` and `symtridiagonal` matrices. See `lagtm` for more details. @endwarning #### Examples From 2fbda5d7d546ed3d9e4feb8b4473bde126838ea6 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 16:48:16 +0200 Subject: [PATCH 24/27] Update doc/specs/stdlib_specialmatrices.md Co-authored-by: Federico Perini --- doc/specs/stdlib_specialmatrices.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialmatrices.md b/doc/specs/stdlib_specialmatrices.md index fbae1f99f..6c86f3bd1 100644 --- a/doc/specs/stdlib_specialmatrices.md +++ b/doc/specs/stdlib_specialmatrices.md @@ -211,5 +211,5 @@ The definition of all standard artihmetic operators have been overloaded to be a `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. +For addition (`+`) and subtraction (`-`), matrices `A`, `B` and `C` all need to be of the same type and kind. For scalar multiplication (`*`), `A` and `B` need to be of the same type and kind, while `alpha` is either `real` or `complex` (with the same kind again) depending on the type being used. @endnote From 1881c44aecb4d8cd47320745544c125f0d1c3470 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 16:48:36 +0200 Subject: [PATCH 25/27] Update src/stdlib_specialmatrices.fypp Co-authored-by: Federico Perini --- src/stdlib_specialmatrices.fypp | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 939eed294..8f1b7617a 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -17,23 +17,6 @@ module stdlib_specialmatrices public :: dense, transpose, hermitian public :: operator(*), operator(+), operator(-) - ! real(sp), parameter :: zero_sp = 0._sp - ! real(sp), parameter :: one_sp = 1._sp - ! real(dp), parameter :: zero_dp = 0._dp - ! real(dp), parameter :: one_dp = 1._dp - ! real(xdp), parameter :: zero_xdp = 0._xdp - ! real(xdp), parameter :: one_xdp = 1._xdp - ! real(qp), parameter :: zero_qp = 0._qp - ! real(qp), parameter :: one_qp = 1._qp - ! complex(sp), parameter :: zero_csp = (0._sp,0._sp) - ! complex(sp), parameter :: one_csp = (1._sp,1._sp) - ! complex(dp), parameter :: zero_cdp = (0._dp,0._dp) - ! complex(dp), parameter :: one_cdp = (1._dp,1._dp) - ! complex(xdp), parameter :: zero_cxdp = (0._xdp,0._xdp) - ! complex(xdp), parameter :: one_cxdp = (1._xdp,1._xdp) - ! complex(qp), parameter :: zero_cqp = (0._qp,0._qp) - ! complex(qp), parameter :: one_cqp = (1._qp,1._qp) - !-------------------------------------- !----- ------ !----- TYPE DEFINITIONS ------ From 9d37af7cc7559a3667ea7afb26ce0534d837c653 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 16:48:49 +0200 Subject: [PATCH 26/27] Update src/stdlib_specialmatrices_tridiagonal.fypp Co-authored-by: Federico Perini --- src/stdlib_specialmatrices_tridiagonal.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 22a42d4df..69e9b12c7 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -89,7 +89,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices ! Prepare Lapack arguments. n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$ - nrhs = #{if rank==1}# 1 #{else}# size(x, 2) #{endif}# + nrhs = #{if rank==1}# 1 #{else}# size(x, dim=2, kind=ilp) #{endif}# #:if rank == 1 ! Pointer trick. From 29a6e04717c94dd463c1d9cd5296ae2357b7c385 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 7 Jul 2025 18:26:27 +0200 Subject: [PATCH 27/27] Improved error handling. --- src/stdlib_specialmatrices.fypp | 29 ++++++- src/stdlib_specialmatrices_tridiagonal.fypp | 90 +++++++++++++++++++-- 2 files changed, 112 insertions(+), 7 deletions(-) diff --git a/src/stdlib_specialmatrices.fypp b/src/stdlib_specialmatrices.fypp index 8f1b7617a..4db497fe1 100644 --- a/src/stdlib_specialmatrices.fypp +++ b/src/stdlib_specialmatrices.fypp @@ -10,6 +10,8 @@ module stdlib_specialmatrices !! applications. ([Specifications](../page/specs/stdlib_specialmatrices.html)) use stdlib_linalg_constants use stdlib_constants + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none private public :: tridiagonal @@ -80,7 +82,7 @@ module stdlib_specialmatrices !! A = Tridiagonal(a, b, c, n) !! ``` #:for k1, t1, s1 in (KINDS_TYPES) - pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A) + pure module function initialize_tridiagonal_pure_${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(:) @@ -89,7 +91,7 @@ module stdlib_specialmatrices !! Corresponding Tridiagonal matrix. end function - pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A) + pure module function initialize_constant_tridiagonal_pure_${s1}$(dl, dv, du, n) result(A) !! Construct a `tridiagonal` matrix with constant elements. ${t1}$, intent(in) :: dl, dv, du !! Tridiagonal matrix elements. @@ -98,6 +100,29 @@ module stdlib_specialmatrices type(tridiagonal_${s1}$_type) :: A !! Corresponding Tridiagonal matrix. end function + + module function initialize_tridiagonal_impure_${s1}$(dl, dv, du, err) 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(linalg_state_type), intent(out) :: err + !! Error handling. + type(tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + end function + + module function initialize_constant_tridiagonal_impure_${s1}$(dl, dv, du, n, err) 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(linalg_state_type), intent(out) :: err + !! Error handling. + type(tridiagonal_${s1}$_type) :: A + !! Corresponding Tridiagonal matrix. + end function #:endfor end interface diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 69e9b12c7..8fabaa5df 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -5,6 +5,8 @@ #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule (stdlib_specialmatrices) tridiagonal_matrices use stdlib_linalg_lapack, only: lagtm + + character(len=*), parameter :: this = "tridiagonal matrices" contains !-------------------------------- @@ -14,21 +16,92 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices !-------------------------------- #:for k1, t1, s1 in (KINDS_TYPES) - pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A) + pure module function initialize_tridiagonal_pure_${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 + type(linalg_state_type) :: err0 + + ! Sanity check. + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + if (size(dl, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector dl does not have the correct length.") + call linalg_error_handling(err0) + endif + if (size(du, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector du does not have the correct length.") + call linalg_error_handling(err0) + endif + + ! 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_pure_${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 + type(linalg_state_type) :: err0 + + ! Description of the matrix. + A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0) + endif + ! Matrix elements. + A%dl = [(dl, i = 1, n-1)] + A%dv = [(dv, i = 1, n)] + A%du = [(du, i = 1, n-1)] + end function + + module function initialize_tridiagonal_impure_${s1}$(dl, dv, du, err) 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(linalg_state_type), intent(out) :: err + !! Error handling. type(tridiagonal_${s1}$_type) :: A !! Corresponding tridiagonal matrix. ! Internal variables. integer(ilp) :: n + type(linalg_state_type) :: err0 ! 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." + n = size(dv, kind=ilp) + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif + if (size(dl, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector dl does not have the correct length.") + call linalg_error_handling(err0, err) + endif + if (size(du, kind=ilp) /= n-1) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Vector du does not have the correct length.") + call linalg_error_handling(err0, err) + endif ! Description of the matrix. A%n = n @@ -36,20 +109,27 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices A%dl = dl ; A%dv = dv ; A%du = du end function - pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A) + module function initialize_constant_tridiagonal_impure_${s1}$(dl, dv, du, n, err) 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(linalg_state_type), intent(out) :: err + !! Error handling type(tridiagonal_${s1}$_type) :: A !! Corresponding tridiagonal matrix. ! Internal variables. integer(ilp) :: i + type(linalg_state_type) :: err0 ! Description of the matrix. A%n = n + if (n <= 0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".") + call linalg_error_handling(err0, err) + endif ! Matrix elements. A%dl = [(dl, i = 1, n-1)] A%dv = [(dv, i = 1, n)]