49 real(kind=
rp),
allocatable :: w(:)
50 real(kind=
rp),
allocatable :: r(:)
51 real(kind=
rp),
allocatable :: p(:)
52 real(kind=
rp),
allocatable :: z(:)
62 subroutine sx_cg_init(this, n, max_iter, M, rel_tol, abs_tol)
63 class(
sx_cg_t),
intent(inout) :: this
64 class(
pc_t),
optional,
intent(inout),
target :: M
65 integer,
intent(in) :: n
66 integer,
intent(in) :: max_iter
67 real(kind=
rp),
optional,
intent(inout) :: rel_tol
68 real(kind=
rp),
optional,
intent(inout) :: abs_tol
81 if (
present(rel_tol) .and.
present(abs_tol))
then
82 call this%ksp_init(max_iter, rel_tol, abs_tol)
83 else if (
present(rel_tol))
then
84 call this%ksp_init(max_iter, rel_tol=rel_tol)
85 else if (
present(abs_tol))
then
86 call this%ksp_init(max_iter, abs_tol=abs_tol)
88 call this%ksp_init(max_iter)
95 class(
sx_cg_t),
intent(inout) :: this
99 if (
allocated(this%w))
then
103 if (
allocated(this%r))
then
107 if (
allocated(this%p))
then
111 if (
allocated(this%z))
then
120 function sx_cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
121 class(
sx_cg_t),
intent(inout) :: this
122 class(
ax_t),
intent(inout) :: ax
123 type(
field_t),
intent(inout) :: x
124 integer,
intent(in) :: n
125 real(kind=
rp),
dimension(n),
intent(inout) :: f
126 type(
coef_t),
intent(inout) :: coef
128 type(
gs_t),
intent(inout) :: gs_h
130 integer,
optional,
intent(in) :: niter
131 real(kind=
rp),
parameter :: one = 1.0
132 real(kind=
rp),
parameter :: zero = 0.0
133 integer :: i, iter, max_iter
134 real(kind=
rp) :: rnorm, rtr, rtr0, rtz2, rtz1
135 real(kind=
rp) :: beta, pap, alpha, alphm, norm_fac
137 if (
present(niter))
then
140 max_iter = this%max_iter
142 norm_fac = one / sqrt(coef%volume)
146 x%x(i,1,1,1) = 0.0_rp
151 rtr =
glsc3(this%r, coef%mult, this%r, n)
152 rnorm = sqrt(rtr)*norm_fac
153 ksp_results%res_start = rnorm
154 ksp_results%res_final = rnorm
156 if(
abscmp(rnorm, zero))
return
158 do iter = 1, max_iter
159 call this%M%solve(this%z, this%r, n)
161 rtz1 =
glsc3(this%r, coef%mult, this%z, n)
164 if (iter .eq. 1) beta = zero
165 call add2s1(this%p, this%z, beta, n)
167 call ax%compute(this%w, this%p, coef, x%msh, x%Xh)
168 call gs_h%op(this%w, n, gs_op_add)
171 pap =
glsc3(this%w, coef%mult, this%p, n)
176 x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * this%p(i)
177 this%r(i) = this%r(i) + alphm * this%w(i)
180 rtr =
glsc3(this%r, coef%mult, this%r, n)
181 if (iter .eq. 1) rtr0 = rtr
182 rnorm = sqrt(rtr) * norm_fac
183 if (rnorm .lt. this%abs_tol)
then
187 ksp_results%res_final = rnorm
188 ksp_results%iter = iter
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)
Initialise a standard PCG solver.
subroutine sx_cg_free(this)
Deallocate a standard PCG solver.
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.