Choral
cg_mod.f90
Go to the documentation of this file.
1 
13 
14 module cg_mod
15 
16  use real_type
17  use basic_tools
18  use abstract_interfaces, only: rntorn
19  use r1d_mod, only: norm_2, axpy, scalprod, xpay, copy
20 
21  implicit none
22  private
23 
24 
25 ! %----------------------------------------%
26 ! | |
27 ! | PUBLIC DATA |
28 ! | |
29 ! %----------------------------------------%
30 
31  !! GENERIC
32  public :: cg !! TESTED
33  public :: pcg !! TESTED
34 
35 contains
36 
37 
59  subroutine cg(x, iter, res, b, A, tol, itmax, verb)
60 
61  real(RP), dimension(:), intent(inout) :: x
62  integer , intent(out) :: iter
63  real(RP) , intent(out) :: res
64  real(RP), dimension(:), intent(in) :: b
65  procedure(rntorn) :: A
66  real(RP) , intent(in) :: tol
67  integer , intent(in) :: itmax, verb
68 
69 
70  real(RP), dimension(size(x,1)) :: r, p, Ap
71  real(RP) :: alpha, beta, nb2
72  real(RP) :: ps_rr, ps_App, ps_rjrj
73  integer :: nn
74 
75  if (verb>1) write(*,*)'cg_mod : cg'
76 
77  iter = 0
78 
79  nb2 = norm_2(b)
80  nn = size(x,1)
81  if (nb2/re(nn)<real_tol) nb2=sqrt(re(nn))
82  nb2 = 1._rp / nb2
83 
84  call a(r, x)
85  ! r = b-r
86  call axpy(-1._rp, r, b)
87 
88  res = norm_2(r) * nb2
89 
90  if (verb>2) write(*,*)' iter', iter,', residual = ', res
91  if (res<tol) return
92 
93  ! p = r
94  p = r
95  ps_rr = scalprod(r, r)
96  do iter=1, itmax
97  ps_rjrj = ps_rr
98 
99  call a(ap, p)
100 
101  ps_app = scalprod(ap, p)
102 
103  alpha = ps_rr / ps_app
104  ! x = x + alpha * p
105  call xpay(x, alpha, p)
106 
107  ! r = r - alpha * Ap
108  call xpay(r, -alpha, ap)
109 
110  res = norm_2(r) * nb2
111 
112  if (verb>2) write(*,*)' iter', iter,', residual = ', res
113  if (res<tol) return
114 
115  ps_rr = scalprod(r, r)
116  beta = ps_rr/ps_rjrj
117  ! p = r + beta * p
118  call axpy(beta, p, r)
119 
120  end do
121 
122  ! Convergence failure flag
123  iter = -1
124  call warning("cg_mod: cg: not converged", verb-1)
125 
126  end subroutine cg
127 
128 
154  subroutine pcg(x, iter, res, b, A, pc, tol, itmax, verb)
156  real(RP), dimension(:), intent(inout) :: x
157  integer , intent(out) :: iter
158  real(RP) , intent(out) :: res
159  real(RP), dimension(:), intent(in) :: b
160  procedure(rntorn) :: A, pc
161  real(RP) , intent(in) :: tol
162  integer , intent(in) :: itmax, verb
163 
164  real(RP), dimension(size(x,1)) :: r, p, Ap, z
165 
166  real(RP) :: alpha, beta, nb2
167  real(RP) :: ps_rz, ps_App, ps_rjzj
168  integer :: nn
169 
170  if (verb>1) write(*,*)'cg_mod : pcg'
171 
172  iter = 0
173 
174  nb2 = norm_2(b)
175  nn = (size(x,1))
176  if (nb2/re(nn)<real_tol) nb2=sqrt(re(nn))
177  nb2 = 1._rp / nb2
178 
179  call a(r, x)
180  ! r = b-r
181  call axpy(-1._rp, r, b)
182  res = norm_2(r) * nb2
183 
184  if (verb>2) write(*,*)' iter', iter,', residual = ', res
185 
186  if (res<tol) return
187 
188  call pc(z, r)
189  ! p = z
190  call copy(p, z)
191  ps_rz = scalprod(r,z)
192  do iter=1, itmax
193  ps_rjzj = ps_rz
194 
195  call a(ap, p)
196 
197  ps_app = scalprod(ap,p)
198 
199  alpha = ps_rz / ps_app
200  ! x = x + alpha * p
201  call xpay(x, alpha, p)
202  ! r = r - alpha * Ap
203  call xpay(r, -alpha, ap)
204 
205  res = norm_2(r) * nb2
206 
207  if (verb>2) write(*,*)' iter', iter, ' residual = ', res
208 
209  if (res<tol) return
210 
211  call pc(z, r)
212  ps_rz = scalprod(r,z)
213  beta = ps_rz/ps_rjzj
214  !p = z + beta * p
215  call axpy(beta, p, z)
216 
217  end do
218 
219  ! Convergence failure flag
220  iter = -1
221  call warning("cg_mod: pcg: not converged", verb-1)
222 
223  end subroutine pcg
224 
225 end module cg_mod
226 
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
subroutine, public pcg(x, iter, res, b, A, pc, tol, itmax, verb)
Preconditioned conjugate gradient.
Definition: cg_mod.f90:155
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine, public cg(x, iter, res, b, A, tol, itmax, verb)
Conjugate gradient (no preconditioning)
Definition: cg_mod.f90:60
L2 vector norm.
Definition: R1d_mod.F90:82
conversion integers or rational to real
Definition: real_type.F90:153
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
CG LINEAR SOLVER = Conjugate Gradient
Definition: cg_mod.f90:14
y = x
Definition: R1d_mod.F90:44
x = a*x + y
Definition: R1d_mod.F90:55
scalar product
Definition: R1d_mod.F90:77
x = x + b*y // OR // z = x + b*y
Definition: R1d_mod.F90:60
real(rp), parameter, public real_tol
Definition: real_type.F90:91
subroutine, public warning(message, verb)
Warning message.