51 real(kind=
rp),
allocatable :: w(:)
52 real(kind=
rp),
allocatable :: r(:)
53 real(kind=
rp),
allocatable :: z(:,:)
54 real(kind=
rp),
allocatable :: v(:,:)
56 real(kind=
xp),
allocatable :: h(:,:)
57 real(kind=
xp),
allocatable :: s(:)
58 real(kind=
xp),
allocatable :: gam(:)
59 real(kind=
xp),
allocatable :: c(:)
71 rel_tol, abs_tol, monitor)
72 class(
gmres_t),
intent(inout) :: this
73 integer,
intent(in) :: n
74 integer,
intent(in) :: max_iter
75 class(
pc_t),
optional,
intent(in),
target :: M
76 integer,
optional,
intent(in) :: lgmres
77 real(kind=
rp),
optional,
intent(in) :: rel_tol
78 real(kind=
rp),
optional,
intent(in) :: abs_tol
79 logical,
optional,
intent(in) :: monitor
81 if (
present(lgmres))
then
97 allocate(this%c(this%lgmres))
98 allocate(this%s(this%lgmres))
99 allocate(this%gam(this%lgmres + 1))
101 allocate(this%z(n, this%lgmres))
102 allocate(this%v(n, this%lgmres))
104 allocate(this%h(this%lgmres, this%lgmres))
106 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
107 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
108 else if (
present(rel_tol) .and.
present(abs_tol))
then
109 call this%ksp_init(max_iter, rel_tol, abs_tol)
110 else if (
present(monitor) .and.
present(abs_tol))
then
111 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
112 else if (
present(rel_tol) .and.
present(monitor))
then
113 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
114 else if (
present(rel_tol))
then
115 call this%ksp_init(max_iter, rel_tol = rel_tol)
116 else if (
present(abs_tol))
then
117 call this%ksp_init(max_iter, abs_tol = abs_tol)
118 else if (
present(monitor))
then
119 call this%ksp_init(max_iter, monitor = monitor)
121 call this%ksp_init(max_iter)
128 class(
gmres_t),
intent(inout) :: this
132 if (
allocated(this%w))
then
136 if (
allocated(this%c))
then
140 if (
allocated(this%r))
then
144 if (
allocated(this%z))
then
148 if (
allocated(this%h))
then
152 if (
allocated(this%v))
then
156 if (
allocated(this%s))
then
161 if (
allocated(this%gam))
then
170 function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
172 class(
gmres_t),
intent(inout) :: this
173 class(
ax_t),
intent(in) :: ax
174 type(
field_t),
intent(inout) :: x
175 integer,
intent(in) :: n
176 real(kind=
rp),
dimension(n),
intent(in) :: f
177 type(
coef_t),
intent(inout) :: coef
179 type(
gs_t),
intent(inout) :: gs_h
181 integer,
optional,
intent(in) :: niter
182 integer :: iter, max_iter
183 integer :: i, j, k, l, ierr
184 real(kind=
xp) :: w_plus(neko_blk_size), x_plus(neko_blk_size)
185 real(kind=
xp) :: alpha, lr, alpha2, norm_fac
186 real(kind=
rp) :: temp, rnorm
192 if (
present(niter))
then
195 max_iter = this%max_iter
198 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
199 v => this%v, s => this%s, gam => this%gam)
201 norm_fac = 1.0_rp / sqrt(coef%volume)
207 call this%monitor_start(
'GMRES')
208 do while (.not. conv .and. iter .lt. max_iter)
210 if (iter .eq. 0)
then
214 call ax%compute(w, x%x, coef, x%msh, x%Xh)
215 call gs_h%op(w, n, gs_op_add)
220 gam(1) = sqrt(
glsc3(r, r, coef%mult, n))
221 if (iter .eq. 0)
then
222 ksp_results%res_start = gam(1) * norm_fac
225 if (
abscmp(gam(1), 0.0_xp))
return
228 temp = 1.0_rp / gam(1)
229 call cmult2(v(1,1), r, temp, n)
230 do j = 1, this%lgmres
233 call this%M%solve(z(1,j), v(1,j), n)
235 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
236 call gs_h%op(w, n, gs_op_add)
243 do i = 0, n, neko_blk_size
244 if (i + neko_blk_size .le. n)
then
246 do k = 1, neko_blk_size
248 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
255 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
261 call mpi_allreduce(mpi_in_place, h(1,j), j, &
265 do i = 0, n, neko_blk_size
266 if (i + neko_blk_size .le. n)
then
267 do k = 1, neko_blk_size
271 do k = 1, neko_blk_size
272 w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
275 do k = 1, neko_blk_size
276 w(i+k) = w(i+k) + w_plus(k)
277 alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
283 w_plus(1) = w_plus(1) - h(l,j) * v(i+k,l)
285 w(i+k) = w(i+k) + w_plus(1)
286 alpha2 = alpha2 + (w(i+k)**2) * coef%mult(i+k,1,1,1)
291 call mpi_allreduce(mpi_in_place,alpha2, 1, &
296 h(i,j) = c(i)*temp + s(i) * h(i+1,j)
297 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
301 if (
abscmp(alpha, 0.0_xp))
then
306 lr = sqrt(h(j,j) * h(j,j) + alpha2)
311 gam(j+1) = -s(j) * gam(j)
312 gam(j) = c(j) * gam(j)
313 rnorm = abs(gam(j+1)) * norm_fac
314 call this%monitor_iter(iter, rnorm)
315 if (rnorm .lt. this%abs_tol)
then
320 if (iter + 1 .gt. max_iter)
exit
322 if (j .lt. this%lgmres)
then
323 temp = 1.0_rp / alpha
324 call cmult2(v(1,j+1), w, temp, n)
329 j = min(j, this%lgmres)
333 temp = temp - h(k,i) * c(i)
338 do i = 0, n, neko_blk_size
339 if (i + neko_blk_size .le. n)
then
340 do k = 1, neko_blk_size
344 do k = 1, neko_blk_size
345 x_plus(k) = x_plus(k) + c(l) * z(i+k,l)
348 do k = 1, neko_blk_size
349 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
355 x_plus(1) = x_plus(1) + c(l) * z(i+k,l)
357 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
364 call this%monitor_stop()
365 ksp_results%res_final = rnorm
366 ksp_results%iter = iter
372 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
373 class(
gmres_t),
intent(inout) :: this
374 class(ax_t),
intent(in) :: ax
375 type(field_t),
intent(inout) :: x
376 type(field_t),
intent(inout) :: y
377 type(field_t),
intent(inout) :: z
378 integer,
intent(in) :: n
379 real(kind=rp),
dimension(n),
intent(in) :: fx
380 real(kind=rp),
dimension(n),
intent(in) :: fy
381 real(kind=rp),
dimension(n),
intent(in) :: fz
382 type(coef_t),
intent(inout) :: coef
383 type(bc_list_t),
intent(in) :: blstx
384 type(bc_list_t),
intent(in) :: blsty
385 type(bc_list_t),
intent(in) :: blstz
386 type(gs_t),
intent(inout) :: gs_h
387 type(ksp_monitor_t),
dimension(3) :: ksp_results
388 integer,
optional,
intent(in) :: niter
390 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
391 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
392 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
Defines a Matrix-vector product.
Defines a boundary condition.
type(mpi_comm) neko_comm
MPI communicator.
type(mpi_datatype) mpi_extra_precision
Defines various GMRES methods.
type(ksp_monitor_t) function, dimension(3) gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard GMRES coupled solve.
subroutine gmres_free(this)
Deallocate a standard GMRES solver.
subroutine gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol, monitor)
Initialise a standard GMRES solver.
type(ksp_monitor_t) function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard GMRES solve.
Implements the base abstract type for Krylov solvers plus helper types.
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
subroutine, public rone(a, n)
Set all elements to one.
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
subroutine, public sub2(a, b, n)
Vector substraction .
integer, parameter, public xp
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Standard preconditioned generalized minimal residual method.
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.