#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
#:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES
#:set RHS_SUFFIX = ["one","many"]
#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]]
#:set RHS_EMPTY = [emptyranksuffix(r) for r in [1,2]]
#:set ALL_RHS = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY))
module stdlib_linalg
!!Provides a support for various linear algebra procedures
!! ([Specification](../page/specs/stdlib_linalg.html))
use stdlib_kinds, only: xdp, int8, int16, int32, int64
use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp
use stdlib_error, only: error_stop
use stdlib_optval, only: optval
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
implicit none
private
public :: chol
public :: cholesky
public :: det
public :: operator(.det.)
public :: diag
public :: eig
public :: eigh
public :: eigvals
public :: eigvalsh
public :: eye
public :: inv
public :: invert
public :: operator(.inv.)
public :: pinv
public :: pseudoinvert
public :: operator(.pinv.)
public :: lstsq
public :: lstsq_space
public :: norm
public :: mnorm
public :: get_norm
public :: solve
public :: solve_lu
public :: solve_lstsq
public :: trace
public :: svd
public :: svdvals
public :: outer_product
public :: kronecker_product
public :: cross_product
public :: qr
public :: qr_space
public :: is_square
public :: is_diagonal
public :: is_symmetric
public :: is_skew_symmetric
public :: hermitian
public :: is_hermitian
public :: is_triangular
public :: is_hessenberg
! Export linalg error handling
public :: linalg_state_type, linalg_error_handling
interface chol
!! version: experimental
!!
!! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \).
!! ([Specification](../page/specs/stdlib_linalg.html#chol-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
!!
!!### Summary
!! Pure function interface for computing the Cholesky triangular factors.
!!
!!### Description
!!
!! This interface provides methods for computing the lower- or upper- triangular matrix from the
!! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
!! Supported data types include `real` and `complex`.
!!
!!@note The solution is based on LAPACK's `*POTRF` methods.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module function stdlib_linalg_${ri}$_cholesky_fun(a,lower,other_zeroed) result(c)
!> Input matrix a[m,n]
${rt}$, intent(in) :: a(:,:)
!> [optional] is the lower or upper triangular factor required? Default = lower
logical(lk), optional, intent(in) :: lower
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
logical(lk), optional, intent(in) :: other_zeroed
!> Output matrix with Cholesky factors c[n,n]
${rt}$ :: c(size(a,1),size(a,2))
end function stdlib_linalg_${ri}$_cholesky_fun
#:endfor
end interface chol
interface cholesky
!! version: experimental
!!
!! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \).
!! ([Specification](../page/specs/stdlib_linalg.html#cholesky-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
!!
!!### Summary
!! Pure subroutine interface for computing the Cholesky triangular factors.
!!
!!### Description
!!
!! This interface provides methods for computing the lower- or upper- triangular matrix from the
!! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
!! Supported data types include `real` and `complex`.
!! The factorization is computed in-place if only one matrix argument is present; or returned into
!! a second matrix argument, if present. The `lower` `logical` flag allows to select between upper or
!! lower factorization; the `other_zeroed` optional `logical` flag allows to choose whether the unused
!! part of the triangular matrix should be filled with zeroes.
!!
!!@note The solution is based on LAPACK's `*POTRF` methods.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine stdlib_linalg_${ri}$_cholesky_inplace(a,lower,other_zeroed,err)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> [optional] is the lower or upper triangular factor required? Default = lower
logical(lk), optional, intent(in) :: lower
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
logical(lk), optional, intent(in) :: other_zeroed
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_cholesky_inplace
pure module subroutine stdlib_linalg_${ri}$_cholesky(a,c,lower,other_zeroed,err)
!> Input matrix a[n,n]
${rt}$, intent(in) :: a(:,:)
!> Output matrix with Cholesky factors c[n,n]
${rt}$, intent(out) :: c(:,:)
!> [optional] is the lower or upper triangular factor required? Default = lower
logical(lk), optional, intent(in) :: lower
!> [optional] should the unused half of the return matrix be zeroed out? Default: yes
logical(lk), optional, intent(in) :: other_zeroed
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_cholesky
#:endfor
end interface cholesky
interface diag
!! version: experimental
!!
!! Creates a diagonal array or extract the diagonal elements of an array
!! ([Specification](../page/specs/stdlib_linalg.html#
!! diag-create-a-diagonal-array-or-extract-the-diagonal-elements-of-an-array))
!
! Vector to matrix
!
#:for k1, t1 in RCI_KINDS_TYPES
module function diag_${t1[0]}$${k1}$(v) result(res)
${t1}$, intent(in) :: v(:)
${t1}$ :: res(size(v),size(v))
end function diag_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
module function diag_${t1[0]}$${k1}$_k(v,k) result(res)
${t1}$, intent(in) :: v(:)
integer, intent(in) :: k
${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
end function diag_${t1[0]}$${k1}$_k
#:endfor
!
! Matrix to vector
!
#:for k1, t1 in RCI_KINDS_TYPES
module function diag_${t1[0]}$${k1}$_mat(A) result(res)
${t1}$, intent(in) :: A(:,:)
${t1}$ :: res(minval(shape(A)))
end function diag_${t1[0]}$${k1}$_mat
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
module function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
${t1}$, intent(in) :: A(:,:)
integer, intent(in) :: k
${t1}$ :: res(minval(shape(A))-abs(k))
end function diag_${t1[0]}$${k1}$_mat_k
#:endfor
end interface
! Matrix trace
interface trace
!! version: experimental
!!
!! Computes the trace of a matrix
!! ([Specification](../page/specs/stdlib_linalg.html#
!! trace-trace-of-a-matrix))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure trace_${t1[0]}$${k1}$
#:endfor
end interface
! Identity matrix
interface eye
!! version: experimental
!!
!! Constructs the identity matrix
!! ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure eye_${t1[0]}$${k1}$
#:endfor
end interface eye
! Outer product (of two vectors)
interface outer_product
!! version: experimental
!!
!! Computes the outer product of two vectors, returning a rank-2 array
!! ([Specification](../page/specs/stdlib_linalg.html#
!! outer_product-computes-the-outer-product-of-two-vectors))
#:for k1, t1 in RCI_KINDS_TYPES
pure module function outer_product_${t1[0]}$${k1}$(u, v) result(res)
${t1}$, intent(in) :: u(:), v(:)
${t1}$ :: res(size(u),size(v))
end function outer_product_${t1[0]}$${k1}$
#:endfor
end interface outer_product
interface kronecker_product
!! version: experimental
!!
!! Computes the Kronecker product of two arrays of size M1xN1, and of M2xN2, returning an (M1*M2)x(N1*N2) array
!! ([Specification](../page/specs/stdlib_linalg.html#
!! kronecker_product-computes-the-kronecker-product-of-two-matrices))
#:for k1, t1 in RCI_KINDS_TYPES
pure module function kronecker_product_${t1[0]}$${k1}$(A, B) result(C)
${t1}$, intent(in) :: A(:,:), B(:,:)
${t1}$ :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
end function kronecker_product_${t1[0]}$${k1}$
#:endfor
end interface kronecker_product
! Cross product (of two vectors)
interface cross_product
!! version: experimental
!!
!! Computes the cross product of two vectors, returning a rank-1 and size-3 array
!! ([Specification](../page/specs/stdlib_linalg.html#cross_product-computes-the-cross-product-of-two-3-d-vectors))
#:for k1, t1 in RCI_KINDS_TYPES
pure module function cross_product_${t1[0]}$${k1}$(a, b) result(res)
${t1}$, intent(in) :: a(3), b(3)
${t1}$ :: res(3)
end function cross_product_${t1[0]}$${k1}$
#:endfor
end interface cross_product
! Check for squareness
interface is_square
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is square
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_square-checks-if-a-matrix-is-square))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_square_${t1[0]}$${k1}$
#:endfor
end interface is_square
! Check for diagonality
interface is_diagonal
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is diagonal
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_diagonal-checks-if-a-matrix-is-diagonal))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_diagonal_${t1[0]}$${k1}$
#:endfor
end interface is_diagonal
! Check for symmetry
interface is_symmetric
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is symmetric
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_symmetric-checks-if-a-matrix-is-symmetric))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_symmetric_${t1[0]}$${k1}$
#:endfor
end interface is_symmetric
! Check for skew-symmetry
interface is_skew_symmetric
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is skew-symmetric
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_skew_symmetric-checks-if-a-matrix-is-skew-symmetric))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_skew_symmetric_${t1[0]}$${k1}$
#:endfor
end interface is_skew_symmetric
! Check for Hermiticity
interface is_hermitian
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is Hermitian
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_hermitian-checks-if-a-matrix-is-hermitian))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_hermitian_${t1[0]}$${k1}$
#:endfor
end interface is_hermitian
interface hermitian
!! version: experimental
!!
!! Computes the Hermitian version of a rank-2 matrix.
!! For complex matrices, this returns `conjg(transpose(a))`.
!! For real or integer matrices, this returns `transpose(a)`.
!!
!! Usage:
!! ```
!! A = reshape([(1, 2), (3, 4), (5, 6), (7, 8)], [2, 2])
!! AH = hermitian(A)
!! ```
!!
!! [Specification](../page/specs/stdlib_linalg.html#hermitian-compute-the-hermitian-version-of-a-rank-2-matrix)
!!
#:for k1, t1 in RCI_KINDS_TYPES
pure module function hermitian_${t1[0]}$${k1}$(a) result(ah)
${t1}$, intent(in) :: a(:,:)
${t1}$ :: ah(size(a, 2), size(a, 1))
end function hermitian_${t1[0]}$${k1}$
#:endfor
end interface hermitian
! Check for triangularity
interface is_triangular
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is triangular
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_triangular-checks-if-a-matrix-is-triangular))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_triangular_${t1[0]}$${k1}$
#:endfor
end interface is_triangular
! Check for matrix being Hessenberg
interface is_hessenberg
!! version: experimental
!!
!! Checks if a matrix (rank-2 array) is Hessenberg
!! ([Specification](../page/specs/stdlib_linalg.html#
!! is_hessenberg-checks-if-a-matrix-is-hessenberg))
#:for k1, t1 in RCI_KINDS_TYPES
module procedure is_Hessenberg_${t1[0]}$${k1}$
#:endfor
end interface is_hessenberg
! Solve linear system system Ax=b.
interface solve
!! version: experimental
!!
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-equations))
!!
!!### Summary
!! Interface for solving a linear system arising from a general matrix.
!!
!!### Description
!!
!! This interface provides methods for computing the solution of a linear matrix system.
!! Supported data types include `real` and `complex`. No assumption is made on the matrix
!! structure.
!! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`)
!! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
!!
!!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x)
!> Input matrix a[n,n]
${rt}$, intent(inout), target :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Result array/matrix x[n] or x[n,nrhs]
${rt}$, allocatable, target :: x${nd}$
end function stdlib_linalg_${ri}$_solve_${ndsuf}$
pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x)
!> Input matrix a[n,n]
${rt}$, intent(in) :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> Result array/matrix x[n] or x[n,nrhs]
${rt}$, allocatable, target :: x${nd}$
end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$
#:endfor
#:endfor
end interface solve
! Solve linear system Ax = b using LU decomposition (subroutine interface).
interface solve_lu
!! version: experimental
!!
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#solve-lu-solves-a-linear-matrix-equation-or-a-linear-system-of-equations-subroutine-interface))
!!
!!### Summary
!! Subroutine interface for solving a linear system using LU decomposition.
!!
!!### Description
!!
!! This interface provides methods for computing the solution of a linear matrix system using
!! a subroutine. Supported data types include `real` and `complex`. No assumption is made on the matrix
!! structure. Preallocated space for the solution vector `x` is user-provided, and it may be provided
!! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal
!! memory allocations take place when using this interface.
!! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`)
!! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
!!
!!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)
!> Input matrix a[n,n]
${rt}$, intent(inout), target :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> Result array/matrix x[n] or x[n,nrhs]
${rt}$, intent(inout), contiguous, target :: x${nd}$
!> [optional] Storage array for the diagonal pivot indices
integer(ilp), optional, intent(inout), target :: pivot(:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
#:endfor
#:endfor
end interface solve_lu
! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
interface lstsq
!! version: experimental
!!
!! Computes the squares solution to system \( A \cdot x = b \).
!! ([Specification](../page/specs/stdlib_linalg.html#lstsq-computes-the-least-squares-solution-to-a-linear-matrix-equation))
!!
!!### Summary
!! Interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution.
!!
!!### Description
!!
!! This interface provides methods for computing the least squares of a linear matrix system.
!! Supported data types include `real` and `complex`.
!!
!!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x)
!> Input matrix a[n,n]
${rt}$, intent(inout), target :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
real(${rk}$), optional, intent(in) :: cond
!> [optional] Can A,b data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Return rank of A
integer(ilp), optional, intent(out) :: rank
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
!> Result array/matrix x[n] or x[n,nrhs]
${rt}$, allocatable, target :: x${nd}$
end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$
#:endfor
#:endfor
end interface lstsq
! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
interface solve_lstsq
!! version: experimental
!!
!! Computes the squares solution to system \( A \cdot x = b \).
!! ([Specification](../page/specs/stdlib_linalg.html#solve-lstsq-compute-the-least-squares-solution-to-a-linear-matrix-equation-subroutine-interface))
!!
!!### Summary
!! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution.
!!
!!### Description
!!
!! This interface provides methods for computing the least squares of a linear matrix system using
!! a subroutine. Supported data types include `real` and `complex`. If pre-allocated work spaces
!! are provided, no internal memory allocations take place when using this interface.
!!
!!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,&
#{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,singvals,overwrite_a,rank,err)
!> Input matrix a[n,n]
${rt}$, intent(inout), target :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> Result array/matrix x[n] or x[n,nrhs]
${rt}$, intent(inout), contiguous, target :: x${nd}$
!> [optional] real working storage space
real(${rk}$), optional, intent(inout), target :: real_storage(:)
!> [optional] integer working storage space
integer(ilp), optional, intent(inout), target :: int_storage(:)
#:if rt.startswith('complex')
!> [optional] complex working storage space
${rt}$, optional, intent(inout), target :: cmpl_storage(:)
#:endif
!> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
real(${rk}$), optional, intent(in) :: cond
!> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD
real(${rk}$), optional, intent(out), target :: singvals(:)
!> [optional] Can A,b data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Return rank of A
integer(ilp), optional, intent(out) :: rank
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$
#:endfor
#:endfor
end interface solve_lstsq
! Return the working array space required by the least squares solver
interface lstsq_space
!! version: experimental
!!
!! Computes the integer, real [, complex] working space required by the least-squares solver
!! ([Specification](../page/specs/stdlib_linalg.html#lstsq-space-compute-internal-working-space-requirements-for-the-least-squares-solver))
!!
!!### Description
!!
!! This interface provides sizes of integer, real [, complex] working spaces required by the
!! least-squares solver. These sizes can be used to pre-allocated working arrays in case several
!! repeated least-squares solutions to a same system are sought. If pre-allocated working arrays
!! are provided, no internal allocations will take place.
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Right hand side vector or array, b[n] or b[n,nrhs]
${rt}$, intent(in) :: b${nd}$
!> Size of the working space arrays
integer(ilp), intent(out) :: lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#
end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$
#:endfor
#:endfor
end interface lstsq_space
! QR factorization of rank-2 array A
interface qr
!! version: experimental
!!
!! Computes the QR factorization of matrix \( A = Q R \).
!! ([Specification](../page/specs/stdlib_linalg.html#qr-compute-the-qr-factorization-of-a-matrix))
!!
!!### Summary
!! Compute the QR factorization of a `real` or `complex` matrix: \( A = Q R \), where \( Q \) is orthonormal
!! and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m\ge n \).
!!
!!### Description
!!
!! This interface provides methods for computing the QR factorization of a matrix.
!! Supported data types include `real` and `complex`. If a pre-allocated work space
!! is provided, no internal memory allocations take place when using this interface.
!!
!! Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \).
!! The user may want the full problem (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) or the reduced
!! problem only: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`).
!!
!!@note The solution is based on LAPACK's QR factorization (`*GEQRF`) and ordered matrix output (`*ORGQR`, `*UNGQR`).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Orthogonal matrix Q ([m,m], or [m,k] if reduced)
${rt}$, intent(out), contiguous, target :: q(:,:)
!> Upper triangular matrix R ([m,n], or [k,n] if reduced)
${rt}$, intent(out), contiguous, target :: r(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Provide pre-allocated workspace, size to be checked with qr_space
${rt}$, intent(out), optional, target :: storage(:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_qr
#:endfor
end interface qr
! Return the working array space required by the QR factorization solver
interface qr_space
!! version: experimental
!!
!! Computes the working array space required by the QR factorization solver
!! ([Specification](../page/specs/stdlib_linalg.html#qr-space-compute-internal-working-space-requirements-for-the-qr-factorization))
!!
!!### Description
!!
!! This interface returns the size of the `real` or `complex` working storage required by the
!! QR factorization solver. The working size only depends on the kind (`real` or `complex`) and size of
!! the matrix being factorized. Storage size can be used to pre-allocate a working array in case several
!! repeated QR factorizations to a same-size matrix are sought. If pre-allocated working arrays
!! are provided, no internal allocations will take place during the factorization.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine get_qr_${ri}$_workspace(a,lwork,err)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Minimum workspace size for both operations
integer(ilp), intent(out) :: lwork
!> State return flag. Returns an error if the query failed
type(linalg_state_type), optional, intent(out) :: err
end subroutine get_qr_${ri}$_workspace
#:endfor
end interface qr_space
interface det
!! version: experimental
!!
!! Computes the determinant of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#det-computes-the-determinant-of-a-square-matrix))
!!
!!### Summary
!! Interface for computing matrix determinant.
!!
!!### Description
!!
!! This interface provides methods for computing the determinant of a matrix.
!! Supported data types include `real` and `complex`.
!!
!!@note The provided functions are intended for square matrices only.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
!!### Example
!!
!!```fortran
!!
!! real(sp) :: a(3,3), d
!! type(linalg_state_type) :: state
!! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!!
!! ! ...
!! d = det(a,err=state)
!! if (state%ok()) then
!! print *, 'Success! det=',d
!! else
!! print *, state%print()
!! endif
!! ! ...
!!```
!!
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module procedure stdlib_linalg_${rt[0]}$${rk}$determinant
module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface det
interface operator(.det.)
!! version: experimental
!!
!! Determinant operator of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#det-determinant-operator-of-a-square-matrix))
!!
!!### Summary
!! Pure operator interface for computing matrix determinant.
!!
!!### Description
!!
!! This pure operator interface provides a convenient way to compute the determinant of a matrix.
!! Supported data types include real and complex.
!!
!!@note The provided functions are intended for square matrices.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
!!### Example
!!
!!```fortran
!!
!! ! ...
!! real(sp) :: matrix(3,3), d
!! matrix = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!! d = .det.matrix
!! ! ...
!!
!!```
!
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface operator(.det.)
interface
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> State return flag.
type(linalg_state_type), intent(out) :: err
!> Matrix determinant
${rt}$ :: det
end function stdlib_linalg_${rt[0]}$${rk}$determinant
pure module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det)
!> Input matrix a[m,n]
${rt}$, intent(in) :: a(:,:)
!> Matrix determinant
${rt}$ :: det
end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface
! Matrix Inverse: Function interface
interface inv
!! version: experimental
!!
!! Inverse of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-of-a-square-matrix))
!!
!!### Summary
!! This interface provides methods for computing the inverse of a square `real` or `complex` matrix.
!! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
!!
!!### Description
!!
!! This function interface provides methods that return the inverse of a square matrix.
!! Supported data types include `real` and `complex`.
!! The inverse matrix \( A^{-1} \) is returned as a function result.
!! Exceptions are raised in case of singular matrix or invalid size, and trigger an `error stop` if
!! the state flag `err` is not provided.
!!
!!@note The provided functions are intended for square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_inverse_${ri}$(a,err) result(inva)
!> Input matrix a[n,n]
${rt}$, intent(in) :: a(:,:)
!> Output matrix inverse
${rt}$, allocatable :: inva(:,:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end function stdlib_linalg_inverse_${ri}$
#:endfor
end interface inv
! Matrix Inverse: Subroutine interface - in-place inversion
interface invert
!! version: experimental
!!
!! Inversion of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#invert-inversion-of-a-square-matrix))
!!
!!### Summary
!! This interface provides methods for inverting a square `real` or `complex` matrix in-place.
!! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
!!
!!### Description
!!
!! This subroutine interface provides a way to compute the inverse of a matrix.
!! Supported data types include `real` and `complex`.
!! The user may provide a unique matrix argument `a`. In this case, `a` is replaced by the inverse matrix.
!! on output. Otherwise, one may provide two separate arguments: an input matrix `a` and an output matrix
!! `inva`. In this case, `a` will not be modified, and the inverse is returned in `inva`.
!! Pre-allocated storage may be provided for the array of pivot indices, `pivot`. If all pre-allocated
!! work spaces are provided, no internal memory allocations take place when using this interface.
!!
!!@note The provided subroutines are intended for square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err)
!> Input matrix a[n,n]
${rt}$, intent(inout) :: a(:,:)
!> [optional] Storage array for the diagonal pivot indices
integer(ilp), optional, intent(inout), target :: pivot(:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_invert_inplace_${ri}$
! Compute the square matrix inverse of a
module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err)
!> Input matrix a[n,n].
${rt}$, intent(in) :: a(:,:)
!> Inverse matrix a[n,n].
${rt}$, intent(out) :: inva(:,:)
!> [optional] Storage array for the diagonal pivot indices
integer(ilp), optional, intent(inout), target :: pivot(:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_invert_split_${ri}$
#:endfor
end interface invert
! Matrix Inverse: Operator interface
interface operator(.inv.)
!! version: experimental
!!
!! Inverse operator of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-operator-of-a-square-matrix))
!!
!!### Summary
!! Operator interface for computing the inverse of a square `real` or `complex` matrix.
!!
!!### Description
!!
!! This operator interface provides a convenient way to compute the inverse of a matrix.
!! Supported data types include `real` and `complex`. On input errors or singular matrix,
!! NaNs will be returned.
!!
!!@note The provided functions are intended for square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva)
!> Input matrix a[n,n]
${rt}$, intent(in) :: a(:,:)
!> Result matrix
${rt}$, allocatable :: inva(:,:)
end function stdlib_linalg_inverse_${ri}$_operator
#:endfor
end interface operator(.inv.)
! Moore-Penrose Pseudo-Inverse: Function interface
interface pinv
!! version: experimental
!!
!! Pseudo-inverse of a matrix
!! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-of-a-matrix))
!!
!!### Summary
!! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a matrix.
!! The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse, computed for square, singular,
!! or rectangular matrices. It is defined such that it satisfies the conditions:
!! - \( A \cdot A^{+} \cdot A = A \)
!! - \( A^{+} \cdot A \cdot A^{+} = A^{+} \)
!! - \( (A \cdot A^{+})^T = A \cdot A^{+} \)
!! - \( (A^{+} \cdot A)^T = A^{+} \cdot A \)
!!
!!### Description
!!
!! This function interface provides methods that return the Moore-Penrose pseudo-inverse of a matrix.
!! Supported data types include `real` and `complex`.
!! The pseudo-inverse \( A^{+} \) is returned as a function result. The computation is based on the
!! singular value decomposition (SVD). An optional relative tolerance `rtol` is provided to control the
!! inclusion of singular values during inversion. Singular values below \( \text{rtol} \cdot \sigma_{\max} \)
!! are treated as zero, where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided,
!! a default threshold is applied.
!!
!! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop`
!! if the state flag `err` is not provided.
!!
!!@note The provided functions are intended for both rectangular and square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] Relative tolerance for singular value cutoff
real(${rk}$), optional, intent(in) :: rtol
!> [optional] State return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
!> Output matrix pseudo-inverse [n,m]
${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp))
end function stdlib_linalg_pseudoinverse_${ri}$
#:endfor
end interface pinv
! Moore-Penrose Pseudo-Inverse: Subroutine interface
interface pseudoinvert
!! version: experimental
!!
!! Computation of the Moore-Penrose pseudo-inverse
!! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-moore-penrose-pseudo-inverse-of-a-matrix))
!!
!!### Summary
!! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a rectangular
!! or square `real` or `complex` matrix.
!! The pseudo-inverse \( A^{+} \) generalizes the matrix inverse and satisfies the properties:
!! - \( A \cdot A^{+} \cdot A = A \)
!! - \( A^{+} \cdot A \cdot A^{+} = A^{+} \)
!! - \( (A \cdot A^{+})^T = A \cdot A^{+} \)
!! - \( (A^{+} \cdot A)^T = A^{+} \cdot A \)
!!
!!### Description
!!
!! This subroutine interface provides a way to compute the Moore-Penrose pseudo-inverse of a matrix.
!! Supported data types include `real` and `complex`.
!! Users must provide two matrices: the input matrix `a` [m,n] and the output pseudo-inverse `pinva` [n,m].
!! The input matrix `a` is used to compute the pseudo-inverse and is not modified. The computed
!! pseudo-inverse is stored in `pinva`. The computation is based on the singular value decomposition (SVD).
!!
!! An optional relative tolerance `rtol` is used to control the inclusion of singular values in the
!! computation. Singular values below \( \text{rtol} \cdot \sigma_{\max} \) are treated as zero,
!! where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided, a default
!! threshold is applied.
!!
!! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop`
!! if the state flag `err` is not provided.
!!
!!@note The provided subroutines are intended for both rectangular and square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err)
!> Input matrix a[m,n]
${rt}$, intent(inout) :: a(:,:)
!> Output pseudo-inverse matrix [n,m]
${rt}$, intent(out) :: pinva(:,:)
!> [optional] Relative tolerance for singular value cutoff
real(${rk}$), optional, intent(in) :: rtol
!> [optional] State return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_pseudoinvert_${ri}$
#:endfor
end interface pseudoinvert
! Moore-Penrose Pseudo-Inverse: Operator interface
interface operator(.pinv.)
!! version: experimental
!!
!! Pseudo-inverse operator of a matrix
!! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-operator))
!!
!!### Summary
!! Operator interface for computing the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix.
!!
!!### Description
!!
!! This operator interface provides a convenient way to compute the Moore-Penrose pseudo-inverse
!! of a matrix. Supported data types include `real` and `complex`. The pseudo-inverse \( A^{+} \)
!! is computed using singular value decomposition (SVD), with singular values below an internal
!! threshold treated as zero.
!!
!! For computational errors or invalid input, the function may return a matrix filled with NaNs.
!!
!!@note The provided functions are intended for both rectangular and square matrices.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Result pseudo-inverse matrix
${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp))
end function stdlib_linalg_pinv_${ri}$_operator
#:endfor
end interface operator(.pinv.)
! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors
interface eig
!! version: experimental
!!
!! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for square matrix \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#eig-eigenvalues-and-eigenvectors-of-a-square-matrix))
!!
!!### Summary
!! Subroutine interface for computing eigenvalues and eigenvectors of a square matrix.
!!
!!### Description
!!
!! This interface provides methods for computing the eigenvalues, and optionally eigenvectors,
!! of a general square matrix. Supported data types include `real` and `complex`, and no assumption is
!! made on the matrix structure. The user may request either left, right, or both
!! eigenvectors to be returned. They are returned as columns of a square matrix with the same size as `A`.
!! Preallocated space for both eigenvalues `lambda` and the eigenvector matrices must be user-provided.
!!
!!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err)
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
!! and optionally right or left eigenvectors.
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Array of eigenvalues
complex(${rk}$), intent(out) :: lambda(:)
!> The columns of RIGHT contain the right eigenvectors of A
complex(${rk}$), optional, intent(out), target :: right(:,:)
!> The columns of LEFT contain the left eigenvectors of A
complex(${rk}$), optional, intent(out), target :: left(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_eig_${ri}$
#:endif
#:endfor
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if rk!="xdp"
module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err)
!! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues,
!! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
!! non-trivial imaginary parts.
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Array of real eigenvalues
real(${rk}$), intent(out) :: lambda(:)
!> The columns of RIGHT contain the right eigenvectors of A
complex(${rk}$), optional, intent(out), target :: right(:,:)
!> The columns of LEFT contain the left eigenvectors of A
complex(${rk}$), optional, intent(out), target :: left(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_real_eig_${ri}$
#:endif
#:endfor
end interface eig
! Eigenvalues of a square matrix
interface eigvals
!! version: experimental
!!
!! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for square matrix \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#eigvals-eigenvalues-of-a-square-matrix))
!!
!!### Summary
!! Function interface for computing the eigenvalues of a square matrix.
!!
!!### Description
!!
!! This interface provides functions for returning the eigenvalues of a general square matrix.
!! Supported data types include `real` and `complex`, and no assumption is made on the matrix structure.
!! An `error stop` is thrown in case of failure; otherwise, error information can be returned
!! as an optional `type(linalg_state_type)` output flag.
!!
!!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Array of singular values
complex(${rk}$), allocatable :: lambda(:)
end function stdlib_linalg_eigvals_${ri}$
module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Array of singular values
complex(${rk}$), allocatable :: lambda(:)
end function stdlib_linalg_eigvals_noerr_${ri}$
#:endif
#:endfor
end interface eigvals
! Eigendecomposition of a real symmetric or complex hermitian matrix
interface eigh
!! version: experimental
!!
!! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for a real symmetric
!! \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix.
!! ([Specification](../page/specs/stdlib_linalg.html#eigh-eigenvalues-and-eigenvectors-of-a-real-symmetric-or-complex-hermitian-square-matrix))
!!
!!### Summary
!! Subroutine interface for computing eigenvalues and eigenvectors of a real symmetric or complex Hermitian square matrix.
!!
!!### Description
!!
!! This interface provides methods for computing the eigenvalues, and optionally eigenvectors,
!! of a real symmetric or complex Hermitian square matrix. Supported data types include `real` and `complex`.
!! The matrix must be symmetric (if `real`) or Hermitian (if `complex`). Only the lower or upper
!! half of the matrix is accessed, and the user can select which using the optional `upper_a`
!! flag (default: use lower half). The vectors are orthogonal, and may be returned as columns of an optional
!! matrix `vectors` with the same kind and size as `A`.
!! Preallocated space for both eigenvalues `lambda` and the eigenvector matrix must be user-provided.
!!
!!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err)
!! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda`
!! of eigenvalues, and optionally right or left eigenvectors.
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Array of eigenvalues
real(${rk}$), intent(out) :: lambda(:)
!> The columns of vectors contain the orthonormal eigenvectors of A
${rt}$, optional, intent(out), target :: vectors(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Should the upper/lower half of A be used? Default: lower
logical(lk), optional, intent(in) :: upper_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_eigh_${ri}$
#:endif
#:endfor
end interface eigh
! Eigenvalues of a real symmetric or complex hermitian matrix
interface eigvalsh
!! version: experimental
!!
!! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for a real
!! symmetric \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix.
!! ([Specification](../page/specs/stdlib_linalg.html#eigvalsh-eigenvalues-of-a-real-symmetric-or-complex-hermitian-square-matrix))
!!
!!### Summary
!! Function interface for computing the eigenvalues of a real symmetric or complex hermitian square matrix.
!!
!!### Description
!!
!! This interface provides functions for returning the eigenvalues of a real symmetric or complex Hermitian
!! square matrix. Supported data types include `real` and `complex`. The matrix must be symmetric
!! (if `real`) or Hermitian (if `complex`). Only the lower or upper half of the matrix is accessed,
!! and the user can select which using the optional `upper_a` flag (default: use lower half).
!! An `error stop` is thrown in case of failure; otherwise, error information can be returned
!! as an optional `type(linalg_state_type)` output flag.
!!
!!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda)
!! Return an array of eigenvalues of real symmetric / complex hermitian A
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] Should the upper/lower half of A be used? Default: lower
logical(lk), optional, intent(in) :: upper_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Array of singular values
real(${rk}$), allocatable :: lambda(:)
end function stdlib_linalg_eigvalsh_${ri}$
module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda)
!! Return an array of eigenvalues of real symmetric / complex hermitian A
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] Should the upper/lower half of A be used? Default: lower
logical(lk), optional, intent(in) :: upper_a
!> Array of singular values
real(${rk}$), allocatable :: lambda(:)
end function stdlib_linalg_eigvalsh_noerr_${ri}$
#:endif
#:endfor
end interface eigvalsh
! Singular value decomposition
interface svd
!! version: experimental
!!
!! Computes the singular value decomposition of a `real` or `complex` 2d matrix.
!! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-rank-2-array-matrix))
!!
!!### Summary
!! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix.
!!
!!### Description
!!
!! This interface provides methods for computing the singular value decomposition of a matrix.
!! Supported data types include `real` and `complex`. The subroutine returns a `real` array of
!! singular values, and optionally, left- and right- singular vector matrices, `U` and `V`.
!! For a matrix `A` with size [m,n], full matrix storage for `U` and `V` should be [m,m] and [n,n].
!! It is possible to use partial storage [m,k] and [k,n], `k=min(m,n)`, choosing `full_matrices=.false.`.
!!
!!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
!!
!!### Example
!!
!!```fortran
!! real(sp) :: a(2,3), s(2), u(2,2), vt(3,3)
!! a = reshape([3,2, 2,3, 2,-2],[2,3])
!!
!! call svd(A,s,u,v)
!! print *, 'singular values = ',s
!!```
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err)
!!### Summary
!! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \)
!!
!!### Description
!!
!! This function computes the singular value decomposition of a `real` or `complex` matrix \( A \),
!! and returns the array of singular values, and optionally the left matrix \( U \) containing the
!! left unitary singular vectors, and the right matrix \( V^T \), containing the right unitary
!! singular vectors.
!!
!! param: a Input matrix of size [m,n].
!! param: s Output `real` array of size [min(m,n)] returning a list of singular values.
!! param: u [optional] Output left singular matrix of size [m,m] or [m,min(m,n)] (.not.full_matrices). Contains singular vectors as columns.
!! param: vt [optional] Output right singular matrix of size [n,n] or [min(m,n),n] (.not.full_matrices). Contains singular vectors as rows.
!! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten.
!! param: full_matrices [optional] If `.true.` (default), matrices \( U \) and \( V^T \) have size [m,m], [n,n]. Otherwise, they are [m,k], [k,n] with `k=min(m,n)`.
!! param: err [optional] State return flag.
!!
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Array of singular values
real(${rk}$), intent(out) :: s(:)
!> The columns of U contain the left singular vectors
${rt}$, optional, intent(out), target :: u(:,:)
!> The rows of V^T contain the right singular vectors
${rt}$, optional, intent(out), target :: vt(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] full matrices have shape(u)==[m,m], shape(vh)==[n,n] (default); otherwise
!> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n)
logical(lk), optional, intent(in) :: full_matrices
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_svd_${ri}$
#:endfor
end interface svd
! Singular values
interface svdvals
!! version: experimental
!!
!! Computes the singular values of a `real` or `complex` 2d matrix.
!! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-rank-2-array-matrix))
!!
!!### Summary
!!
!! Function interface for computing the array of singular values from the singular value decomposition
!! of a `real` or `complex` 2d matrix.
!!
!!### Description
!!
!! This interface provides methods for computing the singular values a 2d matrix.
!! Supported data types include `real` and `complex`. The function returns a `real` array of
!! singular values, with size [min(m,n)].
!!
!!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
!!
!!### Example
!!
!!```fortran
!! real(sp) :: a(2,3), s(2)
!! a = reshape([3,2, 2,3, 2,-2],[2,3])
!!
!! s = svdvals(A)
!! print *, 'singular values = ',s
!!```
!!
#:for rk,rt,ri in RC_KINDS_TYPES
module function stdlib_linalg_svdvals_${ri}$(a,err) result(s)
!!### Summary
!! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \).
!!
!!### Description
!!
!! This function returns the array of singular values from the singular value decomposition of a `real`
!! or `complex` matrix \( A = U \cdot S \cdot V^T \).
!!
!! param: a Input matrix of size [m,n].
!! param: err [optional] State return flag.
!!
!!### Return value
!!
!! param: s `real` array of size [min(m,n)] returning a list of singular values.
!!
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
!> Array of singular values
real(${rk}$), allocatable :: s(:)
end function stdlib_linalg_svdvals_${ri}$
#:endfor
end interface svdvals
#! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2')
#:set NORM_INPUT_TYPE = ["character(len=*)","integer(ilp)"]
#:set NORM_INPUT_SHORT = ["char","int"]
#:set NORM_INPUT_OPTIONS = list(zip(NORM_INPUT_TYPE,NORM_INPUT_SHORT))
! Vector norms: function interface
interface norm
!! version: experimental
!!
!! Computes the vector norm of a generic-rank array \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#norm-computes-the-vector-norm-of-a-generic-rank-array))
!!
!!### Summary
!! Return one of several scalar norm metrics of a `real` or `complex` input array \( A \),
!! that can have any rank. For generic rank-n arrays, the scalar norm over the whole
!! array is returned by default. If `n>=2` and the optional input dimension `dim` is specified,
!! a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D array norms
!! evaluated along dimension `dim` only.
!!
!!
!!### Description
!!
!! This interface provides methods for computing the vector norm(s) of an array.
!! Supported data types include `real` and `complex`.
!! Input arrays may have generic rank from 1 to ${MAXRANK}$.
!!
!! Norm type input is mandatory, and it is provided via the `order` argument.
!! This can be provided as either an `integer` value or a `character` string.
!! Allowed metrics are:
!! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1'
!! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2'
!! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3
!! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf'
!! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf'
!!
!!### Example
!!
!!```fortran
!!
!! real(sp) :: a(3,3), na, rown(3)
!! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!!
!! ! L2 norm: whole matrix
!! na = norm(a, 2)
!!
!! ! Infinity norm of each row
!! rown = norm(a, 'inf', dim=2)
!!
!!```
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:for it,ii in NORM_INPUT_OPTIONS
!> Scalar norms: ${rt}$
#:for rank in range(1, MAXRANK + 1)
pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm)
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
${rt}$, intent(in) :: a${ranksuffix(rank)}$
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> Norm of the matrix.
real(${rk}$) :: nrm
end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$
module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm)
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
${rt}$, intent(in) :: a${ranksuffix(rank)}$
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> Output state return flag.
type(linalg_state_type), intent(out) :: err
!> Norm of the matrix.
real(${rk}$) :: nrm
end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$
#:endfor
!> Array norms: ${rt}$
#:for rank in range(2, MAXRANK + 1)
pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm)
!> Input matrix a[..]
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> Dimension the norm is computed along
integer(ilp), intent(in) :: dim
!> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$
end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm)
!> Input matrix a[..]
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> Dimension the norm is computed along
integer(ilp), intent(in) :: dim
!> Output state return flag.
type(linalg_state_type), intent(out) :: err
!> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$
end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$
#:endfor
#:endfor
#:endfor
end interface norm
!> Vector norm: subroutine interface
interface get_norm
!! version: experimental
!!
!! Computes the vector norm of a generic-rank array \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#get-norm-computes-the-vector-norm-of-a-generic-rank-array))
!!
!!### Summary
!! Subroutine interface that returns one of several scalar norm metrics of a `real` or `complex`
!! input array \( A \), that can have any rank. For generic rank-n arrays, the scalar norm over
!! the whole array is returned by default. If `n>=2` and the optional input dimension `dim` is
!! specified, a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D
!! array norms evaluated along dimension `dim` only.
!!
!!
!!### Description
!!
!! This `pure subroutine `interface provides methods for computing the vector norm(s) of an array.
!! Supported data types include `real` and `complex`.
!! Input arrays may have generic rank from 1 to ${MAXRANK}$.
!!
!! Norm type input is mandatory, and it is provided via the `order` argument.
!! This can be provided as either an `integer` value or a `character` string.
!! Allowed metrics are:
!! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1'
!! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2'
!! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3
!! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf'
!! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf'
!!
!!### Example
!!
!!```fortran
!!
!! real(sp) :: a(3,3), na, rown(3)
!! type(linalg_state_type) :: err
!! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!!
!! ! L2 norm: whole matrix
!! call get_norm(a, na, 2)
!!
!! ! Infinity norms of each row, with error control
!! call get_norm(a, rown, 'inf', dim=2, err=err)
!!
!!```
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:for it,ii in NORM_INPUT_OPTIONS
!> Scalar norms: ${rt}$
#:for rank in range(1, MAXRANK + 1)
pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err)
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
!> Norm of the matrix.
real(${rk}$), intent(out) :: nrm
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out), optional :: err
end subroutine norm_${rank}$D_${ii}$_${ri}$
#:endfor
!> Array norms: ${rt}$
#:for rank in range(2, MAXRANK + 1)
pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err)
!> Input matrix a[..]
${rt}$, intent(in) :: a${ranksuffix(rank)}$
!> Dimension the norm is computed along
integer(ilp), intent(in) :: dim
!> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$
!> Order of the matrix norm being computed.
${it}$, intent(in) :: order
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out), optional :: err
end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
#:endfor
#:endfor
#:endfor
end interface get_norm
!> Matrix norms: function interface
interface mnorm
!! version: experimental
!!
!! Computes the matrix norm of a generic-rank array \( A \).
!! ([Specification](../page/specs/stdlib_linalg.html#mnorm-computes-the-matrix-norm-of-a-generic-rank-array))
!!
!!### Summary
!! Return one of several matrix norm metrics of a `real` or `complex` input array \( A \),
!! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned.
!! If rank>2 and the optional input dimensions `dim` are specified,
!! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all
!! matrix norms evaluated over the specified dimensions only. `dim==[1,2]` are assumed as default
!! dimensions if not specified.
!!
!!### Description
!!
!! This interface provides methods for computing the matrix norm(s) of an array.
!! Supported data types include `real` and `complex`.
!! Input arrays must have rank >= 2.
!!
!! Norm type input is optional, and it is provided via the `order` argument.
!! This can be provided as either an `integer` value or a `character` string.
!! Allowed metrics are:
!! - 1-norm: `order` = 1 or '1'
!! - 2-norm: `order` = 2 or '2'
!! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified
!! - Infinity norm: `order` = huge(0) or 'Inf'
!!
!! If an invalid norm type is provided, the routine returns an error state.
!!
!!### Example
!!
!!```fortran
!! real(sp) :: a(3,3), na
!! real(sp) :: b(3,3,4), nb(4) ! Array of 4 3x3 matrices
!! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!!
!! ! Euclidean/Frobenius norm of single matrix
!! na = mnorm(a)
!! na = mnorm(a, 'Euclidean')
!!
!! ! 1-norm of each 3x3 matrix in b
!! nb = mnorm(b, 1, dim=[1,2])
!!
!! ! Infinity-norm
!! na = mnorm(b, 'inf', dim=[3,2])
!!```
!!
#:for rk,rt,ri in RC_KINDS_TYPES
#:for it,ii in NORM_INPUT_OPTIONS
!> Matrix norms: ${rt}$ rank-2 arrays
module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm)
!> Input matrix a(m,n)
${rt}$, intent(in), target :: a(:,:)
!> Norm of the matrix.
real(${rk}$) :: nrm
!> Order of the matrix norm being computed.
${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out), optional :: err
end function matrix_norm_${ii}$_${ri}$
!> Matrix norms: ${rt}$ higher rank arrays
#:for rank in range(3, MAXRANK + 1)
module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm)
!> Input matrix a(m,n)
${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$
!> Norm of the matrix.
real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$
!> Order of the matrix norm being computed.
${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
!> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2])
integer(ilp), optional, intent(in) :: dim(2)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out), optional :: err
end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$
#:endfor
#:endfor
#:endfor
end interface mnorm
contains
!> Version: experimental
!>
!> Constructs the identity matrix.
!> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))
#:for k1, t1 in RCI_KINDS_TYPES
pure function eye_${t1[0]}$${k1}$(dim1, dim2, mold) result(result)
integer, intent(in) :: dim1
integer, intent(in), optional :: dim2
${t1}$, intent(in) #{if t1 == 'real(dp)'}#, optional #{endif}#:: mold
${t1}$, allocatable :: result(:, :)
integer :: dim2_
integer :: i
dim2_ = optval(dim2, dim1)
allocate(result(dim1, dim2_))
result = 0
do i = 1, min(dim1, dim2_)
result(i, i) = 1
end do
end function eye_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
function trace_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
${t1}$ :: res
integer :: i
res = 0
do i = 1, minval(shape(A))
res = res + A(i,i)
end do
end function trace_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
pure function is_square_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
res = (size(A,1) == size(A,2))
end function is_square_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
pure function is_diagonal_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
${t1}$, parameter :: zero = 0 !zero of relevant type
integer :: m, n, o, i, j
m = size(A,1)
n = size(A,2)
do j = 1, n !loop over all columns
o = min(j-1,m) !index of row above diagonal (or last row)
do i = 1, o !loop over rows above diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
do i = o+2, m !loop over rows below diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
end do
res = .true. !otherwise A is diagonal
end function is_diagonal_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
pure function is_symmetric_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
integer :: n, i, j
if (.not. is_square(A)) then
res = .false.
return !nonsquare matrices cannot be symmetric
end if
n = size(A,1) !symmetric dimension of A
do j = 1, n !loop over all columns
do i = 1, j-1 !loop over all rows above diagonal
if (A(i,j) /= A(j,i)) then
res = .false.
return
end if
end do
end do
res = .true. !otherwise A is symmetric
end function is_symmetric_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
pure function is_skew_symmetric_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
integer :: n, i, j
if (.not. is_square(A)) then
res = .false.
return !nonsquare matrices cannot be skew-symmetric
end if
n = size(A,1) !symmetric dimension of A
do j = 1, n !loop over all columns
do i = 1, j !loop over all rows above diagonal (and diagonal)
if (A(i,j) /= -A(j,i)) then
res = .false.
return
end if
end do
end do
res = .true. !otherwise A is skew-symmetric
end function is_skew_symmetric_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in (REAL_KINDS_TYPES + INT_KINDS_TYPES)
pure function is_hermitian_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
res = is_symmetric(A) !symmetry and Hermiticity are equivalent for real matrices
end function is_hermitian_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
pure function is_hermitian_${t1[0]}$${k1}$(A) result(res)
${t1}$, intent(in) :: A(:,:)
logical :: res
integer :: n, i, j
if (.not. is_square(A)) then
res = .false.
return !nonsquare matrices cannot be Hermitian
end if
n = size(A,1) !symmetric dimension of A
do j = 1, n !loop over all columns
do i = 1, j !loop over all rows above diagonal (and diagonal)
if (A(i,j) /= conjg(A(j,i))) then
res = .false.
return
end if
end do
end do
res = .true. !otherwise A is Hermitian
end function is_hermitian_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
pure module function hermitian_${t1[0]}$${k1}$(a) result(ah)
${t1}$, intent(in) :: a(:,:)
${t1}$ :: ah(size(a, 2), size(a, 1))
#:if t1.startswith('complex')
ah = conjg(transpose(a))
#:else
ah = transpose(a)
#:endif
end function hermitian_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
function is_triangular_${t1[0]}$${k1}$(A,uplo) result(res)
${t1}$, intent(in) :: A(:,:)
character, intent(in) :: uplo
logical :: res
${t1}$, parameter :: zero = 0 !zero of relevant type
integer :: m, n, o, i, j
m = size(A,1)
n = size(A,2)
if ((uplo == 'u') .or. (uplo == 'U')) then !check for upper triangularity
do j = 1, n !loop over all columns
o = min(j-1,m) !index of row above diagonal (or last row)
do i = o+2, m !loop over rows below diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
end do
else if ((uplo == 'l') .or. (uplo == 'L')) then !check for lower triangularity
do j=1,n !loop over all columns
o = min(j-1,m) !index of row above diagonal (or last row)
do i=1,o !loop over rows above diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
end do
else
call error_stop("ERROR (is_triangular): second argument must be one of {'u','U','l','L'}")
end if
res = .true. !otherwise A is triangular of the requested type
end function is_triangular_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in RCI_KINDS_TYPES
function is_hessenberg_${t1[0]}$${k1}$(A,uplo) result(res)
${t1}$, intent(in) :: A(:,:)
character, intent(in) :: uplo
logical :: res
${t1}$, parameter :: zero = 0 !zero of relevant type
integer :: m, n, o, i, j
m = size(A,1)
n = size(A,2)
if ((uplo == 'u') .or. (uplo == 'U')) then !check for upper Hessenberg
do j = 1, n !loop over all columns
o = min(j-2,m) !index of row two above diagonal (or last row)
do i = o+4, m !loop over rows two or more below main diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
end do
else if ((uplo == 'l') .or. (uplo == 'L')) then !check for lower Hessenberg
do j = 1, n !loop over all columns
o = min(j-2,m) !index of row two above diagonal (or last row)
do i = 1, o !loop over rows one or more above main diagonal
if (A(i,j) /= zero) then
res = .false.
return
end if
end do
end do
else
call error_stop("ERROR (is_hessenberg): second argument must be one of {'u','U','l','L'}")
end if
res = .true. !otherwise A is Hessenberg of the requested type
end function is_hessenberg_${t1[0]}$${k1}$
#:endfor
end module stdlib_linalg
