47 use,
intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
52 real(kind=
rp),
allocatable :: w(:)
53 real(kind=
rp),
allocatable :: r(:)
54 real(kind=
rp),
allocatable :: p(:)
55 real(kind=
rp),
allocatable :: z(:)
56 type(c_ptr) :: w_d = c_null_ptr
57 type(c_ptr) :: r_d = c_null_ptr
58 type(c_ptr) :: p_d = c_null_ptr
59 type(c_ptr) :: z_d = c_null_ptr
60 type(c_ptr) :: gs_event = c_null_ptr
73 class(
pc_t),
optional,
intent(in),
target :: M
74 integer,
intent(in) :: n
75 integer,
intent(in) :: max_iter
76 real(kind=
rp),
optional,
intent(in) :: rel_tol
77 real(kind=
rp),
optional,
intent(in) :: abs_tol
78 logical,
optional,
intent(in) :: monitor
96 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
97 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
98 else if (
present(rel_tol) .and.
present(abs_tol))
then
99 call this%ksp_init(max_iter, rel_tol, abs_tol)
100 else if (
present(monitor) .and.
present(abs_tol))
then
101 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
102 else if (
present(rel_tol) .and.
present(monitor))
then
103 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
104 else if (
present(rel_tol))
then
105 call this%ksp_init(max_iter, rel_tol = rel_tol)
106 else if (
present(abs_tol))
then
107 call this%ksp_init(max_iter, abs_tol = abs_tol)
108 else if (
present(monitor))
then
109 call this%ksp_init(max_iter, monitor = monitor)
111 call this%ksp_init(max_iter)
123 if (
allocated(this%w))
then
127 if (
allocated(this%r))
then
131 if (
allocated(this%p))
then
135 if (
allocated(this%z))
then
141 if (c_associated(this%w_d))
then
145 if (c_associated(this%r_d))
then
149 if (c_associated(this%p_d))
then
153 if (c_associated(this%z_d))
then
157 if (c_associated(this%gs_event))
then
164 function cg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
166 class(
ax_t),
intent(in) :: ax
167 type(
field_t),
intent(inout) :: x
168 integer,
intent(in) :: n
169 real(kind=
rp),
dimension(n),
intent(in) :: f
170 type(
coef_t),
intent(inout) :: coef
172 type(
gs_t),
intent(inout) :: gs_h
174 integer,
optional,
intent(in) :: niter
175 real(kind=
rp),
parameter :: one = 1.0
176 real(kind=
rp),
parameter :: zero = 0.0
177 integer :: iter, max_iter
178 real(kind=
rp) :: rnorm, rtr, rtr0, rtz2, rtz1
179 real(kind=
rp) :: beta, pap, alpha, alphm, norm_fac
184 if (
present(niter))
then
187 max_iter = this%max_iter
189 norm_fac = one/sqrt(coef%volume)
197 rnorm = sqrt(rtr)*norm_fac
198 ksp_results%res_start = rnorm
199 ksp_results%res_final = rnorm
201 if(
abscmp(rnorm, zero))
return
202 call this%monitor_start(
'CG')
203 do iter = 1, max_iter
204 call this%M%solve(this%z, this%r, n)
208 if (iter .eq. 1) beta = zero
211 call ax%compute(this%w, this%p, coef, x%msh, x%Xh)
212 call gs_h%op(this%w, n, gs_op_add, this%gs_event)
214 call blst%apply(this%w, n)
224 if (iter .eq. 1) rtr0 = rtr
225 rnorm = sqrt(rtr)*norm_fac
226 call this%monitor_iter(iter, rnorm)
227 if (rnorm .lt. this%abs_tol)
then
231 call this%monitor_stop()
232 ksp_results%res_final = rnorm
233 ksp_results%iter = iter
234 ksp_results%converged = this%is_converged(iter, rnorm)
240 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
242 class(
ax_t),
intent(in) :: ax
243 type(
field_t),
intent(inout) :: x
244 type(
field_t),
intent(inout) :: y
245 type(
field_t),
intent(inout) :: z
246 integer,
intent(in) :: n
247 real(kind=
rp),
dimension(n),
intent(in) :: fx
248 real(kind=
rp),
dimension(n),
intent(in) :: fy
249 real(kind=
rp),
dimension(n),
intent(in) :: fz
250 type(
coef_t),
intent(inout) :: coef
254 type(
gs_t),
intent(inout) :: gs_h
256 integer,
optional,
intent(in) :: niter
258 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
259 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
260 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
Return the device pointer for an associated Fortran array.
Map a Fortran array to a device (allocate and associate)
Defines a Matrix-vector product.
Defines various Conjugate Gradient methods for accelerators.
type(ksp_monitor_t) function cg_device_solve(this, ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
subroutine cg_device_init(this, n, max_iter, m, rel_tol, abs_tol, monitor)
Initialise a device based PCG solver.
subroutine cg_device_free(this)
Deallocate a device based PCG solver.
type(ksp_monitor_t) function, dimension(3) cg_device_solve_coupled(this, ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard PCG coupled solve.
subroutine, public device_add2s1(a_d, b_d, c1, n)
subroutine, public device_rzero(a_d, n)
Zero a real vector.
subroutine, public device_add2s2(a_d, b_d, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
real(kind=rp) function, public device_glsc3(a_d, b_d, c_d, n)
Weighted inner product .
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
Device abstraction, common interface for various accelerators.
subroutine, public device_event_sync(event)
Synchronize an event.
subroutine, public device_free(x_d)
Deallocate memory on the device.
subroutine, public device_event_destroy(event)
Destroy a device event.
subroutine, public device_event_create(event, flags)
Create a device event queue.
Implements the base abstract type for Krylov solvers plus helper types.
integer, parameter, public ksp_max_iter
Maximum number of iters.
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of allocatable `bc_t`. Follows the standard interface of lists.
Device based preconditioned conjugate gradient method.
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 .
Defines a canonical Krylov preconditioner.