diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index ab52dfb71..8a42ae25a 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -595,3 +595,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
index 1a5875502..f515af446 100644
--- 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
new file mode 100644
index 000000000..00ef7026b
--- /dev/null
+++ 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
new file mode 100644
index 000000000..0dcd89c8e
--- /dev/null
+++ 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
index 165259db3..38d0ab569 100644
--- 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
index 3ed905c56..2f837f15a 100644
--- 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
new file mode 100644
index 000000000..0091d596b
--- /dev/null
+++ b/src/stdlib_linalg_determinant.fypp
@@ -0,0 +1,239 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_linalg) 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)
+
+     ! Function interface
+     character(*), parameter :: this = 'determinant'
+
+     contains
+
+     ! BLAS/LAPACK backends do not currently support xdp
+     #:for rk,rt in RC_KINDS_TYPES
+     #:if rk!="xdp" 
+     pure module 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}$
+            ! Process output and return
+            call linalg_error_handling(err0)
+         return
+         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
+         call linalg_error_handling(err0)
+
+     end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
+     
+     module 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}$
+            ! Process output and return
+            call linalg_error_handling(err0,err)
+            return
+         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, 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
+         call linalg_error_handling(err0,err)
+
+     end function stdlib_linalg_${rt[0]}$${rk}$determinant
+
+     #:endif
+     #:endfor
+
+end submodule stdlib_linalg_determinant
diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index 3d590a9d2..92006323e 100644
--- a/test/linalg/CMakeLists.txt
+++ b/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)
diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp
new file mode 100644
index 000000000..6b9310f72
--- /dev/null
+++ b/test/linalg/test_linalg_determinant.fypp
@@ -0,0 +1,163 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test matrix determinant
+module test_linalg_determinant
+    use testdrive, only: error_type, check, new_unittest, unittest_type
+    use stdlib_linalg_constants
+    use stdlib_linalg, only: eye, det, linalg_state_type
+
+    implicit none (type,external)
+    private
+    
+    public :: test_matrix_determinant 
+
+    contains
+    
+
+    !> Matrix inversion tests
+    subroutine test_matrix_determinant(tests)
+        !> Collection of tests
+        type(unittest_type), allocatable, intent(out) :: tests(:)
+        
+        allocate(tests(0))
+
+        #:for rk,rt in RC_KINDS_TYPES
+        #:if rk!="xdp"
+        tests = [tests,new_unittest("$eye_det_${rt[0]}$${rk}$",test_${rt[0]}$${rk}$_eye_determinant)]
+        tests = [tests,new_unittest("$eye_det_multiple_${rt[0]}$${rk}$",test_${rt[0]}$${rk}$_eye_multiple)]
+        #:endif
+        #:endfor
+        #:for ck,ct in CMPLX_KINDS_TYPES
+        #:if ck!="xdp"
+        tests = [tests,new_unittest("$complex_det_${rt[0]}$${rk}$",test_${ct[0]}$${ck}$_complex_determinant)]
+        #:endif
+        #: endfor
+
+    end subroutine test_matrix_determinant
+
+    #:for rk,rt in RC_KINDS_TYPES
+    #:if rk!="xdp"
+    !> Determinant of identity matrix
+    subroutine test_${rt[0]}$${rk}$_eye_determinant(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        type(linalg_state_type) :: state
+
+        integer(ilp) :: i
+        integer(ilp), parameter :: n = 128_ilp
+
+        ${rt}$ :: a(n,n),deta
+
+        a = eye(n)
+
+        !> Determinant function
+        deta = det(a,err=state)
+        
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+        
+        call check(error, abs(deta-1.0_${rk}$)<epsilon(0.0_${rk}$), 'det(eye(n))==1')
+        
+    end subroutine test_${rt[0]}$${rk}$_eye_determinant
+
+    !> Determinant of identity matrix multiplier
+    subroutine test_${rt[0]}$${rk}$_eye_multiple(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        type(linalg_state_type) :: state
+
+        integer(ilp), parameter :: n = 4_ilp
+        real(${rk}$), parameter :: coef = 0.01_${rk}$
+        integer(ilp) :: i,j
+        ${rt}$ :: a(n,n),deta
+
+        !> Multiply eye by a very small number
+        a = eye(n)
+        do concurrent (i=1:n)
+          a(i,i) = coef
+        end do
+
+        !> Determinant: small, but a is not singular, because it is a multiple of the identity.
+        deta = det(a,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+        
+        call check(error, abs(deta-coef**n)<max(tiny(0.0_${rk}$),epsilon(0.0_${rk}$)*coef**n), &
+                         'det(0.01*eye(n))==0.01^n')
+
+    end subroutine test_${rt[0]}$${rk}$_eye_multiple
+    #:endif
+    #:endfor
+
+    !> Determinant of complex identity matrix
+    #:for ck,ct in CMPLX_KINDS_TYPES
+    #:if ck!="xdp"
+    subroutine test_${ct[0]}$${ck}$_complex_determinant(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        type(linalg_state_type) :: state
+
+        integer(ilp) :: i,j,n
+        integer(ilp), parameter :: nmax = 10_ilp
+
+        ${ct}$, parameter :: res(nmax) = [${ct}$::(1,1),(0,2),(-2,2),(-4,0),(-4,-4), &
+                                                  (0,-8),(8,-8),(16,0),(16,16),(0,32)]
+
+        ${ct}$, allocatable :: a(:,:)
+        ${ct}$ :: deta(nmax)
+
+        !> Test determinant for all sizes, 1:nmax
+        matrix_size: do n=1,nmax
+
+           ! Put 1+i on each diagonal element
+           a = eye(n)
+           do concurrent (i=1:n)
+             a(i,i) = (1.0_${ck}$,1.0_${ck}$)
+           end do
+
+           ! Expected result
+           deta(n) = det(a,err=state)
+
+           deallocate(a)
+           if (state%error()) exit matrix_size
+
+        end do matrix_size
+
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+        
+        call check(error, all(abs(res-deta)<=epsilon(0.0_${ck}$)), &
+                         'det((1+i)*eye(n))  does not match result')
+
+    end subroutine test_${ct[0]}$${ck}$_complex_determinant
+
+    #:endif
+    #:endfor
+
+end module test_linalg_determinant
+
+program test_det
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_determinant, only : test_matrix_determinant
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_determinant", test_matrix_determinant) &
+        ]
+
+    do is = 1, size(testsuites)
+        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        call run_testsuite(testsuites(is)%collect, error_unit, stat)
+    end do
+
+    if (stat > 0) then
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+    end if
+end program