#: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 :: det
public :: operator(.det.)
public :: diag
public :: eye
public :: lstsq
public :: lstsq_space
public :: solve
public :: solve_lu
public :: solve_lstsq
public :: trace
public :: outer_product
public :: kronecker_product
public :: cross_product
public :: is_square
public :: is_diagonal
public :: is_symmetric
public :: is_skew_symmetric
public :: is_hermitian
public :: is_triangular
public :: is_hessenberg
! Export linalg error handling
public :: linalg_state_type, linalg_error_handling
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
! 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
! 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`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
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}$
#:endif
#: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`.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
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}$
#:endif
#: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.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
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}$
#:endif
#: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.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
#:for nd,ndsuf,nde in ALL_RHS
#:for rk,rt,ri in RC_KINDS_TYPES
#:if rk!="xdp"
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}$
#:endif
#: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
#:if rk!="xdp"
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}$
#:endif
#:endfor
#:endfor
end interface lstsq_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
contains
!> Version: experimental
!>
!> Constructs the identity matrix.
!> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))
pure function eye(dim1, dim2) result(result)
integer, intent(in) :: dim1
integer, intent(in), optional :: dim2
integer(int8), allocatable :: result(:, :)
integer :: dim2_
integer :: i
dim2_ = optval(dim2, dim1)
allocate(result(dim1, dim2_))
result = 0_int8
do i = 1, min(dim1, dim2_)
result(i, i) = 1_int8
end do
end function eye
#: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
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
