Skip to content

SymTridiagonal matrices and associated spmv kernel #1021

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 24 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
df6483a
Base type and constructor interfaces for symtridiagonal matrices.
loiseaujc Jul 21, 2025
20683cf
Basic implementations
loiseaujc Jul 22, 2025
f421812
Export symtridiagonal
loiseaujc Jul 22, 2025
a42928f
Added basic unit tests for symtridiagonal.
loiseaujc Jul 22, 2025
f6a2c88
Fix conjugate-transpose spmv.
loiseaujc Jul 22, 2025
9a37e13
Turned-on complex-valued tests.
loiseaujc Jul 22, 2025
d64ebef
SymTridiagonal is extended from Tridiagonal instead of being its own …
loiseaujc Jul 22, 2025
41d4a03
Regrouped everything in a single submodule.
loiseaujc Jul 22, 2025
089b297
Fixed CMakeLists.txt
loiseaujc Jul 22, 2025
fbb3474
Implementation of the Hermitian Tridiagonal matrices.
loiseaujc Jul 22, 2025
58c4eaf
Added tests for Hermitian Tridiagonal matrices.
loiseaujc Jul 22, 2025
1644315
Complete algebra (+, -, scalar mult.) for combinations of tridiagonal…
loiseaujc Jul 22, 2025
5221f13
Fix scalar multiplication for Hermitian matrices.
loiseaujc Jul 22, 2025
f26a17d
Fix missing module keyword.
loiseaujc Jul 22, 2025
54113b4
Remove tridiag_arithmetic entirely to pinpoint intel error.
loiseaujc Jul 22, 2025
4d6e526
Remove tridiagonal error handling to pinpoint issue with intel 2024.1
loiseaujc Jul 22, 2025
d1d1845
Remove all tridiag arithmetic tests to pinpoint intel issue.
loiseaujc Jul 22, 2025
6209b80
Remove symtridiag spmv kernel tests to pinpoint intel error.
loiseaujc Jul 22, 2025
ca9bf94
Remove hermitian spmv kernel test to pinpoint intel issue.
loiseaujc Jul 22, 2025
d8bcc82
Intel fix (#5)
loiseaujc Jul 23, 2025
7ea8d29
In-code documentation.
loiseaujc Jul 23, 2025
345f6f9
Added documentation for symmetric and hermitian matrices.
loiseaujc Jul 23, 2025
15b0769
Added example for symtridiagonal and hermtridiagonal.
loiseaujc Jul 23, 2025
fb589f9
Added examples to CMakeLists.txt
loiseaujc Jul 23, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
116 changes: 109 additions & 7 deletions doc/specs/stdlib_specialmatrices.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@ The `stdlib_specialmatrices` module provides derived types and specialized drive
These include:

- Tridiagonal matrices
- Symmetric Tridiagonal matrices (not yet supported)
- Symmetric Tridiagonal matrices
- Hermitian Tridiagonal matrices
- Circulant matrices (not yet supported)
- Toeplitz matrices (not yet supported)
- Hankel matrices (not yet supported)
Expand Down Expand Up @@ -76,6 +77,107 @@ Tridiagonal matrices are available with all supported data types as `tridiagonal
{!example/specialmatrices/example_tridiagonal_dp_type.f90!}
```

### Symmetric tridiagonal matrices {#SymTridiagonal}

#### Status

Experimental

#### Description

Symmetric tridiagonal matrices are ubiquituous in scientific computing and often appear when discretizing 1D differential operators.
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 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`.
Symmetric tridiagonal matrices are available with all supported data types as `symtridiagonal_<kind>_type`, for example:

- `symtridiagonal_sp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`single precision` data.
- `symtridiagonal_dp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`double precision` data.
- `symtridiagonal_xdp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`extended precision` data.
- `symtridiagonal_qp_type` : Symmetric tridiagonal matrix of size `n` with `real`/`quadruple precision` data.
- `symtridiagonal_csp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`single precision` data.
- `symtridiagonal_cdp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`double precision` data.
- `symtridiagonal_cxdp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`extended precision` data.
- `symtridiagonal_cqp_type` : Symmetric tridiagonal matrix of size `n` with `complex`/`quadruple precision` data.

#### Syntax

- To construct a symmetric tridiagonal matrix from already allocated arrays `dv` (main diagonal, size `n`, only its real part is being referenced) and `ev` (upper diagonal, size `n-1`):

`A = ` [[stdlib_specialmatrices(module):symtridiagonal(interface)]] `(dv, ev)`

- To construct a symmetric tridiagonal matrix of size `n x n` with constant diagonal elements `dv`, and `ev`:

`A = ` [[stdlib_specialmatrices(module):symtridiagonal(interface)]] `(dv, ev, n)`

#### Example

```fortran
{!example/specialmatrices/example_symtridiagonal_dp_type.f90!}
```

### Hermitian tridiagonal matrices {#HermTridiagonal}

#### Status

Experimental

#### Description

Hermitian tridiagonal matrices are ubiquituous in scientific computing.
A generic hermitian tridiagonal matrix has the following structure:
$$
A
=
\begin{bmatrix}
a_1 & b_1 \\
\bar{b}_1 & a_2 & b_2 \\
& \ddots & \ddots & \ddots \\
& & \bar{b}_{n-2} & a_{n-1} & b_{n-1} \\
& & & \bar{b}_{n-1} & a_n
\end{bmatrix},
$$
where \( a_i \in \mathbb{R} \), \( b_i \in \mathbb{C} \) and the overbar denotes the complex conjugate operation.
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`.
Hermitian tridiagonal matrices are available with all supported data types as `hermtridiagonal_<kind>_type`, for example:

- `hermtridiagonal_csp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`single precision` data.
- `hermtridiagonal_cdp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`double precision` data.
- `hermtridiagonal_cxdp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`extended precision` data.
- `hermtridiagonal_cqp_type` : Hermitian tridiagonal matrix of size `n` with `complex`/`quadruple precision` data.

#### Syntax

- To construct a hermitian tridiagonal matrix from already allocated arrays `dv` (main diagonal, size `n`) and `ev` (upper diagonal, size `n-1`):

`A = ` [[stdlib_specialmatrices(module):hermtridiagonal(interface)]] `(dv, ev)`

- To construct a hermitian tridiagonal matrix of size `n x n` with constant diagonal elements `dv`, and `ev`:

`A = ` [[stdlib_specialmatrices(module):hermtridiagonal(interface)]] `(dv, ev, n)`

Note that only the real parts of the diagonal elements `dv` are being used to construct the corresponding Hermitian matrix.

#### Example

```fortran
{!example/specialmatrices/example_hermtridiagonal_cdp_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.
Expand Down Expand Up @@ -111,7 +213,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 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.
Due to limitations of the underlying `lapack` driver, currently `alpha` and `beta` can only take one of the values `[-1, 0, 1]` for `tridiagonal`, `symtridiagonal` and `hermtridiagonal` matrices. See `lagtm` for more details.
@endwarning

#### Examples
Expand Down Expand Up @@ -192,17 +294,17 @@ Experimental

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 `+` operator for adding two matrices of the same class and kind.
- Overloading the `-` operator for subtracting two matrices of the same class and kind.
- Overloading the `*` for scalar-matrix multiplication.

#### Syntax

- Adding two matrices of the same type:
- Adding two matrices of the same class:

`C = A` [[stdlib_specialmatrices(module):operator(+)(interface)]] `B`

- Subtracting two matrices of the same type:
- Subtracting two matrices of the same class:

`C = A` [[stdlib_specialmatrices(module):operator(-)(interface)]] `B`

Expand All @@ -211,5 +313,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 (`-`), 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`, and `B` need to be of the same class and kind. For scalar multiplication (`*`), `A` and `B` need to be of the same class and kind, while `alpha` is either `real` or `complex` (with the same kind again) depending on the type being used.
@endnote
Loading
Loading