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(:)
59 real(kind=
rp),
allocatable :: wk1(:)
69 subroutine gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol)
70 class(
gmres_t),
intent(inout) :: this
71 integer,
intent(in) :: n
72 integer,
intent(in) :: max_iter
73 class(
pc_t),
optional,
intent(inout),
target :: M
74 integer,
optional,
intent(inout) :: lgmres
75 real(kind=
rp),
optional,
intent(inout) :: rel_tol
76 real(kind=
rp),
optional,
intent(inout) :: abs_tol
78 if (
present(lgmres))
then
95 allocate(this%c(this%lgmres))
96 allocate(this%s(this%lgmres))
97 allocate(this%gam(this%lgmres + 1))
99 allocate(this%z(n,this%lgmres))
100 allocate(this%v(n,this%lgmres))
102 allocate(this%h(this%lgmres,this%lgmres))
105 if (
present(rel_tol) .and.
present(abs_tol))
then
106 call this%ksp_init(max_iter, rel_tol, abs_tol)
107 else if (
present(rel_tol))
then
108 call this%ksp_init(max_iter, rel_tol=rel_tol)
109 else if (
present(abs_tol))
then
110 call this%ksp_init(max_iter, abs_tol=abs_tol)
112 call this%ksp_init(max_iter)
119 class(
gmres_t),
intent(inout) :: this
123 if (
allocated(this%w))
then
127 if (
allocated(this%c))
then
131 if (
allocated(this%r))
then
135 if (
allocated(this%z))
then
139 if (
allocated(this%h))
then
143 if (
allocated(this%v))
then
147 if (
allocated(this%s))
then
152 if (
allocated(this%gam))
then
156 if (
allocated(this%wk1))
then
165 function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
166 class(
gmres_t),
intent(inout) :: this
167 class(
ax_t),
intent(inout) :: ax
168 type(
field_t),
intent(inout) :: x
169 integer,
intent(in) :: n
170 real(kind=
rp),
dimension(n),
intent(inout) :: f
171 type(
coef_t),
intent(inout) :: coef
173 type(
gs_t),
intent(inout) :: gs_h
175 integer,
optional,
intent(in) :: niter
176 integer :: iter, max_iter
177 integer :: i, j, k, l, ierr
178 real(kind=
rp) :: w_plus(neko_blk_size), x_plus(neko_blk_size)
179 real(kind=
rp) :: rnorm, alpha, temp, lr, alpha2, norm_fac
185 if (
present(niter))
then
188 max_iter = this%max_iter
191 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
192 v => this%v, s => this%s, gam => this%gam, wk1 =>this%wk1)
194 norm_fac = 1.0_rp / sqrt(coef%volume)
196 call rzero(gam, this%lgmres + 1)
197 call rone(s, this%lgmres)
198 call rone(c, this%lgmres)
199 call rzero(h, this%lgmres * this%lgmres)
200 do while (.not. conv .and. iter .lt. max_iter)
206 call ax%compute(w, x%x, coef, x%msh, x%Xh)
207 call gs_h%op(w, n, gs_op_add)
212 gam(1) = sqrt(
glsc3(r, r, coef%mult, n))
214 ksp_results%res_start = gam(1) * norm_fac
217 if (
abscmp(gam(1), 0.0_rp))
return
220 temp = 1.0_rp / gam(1)
221 call cmult2(v(1,1), r, temp, n)
222 do j = 1, this%lgmres
225 call this%M%solve(z(1,j), v(1,j), n)
227 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
228 call gs_h%op(w, n, gs_op_add)
235 do i = 0, n, neko_blk_size
236 if (i + neko_blk_size .le. n)
then
238 do k = 1, neko_blk_size
240 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
247 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
253 call mpi_allreduce(h(1,j), wk1, j, &
255 call copy(h(1,j), wk1, j)
258 do i = 0,n,neko_blk_size
259 if (i + neko_blk_size .le. n)
then
260 do k = 1, neko_blk_size
264 do k = 1, neko_blk_size
265 w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
268 do k = 1, neko_blk_size
269 w(i+k) = w(i+k) + w_plus(k)
270 alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
276 w_plus(1) = w_plus(1) - h(l,j) * v(i+k,l)
278 w(i+k) = w(i+k) + w_plus(1)
279 alpha2 = alpha2 + (w(i+k)**2) * coef%mult(i+k,1,1,1)
284 call mpi_allreduce(alpha2, temp, 1, &
290 h(i ,j) = c(i)*temp + s(i) * h(i+1,j)
291 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
295 if(
abscmp(alpha, 0.0_rp))
then
300 lr = sqrt(h(j,j) * h(j,j) + alpha**2)
305 gam(j+1) = -s(j) * gam(j)
306 gam(j) = c(j) * gam(j)
308 rnorm = abs(gam(j+1)) * norm_fac
309 if (rnorm .lt. this%abs_tol)
then
314 if (iter + 1 .gt. max_iter)
exit
316 if( j .lt. this%lgmres)
then
317 temp = 1.0_rp / alpha
318 call cmult2(v(1,j+1), w, temp, n)
323 j = min(j, this%lgmres)
327 temp = temp - h(k,i) * c(i)
332 do i = 0, n, neko_blk_size
333 if (i + neko_blk_size .le. n)
then
334 do k = 1, neko_blk_size
338 do k = 1, neko_blk_size
339 x_plus(k) = x_plus(k) + c(l) * z(i+k,l)
342 do k = 1, neko_blk_size
343 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
349 x_plus(1) = x_plus(1) + c(l) * z(i+k,l)
351 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
359 ksp_results%res_final = rnorm
360 ksp_results%iter = iter
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.
subroutine gmres_free(this)
Deallocate a standard GMRES solver.
subroutine gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol)
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.