51 real(kind=
rp),
allocatable :: w(:)
52 real(kind=
rp),
allocatable :: r(:)
53 real(kind=
rp),
allocatable :: p(:)
54 real(kind=
rp),
allocatable :: z(:)
55 type(c_ptr) :: w_d = c_null_ptr
56 type(c_ptr) :: r_d = c_null_ptr
57 type(c_ptr) :: p_d = c_null_ptr
58 type(c_ptr) :: z_d = c_null_ptr
59 type(c_ptr) :: gs_event = c_null_ptr
72 class(
pc_t),
optional,
intent(inout),
target :: M
73 integer,
intent(in) :: n
74 integer,
intent(in) :: max_iter
75 real(kind=
rp),
optional,
intent(inout) :: rel_tol
76 real(kind=
rp),
optional,
intent(inout) :: abs_tol
77 logical,
optional,
intent(in) :: monitor
95 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
96 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
97 else if (
present(rel_tol) .and.
present(abs_tol))
then
98 call this%ksp_init(max_iter, rel_tol, abs_tol)
99 else if (
present(monitor) .and.
present(abs_tol))
then
100 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
101 else if (
present(rel_tol) .and.
present(monitor))
then
102 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
103 else if (
present(rel_tol))
then
104 call this%ksp_init(max_iter, rel_tol = rel_tol)
105 else if (
present(abs_tol))
then
106 call this%ksp_init(max_iter, abs_tol = abs_tol)
107 else if (
present(monitor))
then
108 call this%ksp_init(max_iter, monitor = monitor)
110 call this%ksp_init(max_iter)
122 if (
allocated(this%w))
then
126 if (
allocated(this%r))
then
130 if (
allocated(this%p))
then
134 if (
allocated(this%z))
then
140 if (c_associated(this%w_d))
then
144 if (c_associated(this%r_d))
then
148 if (c_associated(this%p_d))
then
152 if (c_associated(this%z_d))
then
156 if (c_associated(this%gs_event))
then
163 function cg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
165 class(
ax_t),
intent(inout) :: ax
166 type(
field_t),
intent(inout) :: x
167 integer,
intent(in) :: n
168 real(kind=
rp),
dimension(n),
intent(inout) :: f
169 type(
coef_t),
intent(inout) :: coef
171 type(
gs_t),
intent(inout) :: gs_h
173 integer,
optional,
intent(in) :: niter
174 real(kind=
rp),
parameter :: one = 1.0
175 real(kind=
rp),
parameter :: zero = 0.0
176 integer :: iter, max_iter
177 real(kind=
rp) :: rnorm, rtr, rtr0, rtz2, rtz1
178 real(kind=
rp) :: beta, pap, alpha, alphm, norm_fac
183 if (
present(niter))
then
186 max_iter = this%max_iter
188 norm_fac = one/sqrt(coef%volume)
196 rnorm = sqrt(rtr)*norm_fac
197 ksp_results%res_start = rnorm
198 ksp_results%res_final = rnorm
200 if(
abscmp(rnorm, zero))
return
201 call this%monitor_start(
'CG')
202 do iter = 1, max_iter
203 call this%M%solve(this%z, this%r, n)
207 if (iter .eq. 1) beta = zero
210 call ax%compute(this%w, this%p, coef, x%msh, x%Xh)
211 call gs_h%op(this%w, n, gs_op_add, this%gs_event)
223 if (iter .eq. 1) rtr0 = rtr
224 rnorm = sqrt(rtr)*norm_fac
225 call this%monitor_iter(iter, rnorm)
226 if (rnorm .lt. this%abs_tol)
then
230 call this%monitor_stop()
231 ksp_results%res_final = rnorm
232 ksp_results%iter = iter
238 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
240 class(
ax_t),
intent(inout) :: ax
241 type(
field_t),
intent(inout) :: x
242 type(
field_t),
intent(inout) :: y
243 type(
field_t),
intent(inout) :: z
244 integer,
intent(in) :: n
245 real(kind=
rp),
dimension(n),
intent(inout) :: fx
246 real(kind=
rp),
dimension(n),
intent(inout) :: fy
247 real(kind=
rp),
dimension(n),
intent(inout) :: fz
248 type(
coef_t),
intent(inout) :: coef
252 type(
gs_t),
intent(inout) :: gs_h
254 integer,
optional,
intent(in) :: niter
256 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
257 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
258 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 a boundary condition.
Defines various Conjugate Gradient methods for accelerators.
subroutine cg_device_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a device based PCG solver.
type(ksp_monitor_t) function cg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
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 boundary conditions.
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.