Skip to content

Commit 9106676

Browse files
jalveszjalvesz
committed
add pccg solver and example
1 parent 9ed419f commit 9106676

6 files changed

+244
-3
lines changed

‎example/linalg/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@ ADD_EXAMPLE(get_norm)
4141
ADD_EXAMPLE(solve1)
4242
ADD_EXAMPLE(solve2)
4343
ADD_EXAMPLE(solve3)
44+
ADD_EXAMPLE(solve_pccg)
4445
ADD_EXAMPLE(sparse_from_ijv)
4546
ADD_EXAMPLE(sparse_data_accessors)
4647
ADD_EXAMPLE(sparse_spmv)

‎example/linalg/example_solve_pccg.f90

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
program example_solve_pccg
2+
use stdlib_kinds, only: dp
3+
use stdlib_sparse
4+
use stdlib_linalg_iterative_solvers, only: solve_pccg
5+
6+
type(CSR_dp_type) :: laplacian_csr
7+
type(COO_dp_type) :: COO
8+
real(dp) :: laplacian(5,5)
9+
real(dp) :: x(5), load(5)
10+
logical(1) :: dirichlet(5)
11+
12+
laplacian = reshape( [1, -1, 0, 0, 0,&
13+
-1, 2, -1, 0, 0,&
14+
0, -1, 2, -1, 0,&
15+
0, 0, -1, 2, -1,&
16+
0, 0, 0, -1, 1] , [5,5])
17+
call dense2coo(laplacian,COO)
18+
call coo2csr(COO,laplacian_csr)
19+
20+
x = 0._dp
21+
load = dble( [0,0,5,0,0] )
22+
23+
dirichlet = .false._1
24+
dirichlet([1,5]) = .true._1
25+
26+
call solve_pccg(laplacian, load, x, tol=1.d-6, di=dirichlet)
27+
print *, x
28+
x = 0._dp
29+
30+
call solve_pccg(laplacian_csr, load, x, tol=1.d-6, di=dirichlet)
31+
print *, x
32+
end program

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,7 @@ set(cppFiles
9494
stdlib_linalg_lapack.fypp
9595
stdlib_linalg_iterative_solvers.fypp
9696
stdlib_linalg_iterative_solvers_cg.fypp
97+
stdlib_linalg_iterative_solvers_pccg.fypp
9798
)
9899

99100
add_subdirectory(blas)

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 42 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,8 +44,9 @@ module stdlib_linalg_iterative_solvers
4444
#:for k, t, s in R_KINDS_TYPES
4545
module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace)
4646
class(linop_${s}$), intent(in) :: A
47-
${t}$, intent(in) :: b(:), tol
47+
${t}$, intent(in) :: b(:)
4848
${t}$, intent(inout) :: x(:)
49+
${t}$, intent(in) :: tol
4950
integer, intent(in) :: maxiter
5051
type(solver_workspace_${s}$), intent(inout) :: workspace
5152
end subroutine
@@ -62,14 +63,53 @@ module stdlib_linalg_iterative_solvers
6263
#:else
6364
type(${matrix}$_${s}$_type), intent(in) :: A
6465
#:endif
65-
${t}$, intent(in) :: b(:), tol
66+
${t}$, intent(in) :: b(:)
6667
${t}$, intent(inout) :: x(:)
68+
${t}$, intent(in), optional :: tol
6769
integer, intent(in), optional :: maxiter
6870
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
6971
end subroutine
7072
#:endfor
7173
#:endfor
7274
end interface
7375
public :: solve_cg
76+
77+
interface solve_pccg_generic
78+
#:for k, t, s in R_KINDS_TYPES
79+
module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
80+
class(linop_${s}$), intent(in) :: A
81+
${t}$, intent(in) :: b(:)
82+
${t}$, intent(inout) :: x(:)
83+
${t}$, intent(in) :: tol
84+
logical(1), intent(in) :: di(:)
85+
integer, intent(in) :: maxiter
86+
logical, intent(in) :: restart
87+
type(solver_workspace_${s}$), intent(inout) :: workspace
88+
end subroutine
89+
#:endfor
90+
end interface
91+
public :: solve_pccg_generic
92+
93+
interface solve_pccg
94+
#:for matrix in MATRIX_TYPES
95+
#:for k, t, s in R_KINDS_TYPES
96+
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
97+
#:if matrix == "dense"
98+
${t}$, intent(in) :: A(:,:)
99+
#:else
100+
type(${matrix}$_${s}$_type), intent(in) :: A
101+
#:endif
102+
${t}$, intent(in) :: b(:)
103+
${t}$, intent(inout) :: x(:)
104+
${t}$, intent(in), optional :: tol
105+
logical(1), intent(in), optional, target :: di(:)
106+
integer, intent(in), optional :: maxiter
107+
logical, intent(in), optional :: restart
108+
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
109+
end subroutine
110+
#:endfor
111+
#:endfor
112+
end interface
113+
public :: solve_pccg
74114

