51 real(kind=
rp),
allocatable :: w(:)
52 real(kind=
rp),
allocatable :: c(:)
53 real(kind=
rp),
allocatable :: r(:)
54 real(kind=
rp),
allocatable :: z(:,:)
55 real(kind=
rp),
allocatable :: h(:,:)
56 real(kind=
rp),
allocatable :: v(:,:)
57 real(kind=
rp),
allocatable :: s(:)
58 real(kind=
rp),
allocatable :: gam(:)
70 rel_tol, abs_tol, monitor)
71 class(
gmres_t),
intent(inout) :: this
72 integer,
intent(in) :: n
73 integer,
intent(in) :: max_iter
74 class(
pc_t),
optional,
intent(inout),
target :: M
75 integer,
optional,
intent(inout) :: lgmres
76 real(kind=
rp),
optional,
intent(inout) :: rel_tol
77 real(kind=
rp),
optional,
intent(inout) :: abs_tol
78 logical,
optional,
intent(in) :: monitor
80 if (
present(lgmres))
then
96 allocate(this%c(this%lgmres))
97 allocate(this%s(this%lgmres))
98 allocate(this%gam(this%lgmres + 1))
100 allocate(this%z(n, this%lgmres))
101 allocate(this%v(n, this%lgmres))
103 allocate(this%h(this%lgmres, this%lgmres))
105 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
106 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
107 else if (
present(rel_tol) .and.
present(abs_tol))
then
108 call this%ksp_init(max_iter, rel_tol, abs_tol)
109 else if (
present(monitor) .and.
present(abs_tol))
then
110 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
111 else if (
present(rel_tol) .and.
present(monitor))
then
112 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
113 else if (
present(rel_tol))
then
114 call this%ksp_init(max_iter, rel_tol = rel_tol)
115 else if (
present(abs_tol))
then
116 call this%ksp_init(max_iter, abs_tol = abs_tol)
117 else if (
present(monitor))
then
118 call this%ksp_init(max_iter, monitor = monitor)
120 call this%ksp_init(max_iter)
127 class(
gmres_t),
intent(inout) :: this
131 if (
allocated(this%w))
then
135 if (
allocated(this%c))
then
139 if (
allocated(this%r))
then
143 if (
allocated(this%z))
then
147 if (
allocated(this%h))
then
151 if (
allocated(this%v))
then
155 if (
allocated(this%s))
then
160 if (
allocated(this%gam))
then
169 function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
171 class(
gmres_t),
intent(inout) :: this
172 class(
ax_t),
intent(inout) :: ax
173 type(
field_t),
intent(inout) :: x
174 integer,
intent(in) :: n
175 real(kind=
rp),
dimension(n),
intent(inout) :: f
176 type(
coef_t),
intent(inout) :: coef
178 type(
gs_t),
intent(inout) :: gs_h
180 integer,
optional,
intent(in) :: niter
181 integer :: iter, max_iter
182 integer :: i, j, k, l, ierr
183 real(kind=
rp) :: w_plus(neko_blk_size), x_plus(neko_blk_size)
184 real(kind=
rp) :: rnorm, alpha, temp, lr, alpha2, norm_fac
190 if (
present(niter))
then
193 max_iter = this%max_iter
196 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
197 v => this%v, s => this%s, gam => this%gam)
199 norm_fac = 1.0_rp / sqrt(coef%volume)
201 call rzero(gam, this%lgmres + 1)
202 call rone(s, this%lgmres)
203 call rone(c, this%lgmres)
204 call rzero(h, this%lgmres * this%lgmres)
205 call this%monitor_start(
'GMRES')
206 do while (.not. conv .and. iter .lt. max_iter)
208 if (iter .eq. 0)
then
212 call ax%compute(w, x%x, coef, x%msh, x%Xh)
213 call gs_h%op(w, n, gs_op_add)
218 gam(1) = sqrt(
glsc3(r, r, coef%mult, n))
219 if (iter .eq. 0)
then
220 ksp_results%res_start = gam(1) * norm_fac
223 if (
abscmp(gam(1), 0.0_rp))
return
226 temp = 1.0_rp / gam(1)
227 call cmult2(v(1,1), r, temp, n)
228 do j = 1, this%lgmres
231 call this%M%solve(z(1,j), v(1,j), n)
233 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
234 call gs_h%op(w, n, gs_op_add)
241 do i = 0, n, neko_blk_size
242 if (i + neko_blk_size .le. n)
then
244 do k = 1, neko_blk_size
246 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
253 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
259 call mpi_allreduce(mpi_in_place, h(1,j), j, &
263 do i = 0, n, neko_blk_size
264 if (i + neko_blk_size .le. n)
then
265 do k = 1, neko_blk_size
269 do k = 1, neko_blk_size
270 w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
273 do k = 1, neko_blk_size
274 w(i+k) = w(i+k) + w_plus(k)
275 alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
281 w_plus(1) = w_plus(1) - h(l,j) * v(i+k,l)
283 w(i+k) = w(i+k) + w_plus(1)
284 alpha2 = alpha2 + (w(i+k)**2) * coef%mult(i+k,1,1,1)
289 call mpi_allreduce(alpha2, temp, 1, &
295 h(i,j) = c(i)*temp + s(i) * h(i+1,j)
296 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
300 if (
abscmp(alpha, 0.0_rp))
then
305 lr = sqrt(h(j,j) * h(j,j) + alpha2)
310 gam(j+1) = -s(j) * gam(j)
311 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(inout) :: 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(inout) :: fx
380 real(kind=rp),
dimension(n),
intent(inout) :: fy
381 real(kind=rp),
dimension(n),
intent(inout) :: fz
382 type(coef_t),
intent(inout) :: coef
383 type(bc_list_t),
intent(inout) :: blstx
384 type(bc_list_t),
intent(inout) :: blsty
385 type(bc_list_t),
intent(inout) :: 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_real_precision
MPI type for working precision of REAL types.
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 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.