Skip to content

Commit e5fa984

Browse files
jalveszjalvesz
committed
change signature of linop apply to matvec and add operator input argument
1 parent 34b83a0 commit e5fa984

File tree

5 files changed

+47
-32
lines changed

5 files changed

+47
-32
lines changed

‎doc/specs/stdlib_linalg_iterative_solvers.md

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -23,13 +23,13 @@ The `linop_<kind>_type` derive type is an auxiliary class enabling to abstract t
2323

2424
The following type-bound procedures pointer enable customization of the solver:
2525

26-
##### `apply`
26+
##### `matvec`
2727

28-
Proxy procedure for the matrix-vector product $y = alpha * M * x + beta * y$.
28+
Proxy procedure for the matrix-vector product $y = alpha * op(M) * x + beta * y$.
2929

3030
#### Syntax
3131

32-
`call ` [[stdlib_iterative_solvers(module):apply(interface)]] ` (x,y,alpha,beta)`
32+
`call ` [[stdlib_iterative_solvers(module):matvec(interface)]] ` (x,y,alpha,beta,op)`
3333

3434
###### Class
3535

@@ -45,6 +45,8 @@ Subroutine
4545

4646
`beta`: scalar of `real(<kind>)`. This argument is `intent(in)`.
4747

48+
`op`: `character(1)` scalar. This argument is `intent(in)`.
49+
4850
##### `inner_product`
4951

5052
Proxy procedure for the `dot_product`.

‎example/linalg/example_solve_custom.f90

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,9 @@ subroutine solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
3434
norm_sq0 = 0.d0
3535
!-------------------------
3636
! internal memory setup
37-
op%apply => my_apply
37+
op%matvec => my_matvec
3838
op%inner_product => my_dot
39-
M%apply => my_jacobi_preconditioner
39+
M%matvec => my_jacobi_preconditioner
4040
if(present(di))then
4141
di_ => di
4242
else
@@ -69,20 +69,22 @@ subroutine solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
6969
end if
7070
workspace_ => null()
7171
contains
72-
73-
subroutine my_apply(x,y,alpha,beta)
72+
73+
subroutine my_matvec(x,y,alpha,beta,op)
7474
real(dp), intent(in) :: x(:)
7575
real(dp), intent(inout) :: y(:)
7676
real(dp), intent(in) :: alpha
7777
real(dp), intent(in) :: beta
78-
call spmv( A , x, y , alpha, beta )
78+
character(1), intent(in) :: op
79+
call spmv( A , x, y , alpha, beta , op)
7980
y = merge( 0._dp, y, di_ )
8081
end subroutine
81-
subroutine my_jacobi_preconditioner(x,y,alpha,beta)
82+
subroutine my_jacobi_preconditioner(x,y,alpha,beta,op)
8283
real(dp), intent(in) :: x(:)
8384
real(dp), intent(inout) :: y(:)
8485
real(dp), intent(in) :: alpha
8586
real(dp), intent(in) :: beta
87+
character(1), intent(in) :: op
8688
y = merge( 0._dp, diagonal * x , di_ )
8789
end subroutine
8890
pure real(dp) function my_dot(x,y) result(r)

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ module stdlib_linalg_iterative_solvers
2424
!! The `linop` type is used to define the linear operator for the iterative solvers.
2525
#:for k, t, s in R_KINDS_TYPES
2626
type, public :: linop_${s}$_type
27-
procedure(vector_sub_${s}$), nopass, pointer :: apply => null()
27+
procedure(vector_sub_${s}$), nopass, pointer :: matvec => null()
2828
procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => default_dot_${s}$
2929
end type
3030
#:endfor
@@ -42,12 +42,13 @@ module stdlib_linalg_iterative_solvers
4242

4343
abstract interface
4444
#:for k, t, s in R_KINDS_TYPES
45-
subroutine vector_sub_${s}$(x,y,alpha,beta)
45+
subroutine vector_sub_${s}$(x,y,alpha,beta,op)
4646
import :: ${k}$
4747
${t}$, intent(in) :: x(:)
4848
${t}$, intent(inout) :: y(:)
4949
${t}$, intent(in) :: alpha
5050
${t}$, intent(in) :: beta
51+
character(1), intent(in) :: op
5152
end subroutine
5253
pure ${t}$ function reduction_sub_${s}$(x,y) result(r)
5354
import :: ${k}$

‎src/stdlib_linalg_iterative_solvers_cg.fypp

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ contains
3434
if(associated(workspace%callback)) call workspace%callback(x, norm_sq0, iter)
3535

3636
R = B
37-
call A%apply(X, R, alpha= -one_${s}$, beta=one_${s}$) ! R = B - A*X
37+
call A%matvec(X, R, alpha= -one_${s}$, beta=one_${s}$, op='N') ! R = B - A*X
3838
norm_sq = A%inner_product(R, R)
3939

4040
P = R
@@ -43,8 +43,8 @@ contains
4343
beta = zero_${s}$
4444
if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter)
4545
do while( norm_sq > tolsq * norm_sq0 .and. iter < maxiter)
46-
call A%apply(P,Ap, alpha= one_${s}$, beta=zero_${s}$) ! Ap = A*P
47-
46+
call A%matvec(P,Ap, alpha= one_${s}$, beta=zero_${s}$, op='N') ! Ap = A*P
47+
4848
alpha = norm_sq / A%inner_product(P, Ap)
4949

5050
X = X + alpha * P
@@ -94,7 +94,7 @@ contains
9494

