51 real(kind=
rp),
allocatable :: p(:)
52 real(kind=
rp),
allocatable :: p_hat(:)
53 real(kind=
rp),
allocatable :: r(:)
54 real(kind=
rp),
allocatable :: s(:)
55 real(kind=
rp),
allocatable :: s_hat(:)
56 real(kind=
rp),
allocatable :: t(:)
57 real(kind=
rp),
allocatable :: v(:)
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
82 allocate(this%p_hat(n))
85 allocate(this%s_hat(n))
92 if (
present(rel_tol) .and.
present(abs_tol))
then
93 call this%ksp_init(max_iter, rel_tol, abs_tol)
94 else if (
present(rel_tol))
then
95 call this%ksp_init(max_iter, rel_tol=rel_tol)
96 else if (
present(abs_tol))
then
97 call this%ksp_init(max_iter, abs_tol=abs_tol)
99 call this%ksp_init(max_iter)
110 if (
allocated(this%v))
then
114 if (
allocated(this%r))
then
118 if (
allocated(this%t))
then
122 if (
allocated(this%p))
then
126 if (
allocated(this%p_hat))
then
127 deallocate(this%p_hat)
130 if (
allocated(this%s))
then
134 if (
allocated(this%s_hat))
then
135 deallocate(this%s_hat)
144 function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
146 class(
ax_t),
intent(inout) :: ax
147 type(
field_t),
intent(inout) :: x
148 integer,
intent(in) :: n
149 real(kind=
rp),
dimension(n),
intent(inout) :: f
150 type(
coef_t),
intent(inout) :: coef
152 type(
gs_t),
intent(inout) :: gs_h
154 integer,
optional,
intent(in) :: niter
155 integer :: iter, max_iter
156 real(kind=
rp) :: rnorm, rtr, norm_fac, gamma
157 real(kind=
rp) :: beta, alpha, omega, rho_1, rho_2
159 if (
present(niter))
then
162 max_iter = this%max_iter
164 norm_fac = 1.0_rp / sqrt(coef%volume)
166 associate(r => this%r, t => this%t, s => this%s, v => this%v, p => this%p, &
167 s_hat => this%s_hat, p_hat => this%p_hat)
172 rtr = sqrt(
glsc3(r, coef%mult, r, n))
173 rnorm = rtr * norm_fac
174 gamma = rnorm * this%rel_tol
175 ksp_results%res_start = rnorm
176 ksp_results%res_final = rnorm
178 if(
abscmp(rnorm, 0.0_rp))
return
179 do iter = 1, max_iter
181 rho_1 =
glsc3(r, coef%mult, f ,n)
187 if (iter .eq. 1)
then
190 beta = (rho_1 / rho_2) * (alpha / omega)
191 call p_update(p, r, v, beta, omega, n)
194 call this%M%solve(p_hat, p, n)
195 call ax%compute(v, p_hat, coef, x%msh, x%Xh)
196 call gs_h%op(v, n, gs_op_add)
198 alpha = rho_1 /
glsc3(f, coef%mult, v, n)
200 call add2s2(s, v, -alpha, n)
201 rtr =
glsc3(s, coef%mult, s, n)
202 rnorm = sqrt(rtr) * norm_fac
203 if (rnorm .lt. this%abs_tol .or. rnorm .lt. gamma)
then
204 call add2s2(x%x, p_hat, alpha,n)
208 call this%M%solve(s_hat, s, n)
209 call ax%compute(t, s_hat, coef, x%msh, x%Xh)
210 call gs_h%op(t, n, gs_op_add)
212 omega =
glsc3(t, coef%mult, s, n) &
213 /
glsc3(t, coef%mult, t, n)
214 call x_update(x%x, p_hat, s_hat, alpha, omega, n)
216 call add2s2(r, t, -omega, n)
218 rtr =
glsc3(r, coef%mult, r, n)
219 rnorm = sqrt(rtr) * norm_fac
220 if (rnorm .lt. this%abs_tol .or. rnorm .lt. gamma)
then
231 ksp_results%res_final = rnorm
232 ksp_results%iter = iter
Defines a Matrix-vector product.
Defines a boundary condition.
Defines various Bi-Conjugate Gradient Stabilized methods.
subroutine bicgstab_init(this, n, max_iter, M, rel_tol, abs_tol)
Constructor.
subroutine bicgstab_free(this)
Deallocate a standard BiCGSTAB solver.
type(ksp_monitor_t) function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Bi-Conjugate Gradient Stabilized method 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 x_update(a, b, c, c1, c2, n)
Returns .
subroutine, public copy(a, b, n)
Copy a vector .
real(kind=rp), parameter, public neko_eps
Machine epsilon .
subroutine, public rzero(a, n)
Zero a real vector.
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
subroutine, public p_update(a, b, c, c1, c2, n)
Returns .
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Standard preconditioned Bi-Conjugate Gradient Stabilized 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.