From 4a4b89971f057dc3b93cc74b6626f5fb18cafbdc Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 09:52:57 +0200 Subject: [PATCH 01/28] import files --- src/CMakeLists.txt | 1 + src/stdlib_linalg_determinant.fypp | 124 ++++++++++++++++++++++ test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_determinant.fypp | 127 +++++++++++++++++++++++ 4 files changed, 254 insertions(+) create mode 100644 src/stdlib_linalg_determinant.fypp create mode 100644 test/linalg/test_linalg_determinant.fypp 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_determinant.fypp b/src/stdlib_linalg_determinant.fypp new file mode 100644 index 000000000..6ac6e9b4f --- /dev/null +++ b/src/stdlib_linalg_determinant.fypp @@ -0,0 +1,124 @@ +#:include "common.fypp" +module stdlib_linalg_determinant + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state + use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + implicit none(type,external) + private + + !> Determinant of a rectangular matrix + public :: det + + character(*), parameter :: this = 'determinant' + + ! Numpy: det(a) + ! Scipy: det(a, overwrite_a=False, check_finite=True) + ! IMSL: DET(a) + + interface det + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_${ri}$determinant + #:endfor + end interface det + + + contains + + #:for rk,rt,ri in ALL_KINDS_TYPES + ! Compute determinant of a square matrix A + function stdlib_linalg_${ri}$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 + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + !> Result: matrix determinant + ${rt}$ :: det + + !> Local variables + type(linalg_state) :: 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(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 + 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(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']') + case (1:) + err0 = linalg_state(this,LINALG_ERROR,'singular matrix') + case default + err0 = linalg_state(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_${ri}$determinant + + #:endfor + +end module stdlib_linalg_determinant diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 3d590a9d2..2437c83e7 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_determinant.fypp" "test_linalg_matrix_property_checks.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -9,3 +10,4 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_matrix_property_checks) ADDTEST(blas_lapack) +ADDTEST(linalg_determinant) diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp new file mode 100644 index 000000000..dbf5d4deb --- /dev/null +++ b/test/linalg/test_linalg_determinant.fypp @@ -0,0 +1,127 @@ +#:include "common.fypp" +! Test matrix determinant +module test_linalg_determinant + use stdlib_linalg_interface + + implicit none (type,external) + + contains + + !> Matrix inversion tests + subroutine test_matrix_determinant(error) + logical, intent(out) :: error + + real :: t0,t1 + + call cpu_time(t0) + + #:for rk,rt,ri in ALL_KINDS_TYPES + call test_${ri}$_eye_determinant(error) + if (error) return + + call test_${ri}$_eye_multiple(error) + if (error) return + #: endfor + + #:for ck,ct,ci in CMPL_KINDS_TYPES + call test_${ci}$_complex_determinant(error) + if (error) return + #: endfor + + call cpu_time(t1) + + print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) + + 1 format('Determinant tests completed in ',f9.4,' milliseconds, result=',a) + + end subroutine test_matrix_determinant + + !> Determinant of identity matrix + #:for rk,rt,ri in ALL_KINDS_TYPES + subroutine test_${ri}$_eye_determinant(error) + logical, intent(out) :: error + + type(linalg_state) :: 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) + error = state%error() .or. .not.abs(deta-1.0_${rk}$) Determinant of identity matrix multiplier + subroutine test_${ri}$_eye_multiple(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp), parameter :: n = 10_ilp + real(${rk}$), parameter :: coef = 0.0001_${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) + error = state%error() .or. .not.abs(deta-coef**n) Determinant of complex identity matrix + #:for ck,ct,ci in CMPL_KINDS_TYPES + subroutine test_${ci}$_complex_determinant(error) + logical, intent(out) :: error + + type(linalg_state) :: 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 + + error = state%error() .or. any(.not.abs(res-deta)<=tiny(0.0_${ck}$)) + + end subroutine test_${ci}$_complex_determinant + + #:endfor + +end module test_linalg_determinant + + From d54cfa3c8217b67ce0bf1b55f83d6dd3c16312f4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 10:23:45 +0200 Subject: [PATCH 02/28] add `pure` interface --- src/stdlib_linalg.fypp | 6 ++ src/stdlib_linalg_determinant.fypp | 123 ++++++++++++++++++++++++----- 2 files changed, 108 insertions(+), 21 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 3ed905c56..4047612c8 100644 --- a/src/stdlib_linalg.fypp +++ b/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 diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index 6ac6e9b4f..1c7897f3b 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -1,14 +1,15 @@ #:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +!> Determinant of a rectangular matrix module stdlib_linalg_determinant use stdlib_linalg_constants - use stdlib_linalg_blas - use stdlib_linalg_lapack - use stdlib_linalg_state - use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + 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 - !> Determinant of a rectangular matrix + !> Function interface public :: det character(*), parameter :: this = 'determinant' @@ -18,28 +19,106 @@ module stdlib_linalg_determinant ! IMSL: DET(a) interface det - #:for rk,rt,ri in ALL_KINDS_TYPES - module procedure stdlib_linalg_${ri}$determinant + #:for rk,rt in RC_KINDS_TYPES + module procedure stdlib_linalg_${rt[0]}$${rk}$determinant + module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endfor end interface det - contains - #:for rk,rt,ri in ALL_KINDS_TYPES - ! Compute determinant of a square matrix A - function stdlib_linalg_${ri}$determinant(a,overwrite_a,err) result(det) + #:for rk,rt in RC_KINDS_TYPES + ! Compute determinant of a square matrix A: pure interface + function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Result: 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 + + ! Compute determinant of a square matrix A, with error control + 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 !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), optional, intent(out) :: err + type(linalg_state_type), intent(out) :: err !> Result: matrix determinant ${rt}$ :: det !> Local variables - type(linalg_state) :: err0 + type(linalg_state_type) :: err0 integer(ilp) :: m,n,info,perm,k integer(ilp), allocatable :: ipiv(:) logical(lk) :: copy_a @@ -50,7 +129,7 @@ module stdlib_linalg_determinant n = size(a,2,kind=ilp) if (m/=n .or. .not.min(m,n)>=0) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']') det = 0.0_${rk}$ goto 1 end if @@ -64,11 +143,13 @@ module stdlib_linalg_determinant select case (m) case (0) + ! Empty array has determinant 1 because math det = 1.0_${rk}$ case (1) - ! Scalar + + ! Scalar input det = a(1,1) case default @@ -85,8 +166,8 @@ module stdlib_linalg_determinant ! Pivot indices allocate(ipiv(n)) - ! Compute determinant from LU factorization, then calculate the product of - ! all diagonal entries of the U factor. + ! 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) @@ -103,11 +184,11 @@ module stdlib_linalg_determinant if (mod(perm,2)/=0) det = -det case (:-1) - err0 = linalg_state(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']') + err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']') case (1:) - err0 = linalg_state(this,LINALG_ERROR,'singular matrix') + err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix') case default - err0 = linalg_state(this,LINALG_INTERNAL_ERROR,'catastrophic error') + err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') end select if (.not.copy_a) deallocate(amat) @@ -117,7 +198,7 @@ module stdlib_linalg_determinant ! Process output and return 1 call linalg_error_handling(err0,err) - end function stdlib_linalg_${ri}$determinant + end function stdlib_linalg_${rt[0]}$${rk}$determinant #:endfor From c987fa10a08cd8b7c9a5377040ec05dc9583c9dc Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 10:24:41 +0200 Subject: [PATCH 03/28] add operator interface --- src/stdlib_linalg_determinant.fypp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index 1c7897f3b..bce5c248b 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -11,6 +11,7 @@ module stdlib_linalg_determinant !> Function interface public :: det + public :: operator(.det.) character(*), parameter :: this = 'determinant' @@ -20,10 +21,19 @@ module stdlib_linalg_determinant interface det #:for rk,rt in RC_KINDS_TYPES + ! Interface with error control module procedure stdlib_linalg_${rt[0]}$${rk}$determinant + ! Pure interface module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endfor end interface det + + ! Pure Operator interface + interface operator(.det.) + #:for rk,rt in RC_KINDS_TYPES + module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant + #:endfor + end interface operator(.det.) contains From b9a3b0e939c49b1a38166dcee744c0c89180fd4d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 11:11:36 +0200 Subject: [PATCH 04/28] add tests add tests add test file --- test/linalg/CMakeLists.txt | 2 +- test/linalg/test_linalg_determinant.fypp | 121 ++++++++++++++--------- 2 files changed, 76 insertions(+), 47 deletions(-) diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 2437c83e7..92006323e 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -8,6 +8,6 @@ set( fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) +ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) ADDTEST(blas_lapack) -ADDTEST(linalg_determinant) diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp index dbf5d4deb..7c3749e2c 100644 --- a/test/linalg/test_linalg_determinant.fypp +++ b/test/linalg/test_linalg_determinant.fypp @@ -1,47 +1,42 @@ #:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES ! Test matrix determinant module test_linalg_determinant - use stdlib_linalg_interface + 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(error) - logical, intent(out) :: error - - real :: t0,t1 - - call cpu_time(t0) - - #:for rk,rt,ri in ALL_KINDS_TYPES - call test_${ri}$_eye_determinant(error) - if (error) return - - call test_${ri}$_eye_multiple(error) - if (error) return - #: endfor - - #:for ck,ct,ci in CMPL_KINDS_TYPES - call test_${ci}$_complex_determinant(error) - if (error) return + 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 + 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)] + #:endfor + #:for ck,ct in CMPLX_KINDS_TYPES + tests = [tests,new_unittest("$complex_det_${rt[0]}$${rk}$",test_${ct[0]}$${ck}$_complex_determinant)] #: endfor - call cpu_time(t1) - - print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) - - 1 format('Determinant tests completed in ',f9.4,' milliseconds, result=',a) - end subroutine test_matrix_determinant !> Determinant of identity matrix - #:for rk,rt,ri in ALL_KINDS_TYPES - subroutine test_${ri}$_eye_determinant(error) - logical, intent(out) :: error + #:for rk,rt in RC_KINDS_TYPES + subroutine test_${rt[0]}$${rk}$_eye_determinant(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: i integer(ilp), parameter :: n = 128_ilp @@ -52,16 +47,19 @@ module test_linalg_determinant !> Determinant function deta = det(a,err=state) - error = state%error() .or. .not.abs(deta-1.0_${rk}$) Determinant of identity matrix multiplier - subroutine test_${ri}$_eye_multiple(error) - logical, intent(out) :: error + subroutine test_${rt[0]}$${rk}$_eye_multiple(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp), parameter :: n = 10_ilp real(${rk}$), parameter :: coef = 0.0001_${rk}$ @@ -76,19 +74,22 @@ module test_linalg_determinant !> Determinant: small, but a is not singular, because it is a multiple of the identity. deta = det(a,err=state) - error = state%error() .or. .not.abs(deta-coef**n) Determinant of complex identity matrix - #:for ck,ct,ci in CMPL_KINDS_TYPES - subroutine test_${ci}$_complex_determinant(error) - logical, intent(out) :: error + #:for ck,ct in CMPLX_KINDS_TYPES + subroutine test_${ct[0]}$${ck}$_complex_determinant(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: i,j,n integer(ilp), parameter :: nmax = 10_ilp @@ -116,12 +117,40 @@ module test_linalg_determinant end do matrix_size - error = state%error() .or. any(.not.abs(res-deta)<=tiny(0.0_${ck}$)) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(res-deta)<=tiny(0.0_${ck}$)), & + 'det((1+i)*eye(n)) does not match result') - end subroutine test_${ci}$_complex_determinant + end subroutine test_${ct[0]}$${ck}$_complex_determinant #: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 From 4a29219d4d8220f25b59f35bb57db20cbb48b41d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 11:21:26 +0200 Subject: [PATCH 05/28] make `pure` --- src/stdlib_linalg_determinant.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index bce5c248b..aba32792c 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -39,9 +39,9 @@ module stdlib_linalg_determinant #:for rk,rt in RC_KINDS_TYPES ! Compute determinant of a square matrix A: pure interface - function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) + pure function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !> Input matrix a[m,n] - ${rt}$, intent(in), target :: a(:,:) + ${rt}$, intent(in) :: a(:,:) !> Result: matrix determinant ${rt}$ :: det From 3d05cc38e594a110cb935f7e5e79c130eb5a6523 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 12:34:37 +0200 Subject: [PATCH 06/28] remove `xdp` --- src/stdlib_linalg_determinant.fypp | 6 ++++++ test/linalg/test_linalg_determinant.fypp | 15 +++++++++++---- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index aba32792c..9e544d143 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -21,23 +21,28 @@ module stdlib_linalg_determinant interface det #:for rk,rt in RC_KINDS_TYPES + #:if rk!="xdp" ! Interface with error control module procedure stdlib_linalg_${rt[0]}$${rk}$determinant ! Pure interface module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant + #:endif #:endfor end interface det ! Pure Operator interface interface operator(.det.) #: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" ! Compute determinant of a square matrix A: pure interface pure function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !> Input matrix a[m,n] @@ -210,6 +215,7 @@ module stdlib_linalg_determinant end function stdlib_linalg_${rt[0]}$${rk}$determinant + #:endif #:endfor end module stdlib_linalg_determinant diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp index 7c3749e2c..d0a323d80 100644 --- a/test/linalg/test_linalg_determinant.fypp +++ b/test/linalg/test_linalg_determinant.fypp @@ -22,17 +22,22 @@ module test_linalg_determinant 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 !> Determinant of identity matrix #:for rk,rt in RC_KINDS_TYPES + #:if rk!="xdp" subroutine test_${rt[0]}$${rk}$_eye_determinant(error) type(error_type), allocatable, intent(out) :: error @@ -61,8 +66,8 @@ module test_linalg_determinant type(linalg_state_type) :: state - integer(ilp), parameter :: n = 10_ilp - real(${rk}$), parameter :: coef = 0.0001_${rk}$ + integer(ilp), parameter :: n = 4_ilp + real(${rk}$), parameter :: coef = 0.01_${rk}$ integer(ilp) :: i,j ${rt}$ :: a(n,n),deta @@ -78,14 +83,15 @@ module test_linalg_determinant if (allocated(error)) return call check(error, abs(deta-coef**n) 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 @@ -125,6 +131,7 @@ module test_linalg_determinant end subroutine test_${ct[0]}$${ck}$_complex_determinant + #:endif #:endfor end module test_linalg_determinant From 4beab1f004a4dcffc9dd75c3ade275f3a12fb171 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 13 Apr 2024 14:42:14 +0200 Subject: [PATCH 07/28] Documentation Documentation --- src/stdlib_linalg_determinant.fypp | 121 ++++++++++++++++++++++++----- 1 file changed, 101 insertions(+), 20 deletions(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index 9e544d143..cb5542f8c 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -!> Determinant of a rectangular matrix 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, & @@ -9,29 +9,69 @@ module stdlib_linalg_determinant implicit none(type,external) private - !> Function interface + ! Function interface public :: det public :: operator(.det.) character(*), parameter :: this = 'determinant' - ! Numpy: det(a) - ! Scipy: det(a, overwrite_a=False, check_finite=True) - ! IMSL: DET(a) - 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" - ! Interface with error control - module procedure stdlib_linalg_${rt[0]}$${rk}$determinant - ! Pure interface + module procedure stdlib_linalg_${rt[0]}$${rk}$determinant module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endif #:endfor end interface det - ! Pure Operator interface 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 @@ -43,20 +83,39 @@ module stdlib_linalg_determinant #:for rk,rt in RC_KINDS_TYPES #:if rk!="xdp" - ! Compute determinant of a square matrix A: pure interface - pure function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) + 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(:,:) - !> Result: matrix determinant + !> Matrix determinant ${rt}$ :: det - !> Local variables + !! Local variables type(linalg_state_type) :: err0 integer(ilp) :: m,n,info,perm,k integer(ilp), allocatable :: ipiv(:) ${rt}$, allocatable :: amat(:,:) - !> Matrix determinant size + ! Matrix determinant size m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) @@ -121,25 +180,47 @@ module stdlib_linalg_determinant end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant - ! Compute determinant of a square matrix A, with error control 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 - !> [optional] state return flag. On error if not requested, the code will stop + !> State return flag. type(linalg_state_type), intent(out) :: err - !> Result: matrix determinant + !> Matrix determinant ${rt}$ :: det - !> Local variables + !! 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 + ! Matrix determinant size m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) From 15023e1ba1fae9a61f6aeafc25db047dcef327fe Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 14 Apr 2024 10:13:23 +0200 Subject: [PATCH 08/28] `submodule version` --- src/stdlib_linalg.fypp | 95 +++++++++++++++++++++++++++++- src/stdlib_linalg_determinant.fypp | 76 ++---------------------- 2 files changed, 96 insertions(+), 75 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 4047612c8..d1428b613 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1,18 +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 - use stdlib_linalg_determinant, only: det implicit none private public :: det + public :: operator(.det.) public :: diag public :: eye public :: trace @@ -220,6 +221,94 @@ module stdlib_linalg #:endfor end interface is_hessenberg + + 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.) + + 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 + 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 index cb5542f8c..be1d4b79f 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -1,89 +1,21 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module stdlib_linalg_determinant +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) - 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) + module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !!### Summary !! Compute determinant of a real square matrix (pure interface). !! @@ -180,7 +112,7 @@ module stdlib_linalg_determinant end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant - function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det) + module function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det) !!### Summary !! Compute determinant of a square matrix (with error control). !! @@ -299,4 +231,4 @@ module stdlib_linalg_determinant #:endif #:endfor -end module stdlib_linalg_determinant +end submodule stdlib_linalg_determinant From 5923a546b011a910ba3e713a94a7dee88655b380 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 08:56:56 +0200 Subject: [PATCH 09/28] Update src/stdlib_linalg.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d1428b613..d6be3ff3d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -231,7 +231,7 @@ module stdlib_linalg !! 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. + !!@note The provided functions are intended for square matrices only. !! !!### Example !! From 31cdc84a843faf09b0902b19a8ac34db58878246 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 08:57:08 +0200 Subject: [PATCH 10/28] Update src/stdlib_linalg.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d6be3ff3d..c856a020f 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -229,7 +229,7 @@ module stdlib_linalg !!### Description !! !! This interface provides methods for computing the determinant of a matrix. - !! Supported data types include real and complex. + !! Supported data types include `real` and `complex`. !! !!@note The provided functions are intended for square matrices only. !! From 45a606fc638be5b36db03b12bcce5f926ab7612c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:05:00 +0200 Subject: [PATCH 11/28] add `det` example --- example/linalg/CMakeLists.txt | 2 ++ example/linalg/example_determinant.f90 | 14 ++++++++++++++ 2 files changed, 16 insertions(+) create mode 100644 example/linalg/example_determinant.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1a5875502..cf6bd75da 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -18,3 +18,5 @@ ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) +ADD_EXAMPLE(determinant) + 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 From cacb585ccca9e62067aee6b6a9e4f3ff0a8d9553 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:07:17 +0200 Subject: [PATCH 12/28] Update src/stdlib_linalg_determinant.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg_determinant.fypp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index be1d4b79f..efe4eceaa 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -54,7 +54,9 @@ submodule (stdlib_linalg) stdlib_linalg_determinant 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 + ! Process output and return + call linalg_error_handling(err0) + return end if select case (m) From 5c16ff8a663704cb88f70a7abed653ddc8f9b518 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:07:24 +0200 Subject: [PATCH 13/28] Update src/stdlib_linalg_determinant.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg_determinant.fypp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index efe4eceaa..664cf7dca 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -161,7 +161,9 @@ submodule (stdlib_linalg) stdlib_linalg_determinant 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 + ! Process output and return + call linalg_error_handling(err0) + return end if ! Can A be overwritten? By default, do not overwrite From ab030c547642654e37ab69e4dfe16ed22635b2ad Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:09:26 +0200 Subject: [PATCH 14/28] Update src/stdlib_linalg_determinant.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg_determinant.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index 664cf7dca..14f98199a 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -190,7 +190,7 @@ submodule (stdlib_linalg) stdlib_linalg_determinant ! Initialize a matrix temporary if (copy_a) then - allocate(amat(m,n),source=a) + allocate(amat, source=a) else amat => a endif From 13bd98a1b33a944d994f93418414a565c090ee1f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:10:32 +0200 Subject: [PATCH 15/28] warn about `xdp` --- src/stdlib_linalg.fypp | 6 +++--- src/stdlib_linalg_determinant.fypp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d1428b613..9113f9381 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -251,7 +251,7 @@ module stdlib_linalg !!``` !! #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" + #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp module procedure stdlib_linalg_${rt[0]}$${rk}$determinant module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endif @@ -280,7 +280,7 @@ module stdlib_linalg !!``` ! #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" + #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endif #:endfor @@ -288,7 +288,7 @@ module stdlib_linalg interface #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" + #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support 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(:,:) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index be1d4b79f..4efb97ea7 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -14,7 +14,7 @@ submodule (stdlib_linalg) stdlib_linalg_determinant contains #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" + #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !!### Summary !! Compute determinant of a real square matrix (pure interface). From 504d90d8132136246e88bef23240ad6f8b2b92a0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 09:14:55 +0200 Subject: [PATCH 16/28] cleanup xdp notes --- src/stdlib_linalg.fypp | 8 +++++--- src/stdlib_linalg_determinant.fypp | 17 +++++++++-------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index e4081de06..732e21f11 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -232,6 +232,7 @@ module stdlib_linalg !! 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 !! @@ -251,7 +252,7 @@ module stdlib_linalg !!``` !! #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp + #:if rk!="xdp" module procedure stdlib_linalg_${rt[0]}$${rk}$determinant module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endif @@ -268,6 +269,7 @@ module stdlib_linalg !! 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 !! @@ -280,7 +282,7 @@ module stdlib_linalg !!``` ! #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp + #:if rk!="xdp" module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant #:endif #:endfor @@ -288,7 +290,7 @@ module stdlib_linalg interface #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp + #: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(:,:) diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index a67d1eef1..2420c12c3 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -13,8 +13,9 @@ submodule (stdlib_linalg) stdlib_linalg_determinant contains + ! BLAS/LAPACK backends do not currently support xdp #:for rk,rt in RC_KINDS_TYPES - #:if rk!="xdp" ! BLAS/LAPACK backends do not currently support xdp + #:if rk!="xdp" module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !!### Summary !! Compute determinant of a real square matrix (pure interface). @@ -54,8 +55,8 @@ submodule (stdlib_linalg) stdlib_linalg_determinant 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) + ! Process output and return + call linalg_error_handling(err0) return end if @@ -110,7 +111,7 @@ submodule (stdlib_linalg) stdlib_linalg_determinant end select ! Process output and return - 1 call linalg_error_handling(err0) + call linalg_error_handling(err0) end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant @@ -161,9 +162,9 @@ submodule (stdlib_linalg) stdlib_linalg_determinant 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 + ! Process output and return + call linalg_error_handling(err0,err) + return end if ! Can A be overwritten? By default, do not overwrite @@ -228,7 +229,7 @@ submodule (stdlib_linalg) stdlib_linalg_determinant end select ! Process output and return - 1 call linalg_error_handling(err0,err) + call linalg_error_handling(err0,err) end function stdlib_linalg_${rt[0]}$${rk}$determinant From e80b50824305d07a4bebe1c60adde36f05388e58 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 10:27:00 +0200 Subject: [PATCH 17/28] spacing --- src/stdlib_linalg.fypp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 732e21f11..d85a80978 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -241,14 +241,15 @@ module stdlib_linalg !! 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 @@ -275,10 +276,12 @@ module stdlib_linalg !! !!```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 From c6076eaf8f9d4f3b9c1fe63824037652b1c5741f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 10:59:28 +0200 Subject: [PATCH 18/28] relax error thresholds --- test/linalg/test_linalg_determinant.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp index d0a323d80..267d395e9 100644 --- a/test/linalg/test_linalg_determinant.fypp +++ b/test/linalg/test_linalg_determinant.fypp @@ -56,7 +56,7 @@ module test_linalg_determinant call check(error,state%ok(),state%print()) if (allocated(error)) return - call check(error, abs(deta-1.0_${rk}$) Date: Mon, 15 Apr 2024 11:01:49 +0200 Subject: [PATCH 19/28] restore `pure` attribute --- src/stdlib_linalg.fypp | 2 +- src/stdlib_linalg_determinant.fypp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d85a80978..7ca72e035 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -304,7 +304,7 @@ module stdlib_linalg !> Matrix determinant ${rt}$ :: det end function stdlib_linalg_${rt[0]}$${rk}$determinant - module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) + pure module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !> Input matrix a[m,n] ${rt}$, intent(in) :: a(:,:) !> Matrix determinant diff --git a/src/stdlib_linalg_determinant.fypp b/src/stdlib_linalg_determinant.fypp index 2420c12c3..0091d596b 100644 --- a/src/stdlib_linalg_determinant.fypp +++ b/src/stdlib_linalg_determinant.fypp @@ -16,7 +16,7 @@ submodule (stdlib_linalg) stdlib_linalg_determinant ! BLAS/LAPACK backends do not currently support xdp #:for rk,rt in RC_KINDS_TYPES #:if rk!="xdp" - module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) + pure module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det) !!### Summary !! Compute determinant of a real square matrix (pure interface). !! From 5d52d4895b4366afd50273e4b9db0a5d795841fa Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 15 Apr 2024 18:33:25 +0200 Subject: [PATCH 20/28] Update test/linalg/test_linalg_determinant.fypp Co-authored-by: Jeremie Vandenplas --- test/linalg/test_linalg_determinant.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_determinant.fypp b/test/linalg/test_linalg_determinant.fypp index 267d395e9..6b9310f72 100644 --- a/test/linalg/test_linalg_determinant.fypp +++ b/test/linalg/test_linalg_determinant.fypp @@ -35,9 +35,9 @@ module test_linalg_determinant end subroutine test_matrix_determinant - !> Determinant of identity matrix #: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 From cff995d60405dd814b32266e74267e8080572008 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 15:08:06 +0200 Subject: [PATCH 21/28] add docs --- doc/specs/stdlib_linalg.md | 73 +++++++++++++++++++++++++ example/linalg/CMakeLists.txt | 1 + example/linalg/example_determinant2.f90 | 13 +++++ 3 files changed, 87 insertions(+) create mode 100644 example/linalg/example_determinant2.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ab52dfb71..e7ef2af5d 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -595,3 +595,76 @@ 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 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. + +`err` (optional): Shall be a `type(linalg_state_type)` return value. + +### 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 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 + +### 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 cf6bd75da..f515af446 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -19,4 +19,5 @@ ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(determinant) +ADD_EXAMPLE(determinant2) 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 From 2fe2428ea9ae50f601085a1d32525771863586d5 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 15:11:06 +0200 Subject: [PATCH 22/28] link to specs --- src/stdlib_linalg.fypp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 7ca72e035..2f837f15a 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -223,6 +223,11 @@ module stdlib_linalg 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. !! @@ -261,6 +266,11 @@ module stdlib_linalg 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. !! From 45447a89026926e108ad0a4df0d7ab2b2f05c26a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:17:32 +0200 Subject: [PATCH 23/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index e7ef2af5d..c2978406f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -611,7 +611,7 @@ allows for more expert control. ### Syntax -`c = ` [[stdlib_linalg(module):det(interface)]] `(a,overwrite_a,err)` +`c = ` [[stdlib_linalg(module):det(interface)]] `(a [, overwrite_a, err])` ### Arguments From 3ee20ad11933e70397ba9318c33b190c94996313 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:18:08 +0200 Subject: [PATCH 24/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index c2978406f..f49b9ad18 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -617,7 +617,8 @@ allows for more expert control. `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. +`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)` return value. From 1e50115acff1122beeeb036f16145c0d43519279 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:18:16 +0200 Subject: [PATCH 25/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index f49b9ad18..4b5070309 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -620,7 +620,7 @@ allows for more expert control. `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)` return value. +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value From fe933ad9e6cf91a1616418b4a549819121f813bf Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:19:27 +0200 Subject: [PATCH 26/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 4b5070309..353f7d31e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -604,7 +604,7 @@ Experimental ### Description -This function computes the determinant of a real square matrix. +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. From 59dae894392a3b2c081bf2ed98e826c4d0e8a8ed Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:19:53 +0200 Subject: [PATCH 27/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 353f7d31e..69339d071 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -624,7 +624,7 @@ allows for more expert control. ### Return value -Returns a real scalar value that represents the determinnt of the matrix. +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. From 3c72b06a3cd70a0657f3763964811dbaf5287e2c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:20:25 +0200 Subject: [PATCH 28/28] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 69339d071..8a42ae25a 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -654,7 +654,7 @@ This interface is equivalent to the `pure` version of determinant [[stdlib_linal ### Arguments -`a`: Shall be a rank-2 square array +`a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument. ### Return value