Merge pull request #1091 from Mahmood-Sinan/symtridiag

Add symmetric tridiagonal matrices module

Author: jvdp1

doc/specs/stdlib_specialmatrices.md

@@ -12,7 +12,7 @@ The `stdlib_specialmatrices` module provides derived types and specialized drive
 These include:
 
 - Tridiagonal matrices
-- Symmetric Tridiagonal matrices (not yet supported)
+- Symmetric tridiagonal matrices
 - Circulant matrices (not yet supported)
 - Toeplitz matrices (not yet supported)
 - Hankel matrices (not yet supported)
@@ -32,7 +32,7 @@ Experimental
 
 #### Description
 
-Tridiagonal matrices are ubiquituous in scientific computing and often appear when discretizing 1D differential operators.
+Tridiagonal matrices are ubiquitous in scientific computing and often appear when discretizing 1D differential operators.
 A generic tridiagonal matrix has the following structure:
 $$
     A
@@ -64,18 +64,127 @@ Tridiagonal matrices are available with all supported data types as `tridiagonal
 
 - 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 [, err])`
 
-- To construct a tridiagonal matrix of size `n x n` with constant diagonal elements `dl`, `dv`, and `du`:
+- To construct a tridiagonal matrix of size `n x n` with scalar diagonal elements `dl` (lower diagonal), `dv` (main diagonal), and `du` (upper diagonal):
 
-`A = ` [[stdlib_specialmatrices(module):tridiagonal(interface)]] `(dl, dv, du, n)`
+`A = ` [[stdlib_specialmatrices(module):tridiagonal(interface)]] `(dl, dv, du, n [, err])`
+
+#### Arguments
+
+##### Constructing from arrays
+
+`A = tridiagonal(dl, dv, du [, err])`
+
+- `dl`: Shall be a rank-1 array of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `dv`: Shall be a rank-1 array of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `du`: Shall be a rank-1 array of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `err` (optional): Shall be a `type(linalg_state_type)` variable. It is an `intent(out)` argument.
+If the input arrays are not of compatible sizes, the routine calls `linalg_error_handling`.
+If `err` is present, it is assigned `LINALG_VALUE_ERROR`.
+If `err` is not present, the program terminates via `error stop`.
+In either case, the elements of the resulting matrix `A` are unassigned.
+
+##### Constructing from constant diagonal values
+
+`A = tridiagonal(dl, dv, du, n [, err])`
+
+- `dl`: Shall be a scalar of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `dv`: Shall be a scalar of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `du`: Shall be a scalar of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `n`: Shall be a positive `integer`. It is an `intent(in)` argument.
+- `err` (optional): Shall be a `type(linalg_state_type)` variable. It is an `intent(out)` argument.
+If a non-positive matrix size is provided, the routine calls `linalg_error_handling`.
+If `err` is present, it is assigned `LINALG_VALUE_ERROR`.
+If `err` is not present, the program terminates via `error stop`.
+In either case, the elements of the resulting matrix `A` are unassigned.
 
 #### Example
 
 ```fortran
 {!example/specialmatrices/example_tridiagonal_dp_type.f90!}
 ```
 
