Skip to content

Add symmetric tridiagonal matrices module#1091

Open
Mahmood-Sinan wants to merge 23 commits intofortran-lang:masterfrom
Mahmood-Sinan:symtridiag
Open

Add symmetric tridiagonal matrices module#1091
Mahmood-Sinan wants to merge 23 commits intofortran-lang:masterfrom
Mahmood-Sinan:symtridiag

Conversation

@Mahmood-Sinan
Copy link
Contributor

@Mahmood-Sinan Mahmood-Sinan commented Jan 21, 2026

This draft PR adds support for symmetric tridiagonal matrices inside specialmatrices in stdlib.
It includes implementation, tests, documentation and examples.
The PR implements the same arithmetic and matrix operations as the existing tridiagonal module.

Remaining tasks:

Feedback is welcome.

@codecov
Copy link

codecov bot commented Jan 21, 2026

Codecov Report

❌ Patch coverage is 0% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 68.50%. Comparing base (843ffb0) to head (296c35a).
⚠️ Report is 11 commits behind head on master.

Files with missing lines Patch % Lines
...pecialmatrices/example_sym_tridiagonal_dp_type.f90 0.00% 7 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1091      +/-   ##
==========================================
- Coverage   68.55%   68.50%   -0.06%     
==========================================
  Files         396      397       +1     
  Lines       12746    12753       +7     
  Branches     1376     1376              
==========================================
- Hits         8738     8736       -2     
- Misses       4008     4017       +9     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@loiseaujc
Copy link
Contributor

It looks pretty nice. One little thing though (even if I overlooked it when implementing the tridiagonal_type ) : you could define a single pure subroutine to construct the symtridiagonal matrices (make it two if you include the constant diagonal cases) and then have it called by the pure and impure constructors. It's a little thing but that would reduce a bit code duplication I guess.

I'll go through the code a bit more thoroughly tomorrow but it looks pretty good to me already.

@Mahmood-Sinan
Copy link
Contributor Author

One little thing though (even if I overlooked it when implementing the tridiagonal_type ) : you could define a single pure subroutine to construct the symtridiagonal matrices (make it two if you include the constant diagonal cases) and then have it called by the pure and impure constructors. It's a little thing but that would reduce a bit code duplication I guess.

Thanks for the suggestion! You’re right, I implemented so. Now the tridiagonal and symtridiagonal construction is now handled by two pure helper routines(one from arrays and one from constants), which is called by both the pure and impure constructors.

Copy link
Contributor

@loiseaujc loiseaujc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll take a deeper look at it by the end of the week, but it looks quite good to me so far.

implicit none
private
public :: tridiagonal
public :: tridiagonal, sym_tridiagonal
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is mostly a question of personal tastes, but I prefer symtridiagonal over sym_tridiagonal. The no underscore naming is also the one used by Julia. What do you guys think @perazz, @jvdp1 and @jalvesz ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Modulo this very style question, I think it's pretty much all good for me !

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also prefer symtridiagonal. However, the style guide says that "Variable and procedure names should be made up of one or more full words separated by an underscore,"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As @jvdp1 mentioned, if I remember correctly I chose sym_tridiagonal for the same reason

@loiseaujc loiseaujc marked this pull request as ready for review February 5, 2026 09:52
Copy link
Contributor

@loiseaujc loiseaujc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM.

@loiseaujc
Copy link
Contributor

@perazz, @jalvesz and @jvdp1 : This PR is ready for merging. The only question I have left is a triviality though. What convention should we use for naming types : sym_tridiagonal_dp_type or symtridiagonal_dp_type ? My preference goes to the first one personally.

@Mahmood-Sinan
Copy link
Contributor Author

@perazz, @jalvesz and @jvdp1 : This PR is ready for merging. The only question I have left is a triviality though. What convention should we use for naming types : sym_tridiagonal_dp_type or symtridiagonal_dp_type ? My preference goes to the first one personally.

@perazz @jalvesz @jvdp1 Please look into this.

Copy link
Member

@jvdp1 jvdp1 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @Mahmood-Sinan for this PR. Please find below some comments and suggestions.

\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.
Interfaces to the most common ones will soon be provided by `stdlib_specialmatrices`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you mean be "soon"? Is it not what is proposed in this PR?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for noticing. I have missed out a sentence before this which was common for both tridiagonal and symmetric tridiagonal. I added that sentence now. The interfances that this sentence talks about were solve, inverse, etc...


`A = ` [[stdlib_specialmatrices(module):sym_tridiagonal(interface)]] `(du, dv)`

- To construct a symmetric tridiagonal matrix of size `n x n` with constant diagonal elements `du` and `dv`:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

by "constant" do you mean "scalar"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I changed it so scalar now.

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

`A = ` [[stdlib_specialmatrices(module):sym_tridiagonal(interface)]] `(du, dv, n)`

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please define (type, kind, intent, ...) the variables du, dv, n.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have added them now under a new heading named arguments.

@@ -0,0 +1,17 @@
program example_sym_tridiagonal_matrix
use stdlib_linalg_constants, only: dp
use stdlib_specialmatrices
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
use stdlib_specialmatrices
use stdlib_specialmatrices, only: sym_tridiagonal_dp_type, sym_tridiagonal

implicit none
private
public :: tridiagonal
public :: tridiagonal, sym_tridiagonal
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also prefer symtridiagonal. However, the style guide says that "Variable and procedure names should be made up of one or more full words separated by an underscore,"

${t1}$, intent(in) :: alpha
type(sym_tridiagonal_${s1}$_type) :: B
B = sym_tridiagonal(A%du, A%dv)
B%du = alpha*B%du; B%dv = alpha*B%dv;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
B%du = alpha*B%du; B%dv = alpha*B%dv;
B%du = alpha*B%du
B%dv = alpha*B%dv

and similarly at other places...

Comment on lines 142 to 143
call random_number(du) ; call random_number(dv)
A = sym_tridiagonal(du, dv) ; Amat = dense(A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
call random_number(du) ; call random_number(dv)
A = sym_tridiagonal(du, dv) ; Amat = dense(A)
call random_number(du)
call random_number(dv)
A = sym_tridiagonal(du, dv)
Amat = dense(A)

Comment on lines 146 to 147
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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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)
allocate(x(n))
call random_number(x)
allocate(y1(n), source = 0.0_wp)
allocate(y2(n), source=0.0_wp)

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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
y1 = matmul(Amat, x) ; call spmv(A, x, y2)
y1 = matmul(Amat, x)
call spmv(A, x, y2)


! Test y = A @ x
y1 = matmul(Amat, x) ; call spmv(A, x, y2)
call check(error, all_close(y1, y2), .true.)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment in each check to facilitate the debugging, in case a test is failing (currently it won't be easy to trace back which one is failing)

@Mahmood-Sinan
Copy link
Contributor Author

Thank you @Mahmood-Sinan for this PR. Please find below some comments and suggestions.

Thank you so much for reviewing. I have made those changes. Please take a look at them.

@Mahmood-Sinan
Copy link
Contributor Author

@jvdp1 Please handle.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants