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 real(kind=
rp),
allocatable :: alpha(:)
66 subroutine cg_init(this, n, max_iter, M, rel_tol, abs_tol)
67 class(
cg_t),
intent(inout),
target :: this
68 integer,
intent(in) :: max_iter
69 class(
pc_t),
optional,
intent(inout),
target :: M
70 integer,
intent(in) :: n
71 real(kind=
rp),
optional,
intent(inout) :: rel_tol
72 real(kind=
rp),
optional,
intent(inout) :: abs_tol
86 if (
present(rel_tol) .and.
present(abs_tol))
then
87 call this%ksp_init(max_iter, rel_tol, abs_tol)
88 else if (
present(rel_tol))
then
89 call this%ksp_init(max_iter, rel_tol=rel_tol)
90 else if (
present(abs_tol))
then
91 call this%ksp_init(max_iter, abs_tol=abs_tol)
93 call this%ksp_init(max_iter)
100 class(
cg_t),
intent(inout) :: this
104 if (
allocated(this%w))
then
108 if (
allocated(this%r))
then
112 if (
allocated(this%p))
then
116 if (
allocated(this%z))
then
120 if (
allocated(this%alpha))
then
121 deallocate(this%alpha)
129 function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
130 class(
cg_t),
intent(inout) :: this
131 class(
ax_t),
intent(inout) :: ax
132 type(
field_t),
intent(inout) :: x
133 integer,
intent(in) :: n
134 real(kind=
rp),
dimension(n),
intent(inout) :: f
135 type(
coef_t),
intent(inout) :: coef
137 type(
gs_t),
intent(inout) :: gs_h
139 integer,
optional,
intent(in) :: niter
140 integer :: iter, max_iter, i, j, k, p_cur, p_prev
141 real(kind=
rp) :: rnorm, rtr, rtz2, rtz1, x_plus(neko_blk_size)
142 real(kind=
rp) :: beta, pap, norm_fac
144 if (
present(niter))
then
147 max_iter = this%max_iter
149 norm_fac = 1.0_rp / sqrt(coef%volume)
151 associate(w => this%w, r => this%r, p => this%p, &
152 z => this%z, alpha => this%alpha)
159 rtr =
glsc3(r, coef%mult, r, n)
160 rnorm = sqrt(rtr) * norm_fac
161 ksp_results%res_start = rnorm
162 ksp_results%res_final = rnorm
166 if(
abscmp(rnorm, 0.0_rp))
return
167 do iter = 1, max_iter
168 call this%M%solve(z, r, n)
170 rtz1 =
glsc3(r, coef%mult, z, n)
173 if (iter .eq. 1) beta = 0.0_rp
175 p(i,p_cur) = z(i) + beta * p(i,p_prev)
178 call ax%compute(w, p(1,p_cur), coef, x%msh, x%Xh)
179 call gs_h%op(w, n, gs_op_add)
182 pap =
glsc3(w, coef%mult, p(1,p_cur), n)
184 alpha(p_cur) = rtz1 / pap
186 rnorm = sqrt(rtr) * norm_fac
189 (rnorm .lt. this%abs_tol) .or. iter .eq. max_iter)
then
190 do i = 0, n, neko_blk_size
191 if (i + neko_blk_size .le. n)
then
192 do k = 1, neko_blk_size
196 do k = 1, neko_blk_size
197 x_plus(k) = x_plus(k) + alpha(j) * p(i+k,j)
200 do k = 1, neko_blk_size
201 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
207 x_plus(1) = x_plus(1) + alpha(j) * p(i+k,j)
209 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
215 if (rnorm .lt. this%abs_tol)
exit
223 ksp_results%res_final = rnorm
224 ksp_results%iter = iter
229 integer,
intent(in) :: n
230 real(kind=rp),
intent(inout) :: r(n), rtr
231 real(kind=rp),
intent(in) ::mult(n), w(n), alpha
236 r(i) = r(i) - alpha*w(i)
237 rtr = rtr + r(i) * r(i) * mult(i)
239 call mpi_allreduce(mpi_in_place, rtr, 1, &
240 mpi_real_precision, mpi_sum, neko_comm, ierr)
Defines a Matrix-vector product.
Defines a boundary condition.
Defines various Conjugate Gradient methods.
subroutine cg_free(this)
Deallocate a standard PCG solver.
integer, parameter cg_p_space
subroutine cg_init(this, n, max_iter, M, rel_tol, abs_tol)
Initialise a standard PCG solver.
type(ksp_monitor_t) function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
subroutine second_cg_part(rtr, r, mult, w, alpha, n)
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.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Standard 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.