49 real(kind=
rp),
allocatable :: w(:)
50 real(kind=
rp),
allocatable :: r(:)
51 real(kind=
rp),
allocatable :: p(:)
52 real(kind=
rp),
allocatable :: z(:)
63 subroutine sx_cg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
64 class(
sx_cg_t),
intent(inout) :: this
65 class(
pc_t),
optional,
intent(inout),
target :: M
66 integer,
intent(in) :: n
67 integer,
intent(in) :: max_iter
68 real(kind=
rp),
optional,
intent(inout) :: rel_tol
69 real(kind=
rp),
optional,
intent(inout) :: abs_tol
70 logical,
optional,
intent(in) :: monitor
83 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
84 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
85 else if (
present(rel_tol) .and.
present(abs_tol))
then
86 call this%ksp_init(max_iter, rel_tol, abs_tol)
87 else if (
present(monitor) .and.
present(abs_tol))
then
88 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
89 else if (
present(rel_tol) .and.
present(monitor))
then
90 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
91 else if (
present(rel_tol))
then
92 call this%ksp_init(max_iter, rel_tol = rel_tol)
93 else if (
present(abs_tol))
then
94 call this%ksp_init(max_iter, abs_tol = abs_tol)
95 else if (
present(monitor))
then
96 call this%ksp_init(max_iter, monitor = monitor)
98 call this%ksp_init(max_iter)
105 class(
sx_cg_t),
intent(inout) :: this
109 if (
allocated(this%w))
then
113 if (
allocated(this%r))
then
117 if (
allocated(this%p))
then
121 if (
allocated(this%z))
then
130 function sx_cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
131 class(
sx_cg_t),
intent(inout) :: this
132 class(
ax_t),
intent(inout) :: ax
133 type(
field_t),
intent(inout) :: x
134 integer,
intent(in) :: n
135 real(kind=
rp),
dimension(n),
intent(inout) :: f
136 type(
coef_t),
intent(inout) :: coef
138 type(
gs_t),
intent(inout) :: gs_h
140 integer,
optional,
intent(in) :: niter
141 real(kind=
rp),
parameter :: one = 1.0
142 real(kind=
rp),
parameter :: zero = 0.0
143 integer :: i, iter, max_iter
144 real(kind=
rp) :: rnorm, rtr, rtr0, rtz2, rtz1
145 real(kind=
rp) :: beta, pap, alpha, alphm, norm_fac
147 if (
present(niter))
then
150 max_iter = this%max_iter
152 norm_fac = one / sqrt(coef%volume)
156 x%x(i,1,1,1) = 0.0_rp
161 rtr =
glsc3(this%r, coef%mult, this%r, n)
162 rnorm = sqrt(rtr)*norm_fac
163 ksp_results%res_start = rnorm
164 ksp_results%res_final = rnorm
166 if(
abscmp(rnorm, zero))
return
167 call this%monitor_start(
'CG')
168 do iter = 1, max_iter
169 call this%M%solve(this%z, this%r, n)
171 rtz1 =
glsc3(this%r, coef%mult, this%z, n)
174 if (iter .eq. 1) beta = zero
175 call add2s1(this%p, this%z, beta, n)
177 call ax%compute(this%w, this%p, coef, x%msh, x%Xh)
178 call gs_h%op(this%w, n, gs_op_add)
181 pap =
glsc3(this%w, coef%mult, this%p, n)
186 x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * this%p(i)
187 this%r(i) = this%r(i) + alphm * this%w(i)
190 rtr =
glsc3(this%r, coef%mult, this%r, n)
191 if (iter .eq. 1) rtr0 = rtr
192 rnorm = sqrt(rtr) * norm_fac
193 call this%monitor_iter(iter, rnorm)
194 if (rnorm .lt. this%abs_tol)
then
198 call this%monitor_stop()
199 ksp_results%res_final = rnorm
200 ksp_results%iter = iter
205 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
206 class(
sx_cg_t),
intent(inout) :: this
207 class(
ax_t),
intent(inout) :: ax
208 type(
field_t),
intent(inout) :: x
209 type(
field_t),
intent(inout) :: y
210 type(
field_t),
intent(inout) :: z
211 integer,
intent(in) :: n
212 real(kind=
rp),
dimension(n),
intent(inout) :: fx
213 real(kind=
rp),
dimension(n),
intent(inout) :: fy
214 real(kind=
rp),
dimension(n),
intent(inout) :: fz
215 type(
coef_t),
intent(inout) :: coef
219 type(
gs_t),
intent(inout) :: gs_h
221 integer,
optional,
intent(in) :: niter
223 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
224 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
225 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
Defines a Matrix-vector product.
Defines a boundary condition.
Defines various Conjugate Gradient methods.
subroutine sx_cg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a standard PCG solver.
subroutine sx_cg_free(this)
Deallocate a standard PCG solver.
type(ksp_monitor_t) function, dimension(3) sx_cg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard PCG coupled solve.
type(ksp_monitor_t) function sx_cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
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 add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
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 (SX version)
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.