Skip to content

Commit 5c817c8

Browse files
perazzperazz
committed
pure interfaces
1 parent d929077 commit 5c817c8

File tree

1 file changed

+91
-27
lines changed

1 file changed

+91
-27
lines changed

‎src/stdlib_linalg_solve.fypp

Lines changed: 91 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,6 @@ module stdlib_linalg_solve
1010
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
1111
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
1212
implicit none(type,external)
13-
private
1413

1514
!> Solve a linear system
1615
public :: solve
@@ -24,13 +23,40 @@ module stdlib_linalg_solve
2423
#:for rk,rt,ri in RC_KINDS_TYPES
2524
#:if rk!="xdp"
2625
module procedure stdlib_linalg_${ri}$solve${ndsuf}$
26+
module procedure stdlib_linalg_${ri}$_pure_solve${ndsuf}$
2727
#:endif
2828
#:endfor
2929
#:endfor
3030
end interface solve
31-
31+
32+
33+
character(*), parameter :: this = 'solve'
3234

3335
contains
36+
37+
elemental subroutine handle_gesv_info(info,lda,n,nrhs,err)
38+
integer(ilp), intent(in) :: info,lda,n,nrhs
39+
type(linalg_state_type), intent(out) :: err
40+
41+
! Process output
42+
select case (info)
43+
case (0)
44+
! Success
45+
case (-1)
46+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
47+
case (-2)
48+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
49+
case (-4)
50+
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n])
51+
case (-7)
52+
err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
53+
case (1:)
54+
err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
55+
case default
56+
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
57+
end select
58+
59+
end subroutine handle_gesv_info
3460

3561
#:for nd,ndsuf,nde in ALL_RHS
3662
#:for rk,rt,ri in RC_KINDS_TYPES
@@ -44,7 +70,7 @@ module stdlib_linalg_solve
4470
!> [optional] Can A data be overwritten and destroyed?
4571
logical(lk), optional, intent(in) :: overwrite_a
4672
!> [optional] state return flag. On error if not requested, the code will stop
47-
type(linalg_state_type), optional, intent(out) :: err
73+
type(linalg_state_type), intent(out) :: err
4874
!> Result array/matrix x[n] or x[n,nrhs]
4975
${rt}$, allocatable, target :: x${nd}$
5076

@@ -53,20 +79,20 @@ module stdlib_linalg_solve
5379
integer(ilp) :: lda,n,ldb,nrhs,info
5480
integer(ilp), allocatable :: ipiv(:)
5581
logical(lk) :: copy_a
56-
${rt}$, pointer :: xmat(:,:),amat(:,:)
57-
character(*), parameter :: this = 'solve'
82+
${rt}$, pointer :: xmat(:,:),amat(:,:)
5883

5984
!> Problem sizes
6085
lda = size(a,1,kind=ilp)
6186
n = size(a,2,kind=ilp)
6287
ldb = size(b,1,kind=ilp)
6388
nrhs = size(b ,kind=ilp)/ldb
6489

65-
if (lda<1 .or. n<1 .or. ldb<1 .or. lda/=n .or. ldb/=n) then
66-
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=[',lda,',',n,'],',&
67-
'b=[',ldb,',',nrhs,']')
90+
if (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then
91+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], &
92+
', b=',[ldb,nrhs])
6893
allocate(x${nde}$)
69-
goto 1
94+
call linalg_error_handling(err0,err)
95+
return
7096
end if
7197

7298
! Can A be overwritten? By default, do not overwrite
@@ -94,30 +120,68 @@ module stdlib_linalg_solve
94120
call gesv(n,nrhs,amat,lda,ipiv,xmat,ldb,info)
95121

96122
! Process output
97-
select case (info)
98-
case (0)
99-
! Success
100-
case (-1)
101-
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
102-
case (-2)
103-
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
104-
case (-4)
105-
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=[',lda,',',n,']')
106-
case (-7)
107-
err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']')
108-
case (1:)
109-
err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
110-
case default
111-
err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
112-
end select
123+
call handle_gesv_info(info,lda,n,nrhs,err0)
113124

114-
if (.not.copy_a) deallocate(amat)
125+
if (copy_a) deallocate(amat)
115126

116127
! Process output and return
117-
1 call linalg_error_handling(err0,err)
128+
call linalg_error_handling(err0,err)
118129

119130
end function stdlib_linalg_${ri}$solve${ndsuf}$
120131

132+
! Compute the solution to a real system of linear equations A * X = B (pure interface)
133+
pure function stdlib_linalg_${ri}$_pure_solve${ndsuf}$(a,b) result(x)
134+
!> Input matrix a[n,n]
135+
${rt}$, intent(in), target :: a(:,:)
136+
!> Right hand side vector or array, b[n] or b[n,nrhs]
137+
${rt}$, intent(in) :: b${nd}$
138+
!> Result array/matrix x[n] or x[n,nrhs]
139+
${rt}$, allocatable, target :: x${nd}$
140+
141+
!> Local variables
142+
type(linalg_state_type) :: err0
143+
integer(ilp) :: lda,n,ldb,nrhs,info
144+
integer(ilp), allocatable :: ipiv(:)
145+
${rt}$, pointer :: xmat(:,:)
146+
${rt}$, allocatable :: amat(:,:)
147+
character(*), parameter :: this = 'solve'
148+
149+
!> Problem sizes
150+
lda = size(a,1,kind=ilp)
151+
n = size(a,2,kind=ilp)
152+
ldb = size(b,1,kind=ilp)
153+
nrhs = size(b ,kind=ilp)/ldb
154+
155+
if (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then
156+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], &
157+
', b=',[ldb,nrhs])
158+
allocate(x${nde}$)
159+
call linalg_error_handling(err0)
160+
return
161+
end if
162+
163+
! Pivot indices
164+
allocate(ipiv(n))
165+
166+
! Initialize a matrix temporary
167+
allocate(amat,source=a)
168+
169+
! Initialize solution with the rhs
170+
allocate(x,source=b)
171+
xmat(1:n,1:nrhs) => x
172+
173+
! Solve system
174+
call gesv(n,nrhs,amat,lda,ipiv,xmat,ldb,info)
175+
176+
! Process output
177+
call handle_gesv_info(info,lda,n,nrhs,err0)
178+
179+
deallocate(amat)
180+
181+
! Process output and return
182+
call linalg_error_handling(err0)
183+
184+
end function stdlib_linalg_${ri}$_pure_solve${ndsuf}$
121185
#:endif
122186
#:endfor
123187
#:endfor

0 commit comments

Comments
 (0)