9595
!-------------------------
9696
! internal memory setup
97-
op%apply => default_matvec
97+
op%matvec => matvec
9898
! op%inner_product => default_dot
9999
if(present(di))then
100100
di_ => di
@@ -126,15 +126,19 @@ contains
126126
workspace_ => null()
127127
contains
128128

129-
subroutine default_matvec(x,y,alpha,beta)
129+
subroutine matvec(x,y,alpha,beta,op)
130+
#:if matrix == "dense"
131+
use stdlib_linalg_blas, only: gemv
132+
#:endif
130133
${t}$, intent(in) :: x(:)
131134
${t}$, intent(inout) :: y(:)
132135
${t}$, intent(in) :: alpha
133136
${t}$, intent(in) :: beta
137+
character(1), intent(in) :: op
134138
#:if matrix == "dense"
135-
y = alpha * matmul(A,x) + beta * y
139+
call gemv(op,m=size(A,1),n=size(A,2),alpha=alpha,a=A,lda=size(A,1),x=x,incx=1,beta=beta,y=y,incy=1)
136140
#:else
137-
call spmv( A , x, y , alpha, beta )
141+
call spmv( A , x, y , alpha, beta , op)
138142
#:endif
139143
y = merge( zero_${s}$, y, di_ )
140144
end subroutine

‎src/stdlib_linalg_iterative_solvers_pcg.fypp

Lines changed: 20 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -43,25 +43,25 @@ contains
4343
if ( norm_sq0 > zero_${s}$ ) then
4444

4545
R = B
46-
call A%apply(X, R, alpha= -one_${s}$, beta=one_${s}$) ! R = B - A*X
46+
call A%matvec(X, R, alpha= -one_${s}$, beta=one_${s}$, op='N') ! R = B - A*X
4747

48-
call M%apply(R,P, alpha= one_${s}$, beta=zero_${s}$) ! P = M^{-1}*R
48+
call M%matvec(R,P, alpha= one_${s}$, beta=zero_${s}$, op='N') ! P = M^{-1}*R
4949

5050
tolsq = tol*tol
5151

5252
zr1 = zero_${s}$
5353
zr2 = one_${s}$
5454
do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) )
55-
56-
call M%apply(R,S, alpha= one_${s}$, beta=zero_${s}$) ! S = M^{-1}*R
55+
56+
call M%matvec(R,S, alpha= one_${s}$, beta=zero_${s}$, op='N') ! S = M^{-1}*R
5757
zr2 = A%inner_product( R, S )
5858

5959
if (iter>0) then
6060
beta = zr2 / zr1
6161
P = S + beta * P
6262
end if
63-
64-
call A%apply(P, Q, alpha= one_${s}$, beta=zero_${s}$) ! Q = A*P
63+
64+
call A%matvec(P, Q, alpha= one_${s}$, beta=zero_${s}$, op='N') ! Q = A*P
6565
zv2 = A%inner_product( P, Q )
6666

6767
alpha = zr2 / zv2
@@ -134,14 +134,14 @@ contains
134134
#:else
135135
call diag(A,diagonal)
136136
#:endif
137-
M_%apply => precond_jacobi
137+
M_%matvec => precond_jacobi
138138
case default
139-
M_%apply => precond_none
139+
M_%matvec => precond_none
140140
end select
141141
where(abs(diagonal)>epsilon(zero_${s}$)) diagonal = one_${s}$/diagonal
142142
end if
143143
! matvec for the operator
144-
op%apply => matvec
144+
op%matvec => matvec
145145

146146
! direchlet boundary conditions mask
147147
if(present(di))then
@@ -176,31 +176,37 @@ contains
176176
workspace_ => null()
177177
contains
178178

179-
subroutine matvec(x,y,alpha,beta)
179+
subroutine matvec(x,y,alpha,beta,op)
180+
#:if matrix == "dense"
181+
use stdlib_linalg_blas, only: gemv
182+
#:endif
180183
${t}$, intent(in) :: x(:)
181184
${t}$, intent(inout) :: y(:)
182185
${t}$, intent(in) :: alpha
183186
${t}$, intent(in) :: beta
187+
character(1), intent(in) :: op
184188
#:if matrix == "dense"
185-
y = alpha * matmul(A,x) + beta * y
189+
call gemv(op,m=size(A,1),n=size(A,2),alpha=alpha,a=A,lda=size(A,1),x=x,incx=1,beta=beta,y=y,incy=1)
186190
#:else
187-
call spmv( A , x, y , alpha, beta )
191+
call spmv( A , x, y , alpha, beta , op)
188192
#:endif
189193
y = merge( zero_${s}$, y, di_ )
190194
end subroutine
191195

192-
subroutine precond_none(x,y,alpha,beta)
196+
subroutine precond_none(x,y,alpha,beta,op)
193197
${t}$, intent(in) :: x(:)
194198
${t}$, intent(inout) :: y(:)
195199
${t}$, intent(in) :: alpha
196200
${t}$, intent(in) :: beta
201+
character(1), intent(in) :: op
197202
y = merge( zero_${s}$, x, di_ )
198203
end subroutine
199-
subroutine precond_jacobi(x,y,alpha,beta)
204+
subroutine precond_jacobi(x,y,alpha,beta,op)
200205
${t}$, intent(in) :: x(:)
201206
${t}$, intent(inout) :: y(:)
202207
${t}$, intent(in) :: alpha
203208
${t}$, intent(in) :: beta
209+
character(1), intent(in) :: op
204210
y = merge( zero_${s}$, diagonal * x, di_ ) ! inverted diagonal
205211
end subroutine
206212
end subroutine

0 commit comments

Comments
 (0)