Skip to content

Commit d54cfa3

Browse files
perazzperazz
committed
add pure interface
1 parent 4a4b899 commit d54cfa3

File tree

2 files changed

+108
-21
lines changed

2 files changed

+108
-21
lines changed

‎src/stdlib_linalg.fypp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,12 @@ module stdlib_linalg
77
int8, int16, int32, int64
88
use stdlib_error, only: error_stop
99
use stdlib_optval, only: optval
10+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
11+
use stdlib_linalg_determinant, only: det
1012
implicit none
1113
private
1214

15+
public :: det
1316
public :: diag
1417
public :: eye
1518
public :: trace
@@ -23,6 +26,9 @@ module stdlib_linalg
2326
public :: is_hermitian
2427
public :: is_triangular
2528
public :: is_hessenberg
29+
30+
! Export linalg error handling
31+
public :: linalg_state_type, linalg_error_handling
2632

2733
interface diag
2834
!! version: experimental

‎src/stdlib_linalg_determinant.fypp‎

Lines changed: 102 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,15 @@
11
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
!> Determinant of a rectangular matrix
24
module stdlib_linalg_determinant
35
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
6+
use stdlib_linalg_lapack, only: getrf
7+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
8+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
89
implicit none(type,external)
910
private
1011

11-
!> Determinant of a rectangular matrix
12+
!> Function interface
1213
public :: det
1314

1415
character(*), parameter :: this = 'determinant'
@@ -18,28 +19,106 @@ module stdlib_linalg_determinant
1819
! IMSL: DET(a)
1920

2021
interface det
21-
#:for rk,rt,ri in ALL_KINDS_TYPES
22-
module procedure stdlib_linalg_${ri}$determinant
22+
#:for rk,rt in RC_KINDS_TYPES
23+
module procedure stdlib_linalg_${rt[0]}$${rk}$determinant
24+
module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
2325
#:endfor
2426
end interface det
2527

26-
2728
contains
2829

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)
30+
#:for rk,rt in RC_KINDS_TYPES
31+
! Compute determinant of a square matrix A: pure interface
32+
function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det)
33+
!> Input matrix a[m,n]
34+
${rt}$, intent(in), target :: a(:,:)
35+
!> Result: matrix determinant
36+
${rt}$ :: det
37+
38+
!> Local variables
39+
type(linalg_state_type) :: err0
40+
integer(ilp) :: m,n,info,perm,k
41+
integer(ilp), allocatable :: ipiv(:)
42+
${rt}$, allocatable :: amat(:,:)
43+
44+
!> Matrix determinant size
45+
m = size(a,1,kind=ilp)
46+
n = size(a,2,kind=ilp)
47+
48+
if (m/=n .or. .not.min(m,n)>=0) then
49+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']')
50+
det = 0.0_${rk}$
51+
goto 1
52+
end if
53+
54+
select case (m)
55+
case (0)
56+
57+
! Empty array has determinant 1 because math
58+
det = 1.0_${rk}$
59+
60+
case (1)
61+
62+
! Scalar input
63+
det = a(1,1)
64+
65+
case default
66+
67+
! Find determinant from LU decomposition
68+
69+
! Initialize a matrix temporary
70+
allocate(amat(m,n),source=a)
71+
72+
! Pivot indices
73+
allocate(ipiv(n))
74+
75+
! Compute determinant from LU factorization, then calculate the
76+
! product of all diagonal entries of the U factor.
77+
call getrf(m,n,amat,m,ipiv,info)
78+
79+
select case (info)
80+
case (0)
81+
! Success: compute determinant
82+
83+
! Start with real 1.0
84+
det = 1.0_${rk}$
85+
perm = 0
86+
do k=1,n
87+
if (ipiv(k)/=k) perm = perm+1
88+
det = det*amat(k,k)
89+
end do
90+
if (mod(perm,2)/=0) det = -det
91+
92+
case (:-1)
93+
err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
94+
case (1:)
95+
err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
96+
case default
97+
err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
98+
end select
99+
100+
deallocate(amat)
101+
102+
end select
103+
104+
! Process output and return
105+
1 call linalg_error_handling(err0)
106+
107+
end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
108+
109+
! Compute determinant of a square matrix A, with error control
110+
function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det)
32111
!> Input matrix a[m,n]
33112
${rt}$, intent(inout), target :: a(:,:)
34113
!> [optional] Can A data be overwritten and destroyed?
35114
logical(lk), optional, intent(in) :: overwrite_a
36115
!> [optional] state return flag. On error if not requested, the code will stop
37-
type(linalg_state), optional, intent(out) :: err
116+
type(linalg_state_type), intent(out) :: err
38117
!> Result: matrix determinant
39118
${rt}$ :: det
40119

41120
!> Local variables
42-
type(linalg_state) :: err0
121+
type(linalg_state_type) :: err0
43122
integer(ilp) :: m,n,info,perm,k
44123
integer(ilp), allocatable :: ipiv(:)
45124
logical(lk) :: copy_a
@@ -50,7 +129,7 @@ module stdlib_linalg_determinant
50129
n = size(a,2,kind=ilp)
51130

52131
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,']')
132+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or non-square matrix: a=[',m,',',n,']')
54133
det = 0.0_${rk}$
55134
goto 1
56135
end if
@@ -64,11 +143,13 @@ module stdlib_linalg_determinant
64143

65144
select case (m)
66145
case (0)
146+
67147
! Empty array has determinant 1 because math
68148
det = 1.0_${rk}$
69149

70150
case (1)
71-
! Scalar
151+
152+
! Scalar input
72153
det = a(1,1)
73154

74155
case default
@@ -85,8 +166,8 @@ module stdlib_linalg_determinant
85166
! Pivot indices
86167
allocate(ipiv(n))
87168

88-
! Compute determinant from LU factorization, then calculate the product of
89-
! all diagonal entries of the U factor.
169+
! Compute determinant from LU factorization, then calculate the
170+
! product of all diagonal entries of the U factor.
90171
call getrf(m,n,amat,m,ipiv,info)
91172

92173
select case (info)
@@ -103,11 +184,11 @@ module stdlib_linalg_determinant
103184
if (mod(perm,2)/=0) det = -det
104185

105186
case (:-1)
106-
err0 = linalg_state(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
187+
err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',m,',',n,']')
107188
case (1:)
108-
err0 = linalg_state(this,LINALG_ERROR,'singular matrix')
189+
err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
109190
case default
110-
err0 = linalg_state(this,LINALG_INTERNAL_ERROR,'catastrophic error')
191+
err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
111192
end select
112193

113194
if (.not.copy_a) deallocate(amat)
@@ -117,7 +198,7 @@ module stdlib_linalg_determinant
117198
! Process output and return
118199
1 call linalg_error_handling(err0,err)
119200

120-
end function stdlib_linalg_${ri}$determinant
201+
end function stdlib_linalg_${rt[0]}$${rk}$determinant
121202

122203
#:endfor
123204

0 commit comments

Comments
 (0)