50 real(kind=
rp),
allocatable :: p(:)
51 real(kind=
rp),
allocatable :: q(:)
52 real(kind=
rp),
allocatable :: r(:)
53 real(kind=
rp),
allocatable :: s(:)
54 real(kind=
rp),
allocatable :: u(:)
55 real(kind=
rp),
allocatable :: w(:)
56 real(kind=
rp),
allocatable :: z(:)
57 real(kind=
rp),
allocatable :: mi(:)
58 real(kind=
rp),
allocatable :: ni(:)
70 class(
pc_t),
optional,
intent(inout),
target :: M
71 integer,
intent(in) :: n
72 integer,
intent(in) :: max_iter
73 real(kind=
rp),
optional,
intent(inout) :: rel_tol
74 real(kind=
rp),
optional,
intent(inout) :: abs_tol
91 if (
present(rel_tol) .and.
present(abs_tol))
then
92 call this%ksp_init(max_iter, rel_tol, abs_tol)
93 else if (
present(rel_tol))
then
94 call this%ksp_init(max_iter, rel_tol=rel_tol)
95 else if (
present(abs_tol))
then
96 call this%ksp_init(max_iter, abs_tol=abs_tol)
98 call this%ksp_init(max_iter)
109 if (
allocated(this%p))
then
112 if (
allocated(this%q))
then
115 if (
allocated(this%r))
then
118 if (
allocated(this%s))
then
121 if (
allocated(this%u))
then
124 if (
allocated(this%w))
then
127 if (
allocated(this%z))
then
130 if (
allocated(this%mi))
then
133 if (
allocated(this%ni))
then
143 function sx_pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
145 class(
ax_t),
intent(inout) :: ax
146 type(
field_t),
intent(inout) :: x
147 integer,
intent(in) :: n
148 real(kind=
rp),
dimension(n),
intent(inout) :: f
149 type(
coef_t),
intent(inout) :: coef
151 type(
gs_t),
intent(inout) :: gs_h
153 integer,
optional,
intent(in) :: niter
154 integer :: iter, max_iter, i, ierr
155 real(kind=
rp) :: rnorm, rtr, reduction(3), norm_fac
156 real(kind=
rp) :: alpha, beta, gamma1, gamma2, delta
157 real(kind=
rp) :: tmp1, tmp2, tmp3
158 type(mpi_request) :: request
159 type(mpi_status) :: status
161 if (
present(niter))
then
164 max_iter = this%max_iter
166 norm_fac = 1.0_rp / sqrt(coef%volume)
169 x%x(i,1,1,1) = 0.0_rp
177 call this%M%solve(this%u, this%r, n)
178 call ax%compute(this%w, this%u, coef, x%msh, x%Xh)
179 call gs_h%op(this%w, n, gs_op_add)
182 rtr =
glsc3(this%r, coef%mult, this%r, n)
183 rnorm = sqrt(rtr)*norm_fac
184 ksp_results%res_start = rnorm
185 ksp_results%res_final = rnorm
187 if(
abscmp(rnorm, 0.0_rp))
return
191 do iter = 1, max_iter
197 tmp1 = tmp1 + this%r(i) * coef%mult(i,1,1,1) * this%u(i)
198 tmp2 = tmp2 + this%w(i) * coef%mult(i,1,1,1) * this%u(i)
199 tmp3 = tmp3 + this%r(i) * coef%mult(i,1,1,1) * this%r(i)
205 call mpi_iallreduce(mpi_in_place, reduction, 3, &
208 call this%M%solve(this%mi, this%w, n)
209 call ax%compute(this%ni, this%mi, coef, x%msh, x%Xh)
210 call gs_h%op(this%ni, n, gs_op_add)
213 call mpi_wait(request, status, ierr)
215 gamma1 = reduction(1)
219 rnorm = sqrt(rtr)*norm_fac
220 if (rnorm .lt. this%abs_tol)
then
224 if (iter .gt. 1)
then
225 beta = gamma1 / gamma2
226 alpha = gamma1 / (delta - (beta * gamma1/alpha))
233 this%z(i) = beta * this%z(i) + this%ni(i)
234 this%q(i) = beta * this%q(i) + this%mi(i)
235 this%s(i) = beta * this%s(i) + this%w(i)
236 this%p(i) = beta * this%p(i) + this%u(i)
240 x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * this%p(i)
241 this%r(i) = this%r(i) - alpha * this%s(i)
242 this%u(i) = this%u(i) - alpha * this%q(i)
243 this%w(i) = this%w(i) - alpha * this%z(i)
248 ksp_results%res_final = rnorm
249 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 .
integer, parameter, public rp
Global precision used in computations.
Defines a pipelined Conjugate Gradient methods SX-Aurora backend.
subroutine sx_pipecg_init(this, n, max_iter, M, rel_tol, abs_tol)
Initialise a pipelined PCG solver.
subroutine sx_pipecg_free(this)
Deallocate a pipelined PCG solver.
type(ksp_monitor_t) function sx_pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Pipelined PCG solve.
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 for SX-Aurora.
Defines a canonical Krylov preconditioner.