linalg: determinant

src/CMakeLists.txt

@@ -24,6 +24,7 @@ set(fppFiles
     stdlib_linalg_outer_product.fypp
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
+    stdlib_linalg_determinant.fypp
     stdlib_linalg_state.fypp 
     stdlib_optval.fypp
     stdlib_selection.fypp

src/stdlib_linalg.fypp

@@ -7,9 +7,12 @@ module stdlib_linalg
     int8, int16, int32, int64
   use stdlib_error, only: error_stop
   use stdlib_optval, only: optval
+  use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
+  use stdlib_linalg_determinant, only: det
   implicit none
   private
 
+  public :: det
   public :: diag
   public :: eye
   public :: trace
@@ -23,6 +26,9 @@ module stdlib_linalg
   public :: is_hermitian
   public :: is_triangular
   public :: is_hessenberg
+
+  ! Export linalg error handling
+  public :: linalg_state_type, linalg_error_handling
 
   interface diag
     !! version: experimental

src/stdlib_linalg_determinant.fypp

@@ -0,0 +1,302 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+module stdlib_linalg_determinant
+!! Determinant of a rectangular matrix
+     use stdlib_linalg_constants
+     use stdlib_linalg_lapack, only: getrf
+     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
+     implicit none(type,external)
+     private
+
+     ! Function interface
+     public :: det
+     public :: operator(.det.)
+
+     character(*), parameter :: this = 'determinant'
+
+     interface det
+     !!### 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.          
+     !! 
+     !!### 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.)
+     !!### 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.
+     !!
+     !!### 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.)
+
+     contains
+
+     #:for rk,rt in RC_KINDS_TYPES
+     #:if rk!="xdp"
+     pure function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det)     
+     !!### Summary
+     !! Compute determinant of a real square matrix (pure interface).
+     !!
+     !!### Description
+     !!
+     !! This function computes the determinant of a real square matrix.
+     !!
+     !! param: a Input matrix of size [m,n].
+     !! return: det Matrix determinant.
+     !!
+     !!### Example
+     !!
+     !!```fortran
+     !!
+     !! ${rt}$ :: matrix(3,3)
+     !! ${rt}$ :: determinant
+     !! matrix = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+     !! determinant = det(matrix)         
+     !!
+     !!```
+         !> Input matrix a[m,n]
+         ${rt}$, intent(in) :: a(:,:)
+         !> Matrix determinant
+         ${rt}$ :: det
+
+         !! Local variables
+         type(linalg_state_type) :: err0
+         integer(ilp) :: m,n,info,perm,k
+         integer(ilp), allocatable :: ipiv(:)
+         ${rt}$, allocatable :: amat(:,:)
+
+         ! Matrix determinant size
+         m = size(a,1,kind=ilp)
+         n = size(a,2,kind=ilp)
+
+         if (m/=n .or. .not.min(m,n)>=0) then
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']')
+            det  = 0.0_${rk}$
+            goto 1
+         end if
+
+         select case (m)
+            case (0)
+
+                ! Empty array has determinant 1 because math
+                det = 1.0_${rk}$
+
+            case (1)
+
+                ! Scalar input
+                det = a(1,1)
+
+            case default
+
+                ! Find determinant from LU decomposition
+
+                ! Initialize a matrix temporary
+                allocate(amat(m,n),source=a)
+
+                ! Pivot indices
+                allocate(ipiv(n))
+
+                ! Compute determinant from LU factorization, then calculate the 
+                ! product of all diagonal entries of the U factor.
+                call getrf(m,n,amat,m,ipiv,info)
+
+                select case (info)
+                   case (0)
+                       ! Success: compute determinant
+
+                       ! Start with real 1.0
+                       det = 1.0_${rk}$
+                       perm = 0
+                       do k=1,n
+                          if (ipiv(k)/=k) perm = perm+1
+                          det = det*amat(k,k)
+                       end do
+                       if (mod(perm,2)/=0) det = -det
+
+                   case (:-1)
+                       err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
+                   case (1:)
+                       err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+                   case default
+                       err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+                end select
+
+                deallocate(amat)
+
+         end select
+
+         ! Process output and return
+         1 call linalg_error_handling(err0)
+
+     end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
+
+     function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det)
+     !!### Summary
+     !! Compute determinant of a square matrix (with error control).
+     !!
+     !!### Description
+     !!
+     !! This function computes the determinant of a square matrix with error control.
+     !!
+     !! param: a Input matrix of size [m,n].
+     !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten.
+     !! param: err State return flag. 
+     !! return: det Matrix determinant.
+     !!
+     !!### Example
+     !!
+     !!```fortran
+     !!
+     !! ${rt}$ :: matrix(3,3)
+     !! ${rt}$ :: determinant
+     !! matrix = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+     !! determinant = det(matrix, err=err)         
+     !!
+     !!```        
+     !     
+         !> 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
+
+         !! Local variables
+         type(linalg_state_type) :: err0
+         integer(ilp) :: m,n,info,perm,k
+         integer(ilp), allocatable :: ipiv(:)
+         logical(lk) :: copy_a
+         ${rt}$, pointer :: amat(:,:)
+
+         ! Matrix determinant size
+         m = size(a,1,kind=ilp)
+         n = size(a,2,kind=ilp)
+
+         if (m/=n .or. .not.min(m,n)>=0) then
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']')
+            det  = 0.0_${rk}$
+            goto 1
+         end if
+
+         ! Can A be overwritten? By default, do not overwrite
+         if (present(overwrite_a)) then
+            copy_a = .not.overwrite_a
+         else
+            copy_a = .true._lk
+         endif
+
+         select case (m)
+            case (0)
+
+                ! Empty array has determinant 1 because math
+                det = 1.0_${rk}$
+
+            case (1)
+
+                ! Scalar input
+                det = a(1,1)
+
+            case default
+
+                ! Find determinant from LU decomposition
+
+                ! Initialize a matrix temporary
+                if (copy_a) then
+                   allocate(amat(m,n),source=a)
+                else
+                   amat => a
+                endif
+
+                ! Pivot indices
+                allocate(ipiv(n))
+
+                ! Compute determinant from LU factorization, then calculate the 
+                ! product of all diagonal entries of the U factor.
+                call getrf(m,n,amat,m,ipiv,info)
+
+                select case (info)
+                   case (0)
+                       ! Success: compute determinant
+
+                       ! Start with real 1.0
+                       det = 1.0_${rk}$
+                       perm = 0
+                       do k=1,n
+                          if (ipiv(k)/=k) perm = perm+1
+                          det = det*amat(k,k)
+                       end do
+                       if (mod(perm,2)/=0) det = -det
+
+                   case (:-1)
+                       err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
+                   case (1:)
+                       err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+                   case default
+                       err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+                end select
+
+                if (.not.copy_a) deallocate(amat)
+
+         end select
+
+         ! Process output and return
+         1 call linalg_error_handling(err0,err)
+
+     end function stdlib_linalg_${rt[0]}$${rk}$determinant
+
+     #:endif
+     #:endfor
+
+end module stdlib_linalg_determinant

test/linalg/CMakeLists.txt

@@ -2,10 +2,12 @@ set(
   fppFiles
   "test_linalg.fypp"
   "test_blas_lapack.fypp"
+  "test_linalg_determinant.fypp"
   "test_linalg_matrix_property_checks.fypp"
 )
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
 ADDTEST(linalg)
+ADDTEST(linalg_determinant)
 ADDTEST(linalg_matrix_property_checks)
 ADDTEST(blas_lapack)