Skip to content

Commit 4a4b899

Browse files
perazzperazz
committed
import files
1 parent c26253a commit 4a4b899

File tree

4 files changed

+254
-0
lines changed

4 files changed

+254
-0
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ set(fppFiles
2424
stdlib_linalg_outer_product.fypp
2525
stdlib_linalg_kronecker.fypp
2626
stdlib_linalg_cross_product.fypp
27+
stdlib_linalg_determinant.fypp
2728
stdlib_linalg_state.fypp
2829
stdlib_optval.fypp
2930
stdlib_selection.fypp

‎src/stdlib_linalg_determinant.fypp

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
#:include "common.fypp"
2+
module stdlib_linalg_determinant
3+
use stdlib_linalg_constants
4+
use stdlib_linalg_blas
5+
use stdlib_linalg_lapack
6+
use stdlib_linalg_state
7+
use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit
8+
implicit none(type,external)
9+
private
10+
11+
!> Determinant of a rectangular matrix
12+
public :: det
13+
14+
character(*), parameter :: this = 'determinant'
15+
16+
! Numpy: det(a)
17+
! Scipy: det(a, overwrite_a=False, check_finite=True)
18+
! IMSL: DET(a)
19+
20+
interface det
21+
#:for rk,rt,ri in ALL_KINDS_TYPES
22+
module procedure stdlib_linalg_${ri}$determinant
23+
#:endfor
24+
end interface det
25+
26+
27+
contains
28+
29+
#:for rk,rt,ri in ALL_KINDS_TYPES
30+
! Compute determinant of a square matrix A
31+
function stdlib_linalg_${ri}$determinant(a,overwrite_a,err) result(det)
32+
!> Input matrix a[m,n]
33+
${rt}$, intent(inout), target :: a(:,:)
34+
!> [optional] Can A data be overwritten and destroyed?
35+
logical(lk), optional, intent(in) :: overwrite_a
36+
!> [optional] state return flag. On error if not requested, the code will stop
37+
type(linalg_state), optional, intent(out) :: err
38+
!> Result: matrix determinant
39+
${rt}$ :: det
40+
41+
!> Local variables
42+
type(linalg_state) :: err0
43+
integer(ilp) :: m,n,info,perm,k
44+
integer(ilp), allocatable :: ipiv(:)
45+
logical(lk) :: copy_a
46+
${rt}$, pointer :: amat(:,:)
47+
48+
!> Matrix determinant size
49+
m = size(a,1,kind=ilp)
50+
n = size(a,2,kind=ilp)
51+
52+
if (m/=n .or. .not.min(m,n)>=0) then
53+
err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']')
54+
det = 0.0_${rk}$
55+
goto 1
56+
end if
57+
58+
! Can A be overwritten? By default, do not overwrite
59+
if (present(overwrite_a)) then
60+
copy_a = .not.overwrite_a
61+
else
62+
copy_a = .true._lk
63+
endif
64+
65+
select case (m)
66+
case (0)
67+
! Empty array has determinant 1 because math
68+
det = 1.0_${rk}$
69+
70+
case (1)
71+
! Scalar
72+
det = a(1,1)
73+
74+
case default
75+
76+
! Find determinant from LU decomposition
77+
78+
! Initialize a matrix temporary
79+
if (copy_a) then
80+
allocate(amat(m,n),source=a)
81+
else
82+
amat => a
83+
endif
84+
85+
! Pivot indices
86+
allocate(ipiv(n))
87+
88+
! Compute determinant from LU factorization, then calculate the product of
89+
! all diagonal entries of the U factor.
90+
call getrf(m,n,amat,m,ipiv,info)
91+
92+
select case (info)
93+
case (0)
94+
! Success: compute determinant
95+
96+
! Start with real 1.0
97+
det = 1.0_${rk}$
98+
perm = 0
99+
do k=1,n
100+
if (ipiv(k)/=k) perm = perm+1
101+
det = det*amat(k,k)
102+
end do
103+
if (mod(perm,2)/=0) det = -det
104+
105+
case (:-1)
106+
err0 = linalg_state(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
107+
case (1:)
108+
err0 = linalg_state(this,LINALG_ERROR,'singular matrix')
109+
case default
110+
err0 = linalg_state(this,LINALG_INTERNAL_ERROR,'catastrophic error')
111+
end select
112+
113+
if (.not.copy_a) deallocate(amat)
114+
115+
end select
116+
117+
! Process output and return
118+
1 call linalg_error_handling(err0,err)
119+
120+
end function stdlib_linalg_${ri}$determinant
121+
122+
#:endfor
123+
124+
end module stdlib_linalg_determinant

‎test/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,12 @@ set(
22
fppFiles
33
"test_linalg.fypp"
44
"test_blas_lapack.fypp"
5+
"test_linalg_determinant.fypp"
56
"test_linalg_matrix_property_checks.fypp"
67
)
78
fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
89

910
ADDTEST(linalg)
1011
ADDTEST(linalg_matrix_property_checks)
1112
ADDTEST(blas_lapack)
13+
ADDTEST(linalg_determinant)

‎test/linalg/test_linalg_determinant.fypp

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
#:include "common.fypp"
2+
! Test matrix determinant
3+
module test_linalg_determinant
4+
use stdlib_linalg_interface
5+
6+
implicit none (type,external)
7+
8+
contains
9+
10+
!> Matrix inversion tests
11+
subroutine test_matrix_determinant(error)
12+
logical, intent(out) :: error
13+
14+
real :: t0,t1
15+
16+
call cpu_time(t0)
17+
18+
#:for rk,rt,ri in ALL_KINDS_TYPES
19+
call test_${ri}$_eye_determinant(error)
20+
if (error) return
21+
22+
call test_${ri}$_eye_multiple(error)
23+
if (error) return
24+
#: endfor
25+
26+
#:for ck,ct,ci in CMPL_KINDS_TYPES
27+
call test_${ci}$_complex_determinant(error)
28+
if (error) return
29+
#: endfor
30+
31+
call cpu_time(t1)
32+
33+
print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error)
34+
35+
1 format('Determinant tests completed in ',f9.4,' milliseconds, result=',a)
36+
37+
end subroutine test_matrix_determinant
38+
39+
!> Determinant of identity matrix
40+
#:for rk,rt,ri in ALL_KINDS_TYPES
41+
subroutine test_${ri}$_eye_determinant(error)
42+
logical, intent(out) :: error
43+
44+
type(linalg_state) :: state
45+
46+
integer(ilp) :: i
47+
integer(ilp), parameter :: n = 128_ilp
48+
49+
${rt}$ :: a(n,n),deta
50+
51+
a = eye(n)
52+
53+
!> Determinant function
54+
deta = det(a,err=state)
55+
error = state%error() .or. .not.abs(deta-1.0_${rk}$)<tiny(0.0_${rk}$)
56+
if (error) return
57+
58+
end subroutine test_${ri}$_eye_determinant
59+
60+
!> Determinant of identity matrix multiplier
61+
subroutine test_${ri}$_eye_multiple(error)
62+
logical, intent(out) :: error
63+
64+
type(linalg_state) :: state
65+
66+
integer(ilp), parameter :: n = 10_ilp
67+
real(${rk}$), parameter :: coef = 0.0001_${rk}$
68+
integer(ilp) :: i,j
69+
${rt}$ :: a(n,n),deta
70+
71+
!> Multiply eye by a very small number
72+
a = eye(n)
73+
do concurrent (i=1:n)
74+
a(i,i) = coef
75+
end do
76+
77+
!> Determinant: small, but a is not singular, because it is a multiple of the identity.
78+
deta = det(a,err=state)
79+
error = state%error() .or. .not.abs(deta-coef**n)<max(tiny(0.0_${rk}$),epsilon(0.0_${rk}$)*coef**n)
80+
if (error) return
81+
82+
end subroutine test_${ri}$_eye_multiple
83+
84+
#:endfor
85+
86+
!> Determinant of complex identity matrix
87+
#:for ck,ct,ci in CMPL_KINDS_TYPES
88+
subroutine test_${ci}$_complex_determinant(error)
89+
logical, intent(out) :: error
90+
91+
type(linalg_state) :: state
92+
93+
integer(ilp) :: i,j,n
94+
integer(ilp), parameter :: nmax = 10_ilp
95+
96+
${ct}$, parameter :: res(nmax) = [${ct}$::(1,1),(0,2),(-2,2),(-4,0),(-4,-4), &
97+
(0,-8),(8,-8),(16,0),(16,16),(0,32)]
98+
99+
${ct}$, allocatable :: a(:,:)
100+
${ct}$ :: deta(nmax)
101+
102+
!> Test determinant for all sizes, 1:nmax
103+
matrix_size: do n=1,nmax
104+
105+
! Put 1+i on each diagonal element
106+
a = eye(n)
107+
do concurrent (i=1:n)
108+
a(i,i) = (1.0_${ck}$,1.0_${ck}$)
109+
end do
110+
111+
! Expected result
112+
deta(n) = det(a,err=state)
113+
114+
deallocate(a)
115+
if (state%error()) exit matrix_size
116+
117+
end do matrix_size
118+
119+
error = state%error() .or. any(.not.abs(res-deta)<=tiny(0.0_${ck}$))
120+
121+
end subroutine test_${ci}$_complex_determinant
122+
123+
#:endfor
124+
125+
end module test_linalg_determinant
126+
127+

0 commit comments

Comments
 (0)