SymTridiagonal matrices and associated spmv kernel

[loiseaujc]

This PR is the second of a series to port what is available from SpecialMatrices to stdlib. It provides the base type and spmv kernel for SymTridiagonal matrices.

Proposed interfaces

• A = symtridiagonal(dv, ev [, is_posdef]) where dv and ev are rank-1 arrays describing the tridiagonal elements, and is_posdef is an optional logical flag letting the user specify that A is positive-definite (useful for eigenvalue computation and linear solver for instance).
• call spmv(A, x, y [, alpha] [, beta] [, op]) for the matrix-vector product.
• dense(A) to convert a symtridiagonal matrix to its standard rank-2 array representation.
• transpose(A) to compute the transpose of a symtridiagonal matrix.
• hermitian(A) to compute the hermitian of a symtridiagonal matrix.
• + and - are being overloaded for addition and subtraction to two symtridiagonal matrices.
• * is being overloaded for scalar-matrix multiplication.

Key facts

• The spmv kernel uses lagtm LAPACK backend under the hood.
• The kind generality of the spmv kernel has the same limitations as stdlib_linalg_lapack.

Progress

• Base implementation
• Tests
• Documentation
• In-code documentation
• Examples

cc @jvdp1, @jalvesz, @perazz

[loiseaujc]

I do have a couple of design questions, mostly regarding the complex-valued cases.

Constructors

At the moment, the complex types are defined as

type, public :: symtridiagonal_cdp
    private
    integer(ilp) :: n
    complex(dp), allocatable :: dv(:), ev(:)
    logical(lk) :: is_posdef
end type

where dv contains the diagonal elements, and ev the ones on the first super/sub-diagonal. A complex-valued matrix necessarily has real elements on its diagonal. Hence, should we define dv as a real-valued array and then make a complex copy wherever needed when calling some complex lapack routines or leave it as it is and add some test in the constructor to catch a dv with non-zero imaginary part ?

Symmetric vs Hermitian matrices

At the moment, symmetric matrix types have been defined for real and complex-valued arithmetic. Yet, I've hardly seen symmetric complex-valued matrices in the wild. Hermitian matrices on the other hand are a big thing. Should we restrict the symtridiagonal to real-valued matrices only and then have a hermtridiagonal for the complex ones, or have both symtridiagonal and hermtridiagonal for complex-valued matrices (with the latter being hardly used in practice I suppose) ?

Own type vs extending from tridiagonal

At the moment, the symmetric matrices have their own types. This was done mainly out of simplicity (copy-pasting-replacing some bits and pieces). We could however define them as an extension of the tridiagonal types and avoid some code duplication. I don't have a strong opinion about this (avoiding code duplication is good, but having a self-contained submodule is IMO good also). What do you think ? Will go the type extension route actually, feels better.

[loiseaujc]

I think it is pretty much ready for review, at least code wise.

Addition to stdlib_specialmatrices

• Type symtridiagonal extended from tridiagonal

  • Both real and complex are being supported, as well as all precisions.
  • Constructors
  • Matrix-vector product (spmv)
  • Scalar matrix multiplication
  • Utilities: dense, transpose, hermitian
  • Operator overloading: + and -

• Type hermtridiagonal extended from tridiagonal

  • Only complex is being supported, with all precisions.
  • Constructors
  • Matrix-vector product (spmv)
  • Scalar matrix multiplication
  • Utilities: dense, transpose, hermitian
  • Operator overloading: + and -

Note that the operator overloading handles all the possible cases (e.g. tridiagonal + tridiagonal = tridiagonal, symtridiagonal + tridiagonal = tridiagonal, symtridiagonal + symtridiagonal = symtridiagonal, etc). The logic for the is_posdef flag is not yet implemented as I cannot test it thoroughly yet (it needs the extension of the eigh interface which will come in a following PR).

What is missing? (July 22nd)

• Documentation
• Examples

Current issues

It seems like there is a problem with the intel compiler, specifically the 2024.1 version. Not sure what's going on as I do not have it installed on my computer at the moment. If any of you has it and want to give it a shot, be my guest. Otherwise, I'll try to install in the upcoming week as I'll be on vacation with no kids.