+### Symmetric tridiagonal matrices {#Sym_tridiagonal}
+
+#### Status
+
+Experimental
+
+#### Description
+
+Symmetric tridiagonal matrices are a special case of tridiagonal matrices in which the subdiagonal and superdiagonal are identical.
+A generic symmetric tridiagonal matrix has the following structure:
+$$
+    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}.
+$$
+Hence, only one vector of size `n` and one vector 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`.
+Symmetric tridiagonal matrices are available with all supported data types as `sym_tridiagonal_<kind>_type`, for example:
+
+- `sym_tridiagonal_sp_type`   : Symmetric tridiagonal matrix of size `n` with `real`/`single precision` data.
+- `sym_tridiagonal_dp_type`   : Symmetric tridiagonal matrix of size `n` with `real`/`double precision` data.
+- `sym_tridiagonal_xdp_type`  : Symmetric tridiagonal matrix of size `n` with `real`/`extended precision` data.
+- `sym_tridiagonal_qp_type`   : Symmetric tridiagonal matrix of size `n` with `real`/`quadruple precision` data.
+- `sym_tridiagonal_csp_type`  : Symmetric tridiagonal matrix of size `n` with `complex`/`single precision` data.
+- `sym_tridiagonal_cdp_type`  : Symmetric tridiagonal matrix of size `n` with `complex`/`double precision` data.
+- `sym_tridiagonal_cxdp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`extended precision` data.
+- `sym_tridiagonal_cqp_type`  : Symmetric tridiagonal matrix of size `n` with `complex`/`quadruple precision` data.
+
+
+#### Syntax
+
+- To construct a symmetric tridiagonal matrix from already allocated arrays `du` (off-diagonal, size `n-1`) and `dv` (main diagonal, size `n`):
+
+`A = ` [[stdlib_specialmatrices(module):sym_tridiagonal(interface)]] `(du, dv [, err])`
+
+- To construct a symmetric tridiagonal matrix of size `n x n` with scalar diagonal elements `du` (off-diagonal) and `dv` (main diagonal):
+
+`A = ` [[stdlib_specialmatrices(module):sym_tridiagonal(interface)]] `(du, dv, n [, err])`
+
+#### Arguments
+
+##### Constructing from arrays
+
+`A = sym_tridiagonal(du, dv [, err])`
+
+- `du`: Shall be a rank-1 array of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `dv`: Shall be a rank-1 array of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `err` (optional): Shall be a `type(linalg_state_type)` variable. It is an `intent(out)` argument.
+If the input arrays are not of compatible sizes, the routine calls `linalg_error_handling`.
+If `err` is present, it is assigned `LINALG_VALUE_ERROR`.
+If `err` is not present, the program terminates via `error stop`.
+In either case, the elements of the resulting matrix `A` are unassigned.
+
+##### Constructing from constant diagonal values
+
+`A = sym_tridiagonal(du, dv, n [, err])`
+
+- `du`: Shall be a scalar of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `dv`: Shall be a scalar of type `real` or `complex` with the same kind as the resulting matrix. It is an `intent(in)` argument.
+- `n`: Shall be a positive `integer`. It is an `intent(in)` argument.
+- `err` (optional): Shall be a `type(linalg_state_type)` variable. It is an `intent(out)` argument.
+If a non-positive matrix size is provided, the routine calls `linalg_error_handling`.
+If `err` is present, it is assigned `LINALG_VALUE_ERROR`.
+If `err` is not present, the program terminates via `error stop`.
+In either case, the elements of the resulting matrix `A` are unassigned.
+
+#### Example
+
+```fortran
+{!example/specialmatrices/example_sym_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.
@@ -90,7 +199,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 backend is either LAPACK `lagtm` or the generalized routine `glagtm`, depending on the values and types of `alpha` and `beta`.
+- For `tridiagonal` and `symmetric tridiagonal` matrices, the backend is either LAPACK `lagtm` or the generalized routine `glagtm`, depending on the values and types of `alpha` and `beta`.
 
 #### Syntax
 
@@ -115,6 +224,9 @@ With the exception of `extended precision` and `quadruple precision`, all the ty
 ```fortran
 {!example/specialmatrices/example_specialmatrices_dp_spmv.f90!}
 ```
+```fortran
+{!example/specialmatrices/example_specialmatrices_cdp_spmv.f90!}
+```
 
 ## Utility functions
 
@@ -186,7 +298,7 @@ Experimental
 
 #### Description
 
-The definition of all standard artihmetic operators have been overloaded to be applicable for the matrix types defined by `stdlib_specialmatrices`:
+The definition of all standard arithmetic 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.

example/specialmatrices/CMakeLists.txt

@@ -1,3 +1,4 @@
 ADD_EXAMPLE(specialmatrices_dp_spmv)
 ADD_EXAMPLE(specialmatrices_cdp_spmv)
 ADD_EXAMPLE(tridiagonal_dp_type)
+ADD_EXAMPLE(sym_tridiagonal_dp_type)

example/specialmatrices/example_sym_tridiagonal_dp_type.f90

@@ -0,0 +1,17 @@
+program example_sym_tridiagonal_matrix
+   use stdlib_linalg_constants, only: dp
+   use stdlib_specialmatrices, only: sym_tridiagonal_dp_type, sym_tridiagonal
+   implicit none
+
+   integer, parameter :: n = 5
+   type(sym_tridiagonal_dp_type) :: A
+   real(dp) :: du(n - 1), dv(n)
+
+   ! Generate random symmteric tridiagonal elements.
+   call random_number(du)
+   call random_number(dv)
+
+   ! Create the corresponding symmetric tridiagonal matrix.
+   A = sym_tridiagonal(du, dv)
+
+end program example_sym_tridiagonal_matrix

example/specialmatrices/example_tridiagonal_dp_type.f90

@@ -1,6 +1,6 @@
 program example_tridiagonal_matrix
    use stdlib_linalg_constants, only: dp
-   use stdlib_specialmatrices
+   use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal
    implicit none
 
    integer, parameter :: n = 5

src/specialmatrices/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(specialmatrices_fppFiles
     stdlib_specialmatrices.fypp
     stdlib_specialmatrices_tridiagonal.fypp
+    stdlib_specialmatrices_sym_tridiagonal.fypp
     )
 
 set(specialmatrices_cppFiles