Merge pull request #798 from perazz/determinant

jvdp1 · web-flow · commit 28ae6e03398b · 2024-04-25T12:48:40.000-04:00
linalg: determinant
diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
@@ -599,3 +599,77 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l
 ```fortran
 {!example/linalg/example_is_hessenberg.f90!}
 ```
+
+## `det` - Computes the determinant of a square matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the determinant of a `real` or `complex` square matrix.
+
+This interface comes with a `pure` version `det(a)`, and a non-pure version `det(a,overwrite_a,err)` that
+allows for more expert control.
+
+### Syntax
+
+`c = ` [[stdlib_linalg(module):det(interface)]] `(a [, overwrite_a, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2 square array
+
+`overwrite_a` (optional): Shall be an input `logical` flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation.
+ This is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value.  This is an `intent(out)` argument.
+
+### Return value
+
+Returns a `real` scalar value of the same kind of `a` that represents the determinant of the matrix.
+
+Raises `LINALG_ERROR` if the matrix is singular.
+Raises `LINALG_VALUE_ERROR` if the matrix is non-square.
+Exceptions are returned to the `err` argument if provided; an `error stop` is triggered otherwise.
+
+### Example
+
+```fortran
+{!example/linalg/example_determinant.f90!}
+```
+
+## `.det.` - Determinant operator of a square matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This operator returns the determinant of a real square matrix.
+
+This interface is equivalent to the `pure` version of determinant [[stdlib_linalg(module):det(interface)]]. 
+
+### Syntax
+
+`c = ` [[stdlib_linalg(module):operator(.det.)(interface)]] `(a)`
+
+### Arguments
+
+`a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument.
+
+### Return value
+
+Returns a real scalar value that represents the determinnt of the matrix.
+
+Raises `LINALG_ERROR` if the matrix is singular.
+Raises `LINALG_VALUE_ERROR` if the matrix is non-square.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_determinant2.f90!}
+```
diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt
@@ -18,3 +18,6 @@ ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
+ADD_EXAMPLE(determinant)
+ADD_EXAMPLE(determinant2)
+
diff --git a/example/linalg/example_determinant.f90 b/example/linalg/example_determinant.f90
@@ -0,0 +1,14 @@
+program example_determinant
+  use stdlib_kinds, only: dp
+  use stdlib_linalg, only: det, linalg_state_type
+  implicit none
+  type(linalg_state_type) :: err
+
+  real(dp) :: d
+
+  ! Compute determinate of a real matrix
+  d = det(reshape([real(dp)::1,2,3,4],[2,2]))
+
+  print *, d ! a*d-b*c = -2.0
+
+end program example_determinant
diff --git a/example/linalg/example_determinant2.f90 b/example/linalg/example_determinant2.f90
@@ -0,0 +1,13 @@
+program example_determinant2
+  use stdlib_kinds, only: dp
+  use stdlib_linalg, only: operator(.det.)
+  implicit none
+
+  real(dp) :: d
+
+  ! Compute determinate of a real matrix
+  d = .det.reshape([real(dp)::1,2,3,4],[2,2])
+
+  print *, d ! a*d-b*c = -2.0
+
+end program example_determinant2
diff --git a/src/CMakeLists.txt b/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
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
@@ -1,15 +1,19 @@
 #:include "common.fypp"
-#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set RCI_KINDS_TYPES = RC_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, &
+  use stdlib_kinds, only: sp, dp, xdp, qp, lk, &
     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
   implicit none
   private
 
+  public :: det
+  public :: operator(.det.)
   public :: diag
   public :: eye
   public :: trace
@@ -23,6 +27,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
@@ -214,6 +221,109 @@ module stdlib_linalg
     #:endfor
   end interface is_hessenberg
 
+
+  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
 
 
diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp
diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp
