52 real(kind=
rp),
allocatable :: p(:)
53 real(kind=
rp),
allocatable :: q(:)
54 real(kind=
rp),
allocatable :: r(:)
55 real(kind=
rp),
allocatable :: s(:)
56 real(kind=
rp),
allocatable :: u(:,:)
57 real(kind=
rp),
allocatable :: w(:)
58 real(kind=
rp),
allocatable :: z(:)
59 real(kind=
rp),
allocatable :: mi(:)
60 real(kind=
rp),
allocatable :: ni(:)
73 subroutine pipecg_init(this, n, max_iter, M, rel_tol, abs_tol)
74 class(
pipecg_t),
intent(inout) :: this
75 integer,
intent(in) :: max_iter
76 class(
pc_t),
optional,
intent(inout),
target :: M
77 integer,
intent(in) :: n
78 real(kind=
rp),
optional,
intent(inout) :: rel_tol
79 real(kind=
rp),
optional,
intent(inout) :: abs_tol
96 if (
present(rel_tol) .and.
present(abs_tol))
then
97 call this%ksp_init(max_iter, rel_tol, abs_tol)
98 else if (
present(rel_tol))
then
99 call this%ksp_init(max_iter, rel_tol=rel_tol)
100 else if (
present(abs_tol))
then
101 call this%ksp_init(max_iter, abs_tol=abs_tol)
103 call this%ksp_init(max_iter)
110 class(
pipecg_t),
intent(inout) :: this
114 if (
allocated(this%p))
then
117 if (
allocated(this%q))
then
120 if (
allocated(this%r))
then
123 if (
allocated(this%s))
then
126 if (
allocated(this%u))
then
129 if (
allocated(this%w))
then
132 if (
allocated(this%z))
then
135 if (
allocated(this%mi))
then
138 if (
allocated(this%ni))
then
148 function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
149 class(
pipecg_t),
intent(inout) :: this
150 class(
ax_t),
intent(inout) :: ax
151 type(
field_t),
intent(inout) :: x
152 integer,
intent(in) :: n
153 real(kind=
rp),
dimension(n),
intent(inout) :: f
154 type(
coef_t),
intent(inout) :: coef
156 type(
gs_t),
intent(inout) :: gs_h
158 integer,
optional,
intent(in) :: niter
159 integer :: iter, max_iter, i, j, k, ierr, p_cur, p_prev, u_prev
160 real(kind=
rp) :: rnorm, rtr, reduction(3), norm_fac
162 real(kind=
rp) :: gamma1, gamma2, delta
163 real(kind=
rp) :: tmp1, tmp2, tmp3, x_plus(neko_blk_size)
164 type(mpi_request) :: request
165 type(mpi_status) :: status
167 if (
present(niter))
then
170 max_iter = this%max_iter
172 norm_fac = 1.0_rp / sqrt(coef%volume)
174 associate(p => this%p, q => this%q, r => this%r, s => this%s, &
175 u => this%u, w => this%w, z => this%z, mi => this%mi, ni => this%ni)
186 call this%M%solve(u(1,u_prev), r, n)
187 call ax%compute(w, u(1,u_prev), coef, x%msh, x%Xh)
188 call gs_h%op(w, n, gs_op_add)
191 rtr =
glsc3(r, coef%mult, r, n)
192 rnorm = sqrt(rtr)*norm_fac
193 ksp_results%res_start = rnorm
194 ksp_results%res_final = rnorm
196 if(
abscmp(rnorm, 0.0_rp))
return
203 tmp1 = tmp1 + r(i) * coef%mult(i,1,1,1) * u(i,u_prev)
204 tmp2 = tmp2 + w(i) * coef%mult(i,1,1,1) * u(i,u_prev)
205 tmp3 = tmp3 + r(i) * coef%mult(i,1,1,1) * r(i)
211 do iter = 1, max_iter
212 call mpi_iallreduce(mpi_in_place, reduction, 3, &
215 call this%M%solve(mi, w, n)
216 call ax%compute(ni, mi, coef, x%msh, x%Xh)
217 call gs_h%op(ni, n, gs_op_add)
220 call mpi_wait(request, status, ierr)
222 gamma1 = reduction(1)
226 rnorm = sqrt(rtr)*norm_fac
227 if (rnorm .lt. this%abs_tol)
exit
229 if (iter .gt. 1)
then
230 beta(p_cur) = gamma1 / gamma2
231 alpha(p_cur) = gamma1 / (delta - (beta(p_cur) * gamma1/alpha(p_prev)))
234 alpha(p_cur) = gamma1/delta
240 do i = 0, n, neko_blk_size
241 if (i + neko_blk_size .le. n)
then
242 do k = 1, neko_blk_size
243 z(i+k) = beta(p_cur) * z(i+k) + ni(i+k)
244 q(i+k) = beta(p_cur) * q(i+k) + mi(i+k)
245 s(i+k) = beta(p_cur) * s(i+k) + w(i+k)
246 r(i+k) = r(i+k) - alpha(p_cur) * s(i+k)
247 u(i+k,p_cur) = u(i+k,u_prev) - alpha(p_cur) * q(i+k)
248 w(i+k) = w(i+k) - alpha(p_cur) * z(i+k)
249 tmp1 = tmp1 + r(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
250 tmp2 = tmp2 + w(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
251 tmp3 = tmp3 + r(i+k) * coef%mult(i+k,1,1,1) * r(i+k)
255 z(i+k) = beta(p_cur) * z(i+k) + ni(i+k)
256 q(i+k) = beta(p_cur) * q(i+k) + mi(i+k)
257 s(i+k) = beta(p_cur) * s(i+k) + w(i+k)
258 r(i+k) = r(i+k) - alpha(p_cur) * s(i+k)
259 u(i+k,p_cur) = u(i+k,u_prev) - alpha(p_cur) * q(i+k)
260 w(i+k) = w(i+k) - alpha(p_cur) * z(i+k)
261 tmp1 = tmp1 + r(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
262 tmp2 = tmp2 + w(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
263 tmp3 = tmp3 + r(i+k) * coef%mult(i+k,1,1,1) * r(i+k)
273 do i = 0, n, neko_blk_size
274 if (i + neko_blk_size .le. n)
then
275 do k = 1, neko_blk_size
280 do k = 1, neko_blk_size
281 p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
282 x_plus(k) = x_plus(k) + alpha(j) * p(i+k)
286 do k = 1, neko_blk_size
287 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
295 p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
296 x_plus(1) = x_plus(1) + alpha(j) * p(i+k)
299 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
306 alpha(1) = alpha(p_cur)
307 beta(1) = beta(p_cur)
316 if ( p_cur .ne. 1)
then
317 do i = 0, n, neko_blk_size
318 if (i + neko_blk_size .le. n)
then
319 do k = 1, neko_blk_size
324 do k = 1, neko_blk_size
325 p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
326 x_plus(k) = x_plus(k) + alpha(j) * p(i+k)
330 do k = 1, neko_blk_size
331 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
339 p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
340 x_plus(1) = x_plus(1) + alpha(j) * p(i+k)
343 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
350 ksp_results%res_final = rnorm
351 ksp_results%iter = iter
Defines a Matrix-vector product.
Defines a boundary condition.
type(mpi_comm) neko_comm
MPI communicator.
type(mpi_datatype) mpi_real_precision
MPI type for working precision of REAL types.
Implements the base abstract type for Krylov solvers plus helper types.
integer, parameter, public ksp_max_iter
Maximum number of iters.
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer, parameter, public rp
Global precision used in computations.
Defines a pipelined Conjugate Gradient methods.
subroutine pipecg_free(this)
Deallocate a pipelined PCG solver.
integer, parameter pipecg_p_space
type(ksp_monitor_t) function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Pipelined PCG solve.
subroutine pipecg_init(this, n, max_iter, M, rel_tol, abs_tol)
Initialise a pipelined PCG solver.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Type for storing initial and final residuals in a Krylov solver.
Base abstract type for a canonical Krylov method, solving .
Pipelined preconditioned conjugate gradient method.
Defines a canonical Krylov preconditioner.