Tridagonal matrices and associated spmv kernel (fortran-lang#957)

* Added the Tridiagonal type and associated spmv kernel.

* Rename dense to dense_mat because of a naming conflict.

* Added test for Tridiagonal spmv.

* Changed comparison to account for floating point errors.

* 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.

* Implementation of basic operations for Tridiagonal matrices is done.

* In-code documentation.

* Added examples.

* Specifications page.

* Fix typos in module header.

* Fix errors in examples.

* Enabled `xdp` and `qp` for `spmv` kernels.

* Update test/linalg/test_linalg_sparse.fypp

Remove trailing white space.

Co-authored-by: Jeremie Vandenplas <[email protected]>

* Update example/specialmatrices/example_specialmatrices_dp_spmv.f90

Import only required elements from the module.

Co-authored-by: Jeremie Vandenplas <[email protected]>

* Consistent lower-casing and imports.

* Fixed allocation from source for complex types.

* Make use of stdlib_constants to define the various zeros and ones.

* Prevent temporary arrays + lower casing every names.

* Import everything from stdlib_kinds

* Fix imports.

* Fix missing argument in check.

* Update doc/specs/stdlib_specialmatrices.md

Co-authored-by: Federico Perini <[email protected]>

* Update doc/specs/stdlib_specialmatrices.md

Co-authored-by: Federico Perini <[email protected]>

* Update doc/specs/stdlib_specialmatrices.md

Co-authored-by: Federico Perini <[email protected]>

* Update src/stdlib_specialmatrices.fypp

Co-authored-by: Federico Perini <[email protected]>

* Update src/stdlib_specialmatrices_tridiagonal.fypp

Co-authored-by: Federico Perini <[email protected]>

* Improved error handling.

* Added test for tridiagonal error handling.

---------

Co-authored-by: Jeremie Vandenplas <[email protected]>
Co-authored-by: Federico Perini <[email protected]>

Author: 3 people

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`.
+Tridiagonal matrices are available with all supported data types as `tridiagonal_<kind>_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.
+- `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 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
+
+```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 (`-`), 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

example/CMakeLists.txt

@@ -25,6 +25,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)

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
+   ! 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

example/specialmatrices/CMakeLists.txt

@@ -0,0 +1,2 @@
+ADD_EXAMPLE(specialmatrices_dp_spmv)
+ADD_EXAMPLE(tridiagonal_dp_type)

example/specialmatrices/example_specialmatrices_dp_spmv.f90

@@ -0,0 +1,26 @@
+program example_tridiagonal_matrix
+   use stdlib_linalg_constants, only: dp
+   use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal, dense, spmv
+   implicit none
+
+   integer, parameter :: n = 5
+   type(tridiagonal_dp_type) :: A
+   real(dp) :: dl(n - 1), dv(n), du(n - 1)
+   real(dp) :: x(n), y(n), y_dense(n)
+   integer :: i
+
+   ! Create an arbitrary tridiagonal matrix.
+   dl = [(i, i=1, n - 1)]; dv = [(2*i, i=1, n)]; du = [(3*i, i=1, n - 1)]
+   A = tridiagonal(dl, dv, du)
+
+   ! Initialize vectors.
+   x = 1.0_dp; y = 0.0_dp; y_dense = 0.0_dp
+
+   ! Perform matrix-vector products.
+   call spmv(A, x, y)
+   y_dense = matmul(dense(A), x)
+
+   print *, 'dense       :', y_dense
+   print *, 'Tridiagonal :', y
+
+end program example_tridiagonal_matrix

example/specialmatrices/example_tridiagonal_dp_type.f90

@@ -0,0 +1,18 @@
+program example_tridiagonal_matrix
+   use stdlib_linalg_constants, only: dp
+   use stdlib_specialmatrices
+   implicit none
+
+   integer, parameter :: n = 5
+   type(Tridiagonal_dp_type) :: A
+   real(dp) :: dl(n - 1), dv(n), du(n - 1)
+
+   ! Generate random tridiagonal elements.
+   call random_number(dl)
+   call random_number(dv)
+   call random_number(du)
+
+   ! Create the corresponding Tridiagonal matrix.
+   A = Tridiagonal(dl, dv, du)
+
+end program example_tridiagonal_matrix

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
@@ -55,6 +55,8 @@ set(fppFiles
     stdlib_specialfunctions_activations.fypp
     stdlib_specialfunctions_gamma.fypp
     stdlib_specialfunctions.fypp
+    stdlib_specialmatrices.fypp
+    stdlib_specialmatrices_tridiagonal.fypp    
     stdlib_stats.fypp
     stdlib_stats_corr.fypp
     stdlib_stats_cov.fypp