Skip to content

Commit 367987a

Browse files
jalveszjalvesz
committed
add preconditionners
1 parent 2c2196a commit 367987a

File tree

2 files changed

+67
-12
lines changed

2 files changed

+67
-12
lines changed

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,7 @@ module stdlib_linalg_iterative_solvers
111111
interface solve_pccg
112112
#:for matrix in MATRIX_TYPES
113113
#:for k, t, s in R_KINDS_TYPES
114-
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace)
114+
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
115115
#:if matrix == "dense"
116116
${t}$, intent(in) :: A(:,:)
117117
#:else
@@ -123,6 +123,7 @@ module stdlib_linalg_iterative_solvers
123123
logical(1), intent(in), optional, target :: di(:)
124124
integer, intent(in), optional :: maxiter
125125
logical, intent(in), optional :: restart
126+
integer, intent(in), optional :: precond !> preconditionner method flag
126127
class(linop_${s}$), optional , intent(in), target :: M !> preconditionner
127128
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
128129
end subroutine

‎src/stdlib_linalg_iterative_solvers_pccg.fypp

Lines changed: 65 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,18 @@
77
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg
88
use stdlib_kinds
99
use stdlib_sparse
10+
use stdlib_linalg, only: diag
1011
use stdlib_constants
1112
use stdlib_linalg_iterative_solvers
1213
implicit none
1314

15+
enum, bind(c)
16+
enumerator :: none = 0
17+
enumerator :: jacobi
18+
enumerator :: ssor
19+
enumerator :: ichol
20+
end enum
21+
1422
contains
1523

1624
#:for k, t, s in R_KINDS_TYPES
@@ -78,7 +86,7 @@ contains
7886

7987
#:for matrix in MATRIX_TYPES
8088
#:for k, t, s in R_KINDS_TYPES
81-
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace)
89+
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
8290
#:if matrix == "dense"
8391
${t}$, intent(in) :: A(:,:)
8492
#:else
@@ -90,6 +98,7 @@ contains
9098
logical(1), intent(in), optional, target :: di(:)
9199
integer, intent(in), optional :: maxiter
92100
logical, intent(in), optional :: restart
101+
integer, intent(in), optional :: precond !> preconditionner method flag
93102
class(linop_${s}$), optional , intent(in), target :: M !> preconditionner
94103
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
95104
!-------------------------
@@ -101,25 +110,53 @@ contains
101110
logical :: restart_
102111
logical(1), pointer :: di_(:)
103112
!-------------------------
113+
! working data for preconditionner
114+
integer :: precond_
115+
${t}$, allocatable :: diagonal(:)
116+
#:if matrix == "dense"
117+
${t}$, allocatable:: L(:,:) !> lower triangular
118+
#:else
119+
type(${matrix}$_${s}$_type) :: L !> lower triangular
120+
#:endif
121+
!-------------------------
104122
n = size(b)
105123
maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter
106124
restart_ = .true.; if(present(restart)) restart_ = restart
107125
tol_ = 1.e-4_${s}$; if(present(tol)) tol_ = tol
108-
126+
precond_ = none ; if(present(precond)) precond_ = precond
109127
!-------------------------
110128
! internal memory setup
111-
op%apply => default_matvec
129+
op%apply => matvec
112130
if(present(M)) then
113131
M_ => M
114132
else
115133
allocate( M_ )
116-
M_%apply => default_preconditionner
134+
select case(precond_)
135+
case(jacobi)
136+
#:if matrix == "dense"
137+
diagonal = diag(A)
138+
#:else
139+
call diag(A,diagonal)
140+
#:endif
141+
where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal
142+
M_%apply => precond_jacobi
143+
case(ssor)
144+
#:if matrix == "dense"
145+
diagonal = diag(A)
146+
#:else
147+
call diag(A,diagonal)
148+
#:endif
149+
where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal
150+
! TODO: Extract the lower triangular part of A in L
151+
M_%apply => precond_ldl
152+
case(ichol)
153+
! TODO: Add LDL^t factorization
154+
M_%apply => precond_ldl
155+
case default
156+
M_%apply => precond_none
157+
end select
117158
end if
118-
!> TODO: factorize the preconditionner
119-
!> Jacobi: M = diag(A)^{-1}
120-
!> SOR: M = (w/(2-w)) (D/w+L) D^{-1} (D/w+U) ; A = D + L + U
121-
!> ILU: M = ( LU )^{-1} ; A = LU
122-
!> ICHOL: M = (L*D*U)^{-1} ; A = L*D*U
159+
123160
if(present(di))then
124161
di_ => di
125162
else
@@ -151,7 +188,7 @@ contains
151188
workspace_ => null()
152189
contains
153190

154-
subroutine default_matvec(x,y,alpha,beta)
191+
subroutine matvec(x,y,alpha,beta)
155192
${t}$, intent(in) :: x(:)
156193
${t}$, intent(inout) :: y(:)
157194
${t}$, intent(in) :: alpha
@@ -163,13 +200,30 @@ contains
163200
#:endif
164201
y = merge( zero_${s}$, y, di_ )
165202
end subroutine
166-
subroutine default_preconditionner(x,y,alpha,beta)
203+
204+
subroutine precond_none(x,y,alpha,beta)
167205
${t}$, intent(in) :: x(:)
168206
${t}$, intent(inout) :: y(:)
169207
${t}$, intent(in) :: alpha
170208
${t}$, intent(in) :: beta
171209
y = merge( zero_${s}$, x, di_ )
172210
end subroutine
211+
subroutine precond_jacobi(x,y,alpha,beta)
212+
${t}$, intent(in) :: x(:)
213+
${t}$, intent(inout) :: y(:)
214+
${t}$, intent(in) :: alpha
215+
${t}$, intent(in) :: beta
216+
y = merge( zero_${s}$, diagonal * x, di_ ) !> inverted diagonal
217+
end subroutine
218+
subroutine precond_ldl(x,y,alpha,beta)
219+
${t}$, intent(in) :: x(:)
220+
${t}$, intent(inout) :: y(:)
221+
${t}$, intent(in) :: alpha
222+
${t}$, intent(in) :: beta
223+
call solve_forward_triangular( L , x , y )
224+
y = merge( zero_${s}$, diagonal * y, di_ ) !> inverted diagonal
225+
call solve_backward_triangular( L , y , y )
226+
end subroutine
173227
end subroutine
174228

175229
#:endfor

0 commit comments

Comments
 (0)