#:include "common.fypp"
#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
module stdlib_linalg
!!Provides a support for various linear algebra procedures
!! ([Specification](../page/specs/stdlib_linalg.html))
use stdlib_kinds, only: sp, dp, xdp, qp, &
int8, int16, int32, int64
use stdlib_error, only: error_stop
use stdlib_optval, only: optval
implicit none
private
public :: diag
public :: eye
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
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
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