75115
end module stdlib_linalg_iterative_solvers

‎src/stdlib_linalg_iterative_solvers_cg.fypp

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,22 +57,28 @@ contains
5757
#:else
5858
type(${matrix}$_${s}$_type), intent(in) :: A
5959
#:endif
60-
${t}$, intent(in) :: b(:), tol
60+
${t}$, intent(in) :: b(:)
6161
${t}$, intent(inout) :: x(:)
62+
${t}$, intent(in), optional :: tol
6263
integer, intent(in), optional :: maxiter
6364
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
6465
!-------------------------
6566
type(linop_${s}$) :: op
6667
type(solver_workspace_${s}$), pointer :: workspace_
6768
integer :: n, maxiter_
69+
${t}$ :: tol_
6870
!-------------------------
6971
n = size(b)
7072
op%matvec => default_matvec
7173
op%inner_product => default_dot
7274

7375
maxiter_ = n
7476
if(present(maxiter)) maxiter_ = maxiter
77+
tol_ = 1.e-4_${s}$
78+
if(present(tol)) tol_ = tol
79+
7580
if(present(workspace)) then
81+
if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3) )
7682
workspace_ => workspace
7783
else
7884
allocate( workspace_%tmp(n,3) )

‎src/stdlib_linalg_iterative_solvers_pccg.fypp

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
5+
#:set RANKS = range(1, 2+1)
6+
7+
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg
8+
use stdlib_kinds
9+
use stdlib_sparse
10+
use stdlib_constants
11+
use stdlib_intrinsics, only: dot_product => stdlib_dot_product
12+
use stdlib_linalg_iterative_solvers
13+
implicit none
14+
15+
contains
16+
17+
#:for k, t, s in R_KINDS_TYPES
18+
module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
19+
class(linop_${s}$), intent(in) :: A
20+
${t}$, intent(in) :: b(:), tol
21+
${t}$, intent(inout) :: x(:)
22+
logical(1), intent(in) :: di(:)
23+
integer, intent(in) :: maxiter
24+
logical, intent(in) :: restart
25+
type(solver_workspace_${s}$), intent(inout) :: workspace
26+
!-------------------------
27+
integer :: iter
28+
${t}$ :: norm_sq, norm_sq0, norm_sq_old, residual0
29+
${t}$ :: zr1, zr2, zv2, alpha, beta, tolsq
30+
!-------------------------
31+
associate( R => workspace%tmp(:,1), &
32+
S => workspace%tmp(:,2), &
33+
P => workspace%tmp(:,3), &
34+
Q => workspace%tmp(:,4))
35+
norm_sq = A%inner_product( b, b )
36+
norm_sq0 = norm_sq
37+
residual0 = sqrt(norm_sq0)
38+
if ( norm_sq0 > zero_${s}$ ) then
39+
if(restart) X = zero_${s}$
40+
where( di ) X = B
41+
42+
call A%matvec(X, R)
43+
R = B - R
44+
where( di ) R = zero_${s}$
45+
46+
call A%preconditionner(R,P)
47+
where( di ) P = zero_${s}$
48+
49+
tolsq = tol*tol
50+
iter = 0
51+
zr1 = zero_${s}$
52+
zr2 = one_${s}$
53+
do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) )
54+
55+
call A%preconditionner(R,S)
56+
where ( di ) S = zero_${s}$
57+
zr2 = A%inner_product( R, S )
58+
59+
if (iter>0) then
60+
beta = zr2 / zr1
61+
P = S + beta * P
62+
end if
63+
64+
call A%matvec(P, Q)
65+
where( di ) Q = zero_${s}$
66+
zv2 = A%inner_product( P, Q )
67+
68+
alpha = zr2 / zv2
69+
70+
X = X + alpha * P
71+
R = R - alpha * Q
72+
norm_sq = A%inner_product( R, R )
73+
norm_sq_old = norm_sq
74+
zr1 = zr2
75+
iter = iter + 1
76+
end do
77+
end if
78+
end associate
79+
end subroutine
80+
#:endfor
81+
82+
#:for matrix in MATRIX_TYPES
83+
#:for k, t, s in R_KINDS_TYPES
84+
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
85+
#:if matrix == "dense"
86+
${t}$, intent(in) :: A(:,:)
87+
#:else
88+
type(${matrix}$_${s}$_type), intent(in) :: A
89+
#:endif
90+
${t}$, intent(in) :: b(:)
91+
${t}$, intent(inout) :: x(:)
92+
${t}$, intent(in), optional :: tol
93+
logical(1), intent(in), optional, target :: di(:)
94+
integer, intent(in), optional :: maxiter
95+
logical, intent(in), optional :: restart
96+
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
97+
!-------------------------
98+
type(linop_${s}$) :: op
99+
type(solver_workspace_${s}$), pointer :: workspace_
100+
integer :: n, maxiter_
101+
${t}$ :: tol_
102+
logical :: restart_
103+
logical(1), pointer :: di_(:)
104+
!-------------------------
105+
n = size(b)
106+
107+
op%matvec => default_matvec
108+
op%inner_product => default_dot
109+
op%preconditionner => default_preconditionner
110+
111+
maxiter_ = n
112+
if(present(maxiter)) maxiter_ = maxiter
113+
restart_ = .true.
114+
if(present(restart)) restart_ = restart
115+
tol_ = 1.e-4_${s}$
116+
if(present(tol)) tol_ = tol
117+
118+
if(present(di))then
119+
di_ => di
120+
else
121+
allocate(di_(n),source=.false._1)
122+
end if
123+
124+
if(present(workspace)) then
125+
if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) )
126+
workspace_ => workspace
127+
else
128+
allocate( workspace_ )
129+
allocate( workspace_%tmp(n,4) )
130+
end if
131+
132+
call solve_pccg_generic(op,b,x,di_,tol,maxiter_,restart_,workspace_)
133+
134+
if(.not.present(di)) deallocate(di_)
135+
contains
136+
137+
subroutine default_matvec(x,y)
138+
${t}$, intent(in) :: x(:)
139+
${t}$, intent(inout) :: y(:)
140+
#:if matrix == "dense"
141+
y = matmul(A,x)
142+
#:else
143+
call spmv( A , x, y )
144+
#:endif
145+
end subroutine
146+
subroutine default_preconditionner(x,y)
147+
${t}$, intent(in) :: x(:)
148+
${t}$, intent(inout) :: y(:)
149+
y = x
150+
end subroutine
151+
pure ${t}$ function default_dot(x,y) result(r)
152+
${t}$, intent(in) :: x(:)
153+
${t}$, intent(in) :: y(:)
154+
r = dot_product(x,y)
155+
end function
156+
end subroutine
157+
158+
#:endfor
159+
#:endfor
160+
161+
end submodule stdlib_linalg_iterative_pccg

0 commit comments

Comments
 (0)