Base type and constructor interfaces for symtridiagonal matrices.

src/stdlib_specialmatrices.fypp

@@ -35,6 +35,16 @@ module stdlib_specialmatrices
  end type
  #:endfor
 
+ !--> Symmetric Tridiagonal matrices
+ #:for k1, t1, s1 in (KINDS_TYPES)
+ type, public :: symtridiagonal_${s1}$_type
+ !! Base type to define a `symtridiagonal` matrix.
+ private
+ ${t1}$, allocatable :: dv(:), ev(:)
+ integer(ilp) :: n
+ end type
+ #:endfor
+
  !--------------------------------
  !----- -----
  !----- CONSTRUCTORS -----
@@ -126,6 +136,91 @@ module stdlib_specialmatrices
  #:endfor
  end interface
 
+ interface symtridiagonal
+ !! ([Specifications](../page/specs/stdlib_specialmatrices.html#SymTridiagonal)) This
+ !! interface provides different methods to construct a `sylmtridiagonal`
+ !! matrix. Only the non-zero elements of \( A \) are stored, i.e.
+ !!
+ !! \[
+ !! A
+ !! =
+ !! \begin{bmatrix}
+ !! a_1 & b_1 \\
+ !! b_1 & a_2 & b_2 \\
+ !! & \ddots & \ddots & \ddots \\
+ !! & & b_{n-2} & a_{n-1} & b_{n-1} \\
+ !! & & & b_{n-1} & a_n
+ !! \end{bmatrix}.
+ !! \]
+ !!
+ !! #### Syntax
+ !!
+ !! - Construct a real `symtridiagonal` matrix from rank-1 arrays:
+ !!
+ !! ```fortran
+ !! integer, parameter :: n
+ !! real(dp), allocatable :: dv(:), ev(:)
+ !! type(symtridiagonal_rdp_type) :: A
+ !! integer :: i
+ !!
+ !! ev = [(i, i=1, n-1)]; dv = [(2*i, i=1, n)]
+ !! A = SymTridiagonal(dv, ev)
+ !! ```
+ !!
+ !! - Construct a real `symtridiagonal` matrix with constant diagonals:
+ !!
+ !! ```fortran
+ !! integer, parameter :: n
+ !! real(dp), parameter :: a = 1.0_dp, b = 1.0_dp
+ !! type(symtridiagonal_rdp_type) :: A
+ !!
+ !! A = SymTridiagonal(a, b, n)
+ !! ```
+ #:for k1, t1, s1 in (KINDS_TYPES)
+ pure module function initialize_symtridiagonal_pure_${s1}$(dv, ev) result(A)
+ !! Construct a `tridiagonal` matrix from the rank-1 arrays
+ !! `dl`, `dv` and `du`.
+ ${t1}$, intent(in) :: dv(:), ev(:)
+ !! SymTridiagonal matrix elements.
+ type(symtridiagonal_${s1}$_type) :: A
+ !! Corresponding SymTridiagonal matrix.
+ end function
+
+ pure module function initialize_constant_symtridiagonal_pure_${s1}$(dv, ev, n) result(A)
+ !! Construct a `symtridiagonal` matrix with constant elements.
+ ${t1}$, intent(in) :: dv, ev
+ !! SymTridiagonal matrix elements.
+ integer(ilp), intent(in) :: n
+ !! Matrix dimension.
+ type(symtridiagonal_${s1}$_type) :: A
+ !! Corresponding SymTridiagonal matrix.
+ end function 
+
+ module function initialize_symtridiagonal_impure_${s1}$(dv, ev, err) result(A)
+ !! Construct a `symtridiagonal` matrix from the rank-1 arrays
+ !! `dl`, `dv` and `du`.
+ ${t1}$, intent(in) :: dv(:), ev(:)
+ !! Tridiagonal matrix elements.
+ type(linalg_state_type), intent(out) :: err
+ !! Error handling.
+ type(symtridiagonal_${s1}$_type) :: A
+ !! Corresponding SymTridiagonal matrix.
+ end function
+
+ module function initialize_constant_symtridiagonal_impure_${s1}$(dv, ev, n, err) result(A)
+ !! Construct a `symtridiagonal` matrix with constant elements.
+ ${t1}$, intent(in) :: dv, ev
+ !! Tridiagonal matrix elements.
+ integer(ilp), intent(in) :: n
+ !! Matrix dimension.
+ type(linalg_state_type), intent(out) :: err
+ !! Error handling.
+ type(symtridiagonal_${s1}$_type) :: A
+ !! Corresponding SymTridiagonal matrix.
+ end function 
+ #:endfor
+ end interface
+
  !----------------------------------
  !----- -----
  !----- LINEAR ALGEBRA -----

src/stdlib_specialmatrices_symtridiagonal.fypp

@@ -0,0 +1,277 @@
+#: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
+
+ character(len=*), parameter :: this = "tridiagonal matrices"
+ contains
+
+ !--------------------------------
+ !----- -----
+ !----- CONSTRUCTORS -----
+ !----- -----
+ !--------------------------------
+
+ #:for k1, t1, s1 in (KINDS_TYPES)
+ 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, 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
+ ! Matrix elements.
+ A%dl = dl ; A%dv = dv ; A%du = du
+ 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.
+
+ ! 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)]
+ A%du = [(du, i = 1, n-1)]
+ end function
+ #:endfor
+
+ !-----------------------------------------
+ !----- -----
+ !----- MATRIX-VECTOR PRODUCT -----
+ !----- -----
+ !-----------------------------------------
+
+ !! spmv_tridiag
+ #:for k1, t1, s1 in (KINDS_TYPES)
+ #:for rank in RANKS
+ module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op)
+ type(tridiagonal_${s1}$_type), intent(in) :: A
+ ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
+ ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
+ real(${k1}$), intent(in), optional :: alpha
+ real(${k1}$), intent(in), optional :: beta
+ character(1), intent(in), optional :: op
+
+ ! Internal variables.
+ real(${k1}$) :: alpha_, beta_
+ integer(ilp) :: n, nrhs, ldx, ldy
+ character(1) :: op_
+ #:if rank == 1
+ ${t1}$, pointer :: xmat(:, :), ymat(:, :)
+ #:endif
+
+ ! Deal with optional arguments.
+ alpha_ = 1.0_${k1}$ ; if (present(alpha)) alpha_ = alpha
+ beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta
+ op_ = "N" ; if (present(op)) op_ = op
+
+ ! Prepare Lapack arguments.
+ n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$
+ nrhs = #{if rank==1}# 1 #{else}# size(x, dim=2, kind=ilp) #{endif}#
+
+ #:if rank == 1
+ ! Pointer trick.
+ xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y
+ call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
+ #:else
+ call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
+ #:endif
+ end subroutine
+ #:endfor
+ #:endfor
+
+ !-------------------------------------
+ !----- -----
+ !----- UTILITY FUNCTIONS -----
+ !----- -----
+ !-------------------------------------
+
+ #:for k1, t1, s1 in (KINDS_TYPES)
+ pure module function tridiagonal_to_dense_${s1}$(A) result(B)
+ !! Convert a `tridiagonal` matrix to its dense representation.
+ type(tridiagonal_${s1}$_type), intent(in) :: A
+ !! Input tridiagonal matrix.
+ ${t1}$, allocatable :: B(:, :)
+ !! Corresponding dense matrix.
+
+ ! Internal variables.
+ integer(ilp) :: i
+
+ associate (n => A%n)
+ #: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)
+ 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 = 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
+ !! Input matrix.
+ type(tridiagonal_${s1}$_type) :: B
+ #:if t1.startswith("complex")
+ B = tridiagonal(conjg(A%du), conjg(A%dv), conjg(A%dl))
+ #:else
+ B = tridiagonal(A%du, A%dv, A%dl)
+ #: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(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
+ ${t1}$, intent(in) :: alpha
+ 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, 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, 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
+
+end submodule