[stdlib_linalg] matrix property checks

src/stdlib_linalg.fypp

@@ -5,13 +5,21 @@ module stdlib_linalg
   !! ([Specification](../page/specs/stdlib_linalg.html))
   use stdlib_kinds, only: sp, dp, qp, &
     int8, int16, int32, int64
+  use stdlib_error, only: error_stop
   implicit none
   private
 
   public :: diag
   public :: eye
   public :: trace
   public :: outer_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
@@ -80,8 +88,93 @@ module stdlib_linalg
     #:endfor
   end interface outer_product
 
+
+  ! Check for squareness
+  interface is_square
+    !! version: experimental
+    !!
+    !! Checks if a matrix (rank-2 array) is square
+    !! ([Specification](../page/specs/stdlib_linalg.html#description_4))
+    #: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#description_5))
+    #: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#description_6))
+    #: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#description_7))
+    #: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#description_8))
+    #:for k1, t1 in CMPLX_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#description_9))
+    #: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#description_10))
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module procedure is_Hessenberg_${t1[0]}$${k1}$
+    #:endfor
+  end interface is_hessenberg
+
 contains
 
+
     function eye(n) result(res)
       !! version: experimental
       !!
@@ -108,4 +201,197 @@ contains
         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
+        integer :: A_shape(2)
+        A_shape = shape(A)
+        res = (A_shape(1) .eq. A_shape(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}$ :: zero
+        integer :: A_shape(2), m, n, o, i, j
+        zero = 0 !zero of relevant type
+        A_shape = shape(A)
+        m = A_shape(1)
+        n = A_shape(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 (.not. (A(i,j) .eq. zero)) then
+                  res = .false.
+                  return
+                end if
+            end do
+            do i = o+2, m !loop over rows below diagonal
+                if (.not. (A(i,j) .eq. 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 :: A_shape(2), n, i, j
+        if (.not. is_square(A)) then
+           res = .false.
+           return !nonsquare matrices cannot be symmetric
+        end if
+        A_shape = shape(A)
+        n = A_shape(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 (.not. (A(i,j) .eq. 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 :: A_shape(2), n, i, j
+        if (.not. is_square(A)) then
+           res = .false.
+           return !nonsquare matrices cannot be skew-symmetric
+        end if
+        A_shape = shape(A)
+        n = A_shape(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 (.not. (A(i,j) .eq. -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 CMPLX_KINDS_TYPES
+      pure function is_hermitian_${t1[0]}$${k1}$(A) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        logical :: res
+        integer :: A_shape(2), n, i, j
+        if (.not. is_square(A)) then
+           res = .false.
+           return !nonsquare matrices cannot be Hermitian
+        end if
+        A_shape = shape(A)
+        n = A_shape(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 (.not. (A(i,j) .eq. 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}$ :: zero
+        integer :: A_shape(2), m, n, o, i, j
+        zero = 0 !zero of relevant type
+        A_shape = shape(A)
+        m = A_shape(1)
+        n = A_shape(2)
+        if ((uplo .eq. 'u') .or. (uplo .eq. '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 (.not. (A(i,j) .eq. zero)) then
+                    res = .false.
+                    return
+                  end if
+              end do
+          end do
+        else if ((uplo .eq. 'l') .or. (uplo .eq. '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 (.not. (A(i,j) .eq. 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}$ :: zero
+        integer :: A_shape(2), m, n, o, i, j
+        zero = 0 !zero of relevant type
+        A_shape = shape(A)
+        m = A_shape(1)
+        n = A_shape(2)
+        if ((uplo .eq. 'u') .or. (uplo .eq. '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 (.not. (A(i,j) .eq. zero)) then
+                    res = .false.
+                    return
+                  end if
+              end do
+          end do
+        else if ((uplo .eq. 'l') .or. (uplo .eq. '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 (.not. (A(i,j) .eq. 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