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 :: ml(:)
57 real(kind=
rp),
allocatable :: v(:,:)
58 real(kind=
rp),
allocatable :: s(:)
59 real(kind=
rp),
allocatable :: mu(:)
60 real(kind=
rp),
allocatable :: gam(:)
61 real(kind=
rp),
allocatable :: wk1(:)
62 real(kind=
rp) :: rnorm
74 rel_tol, abs_tol, monitor)
76 integer,
intent(in) :: n
77 integer,
intent(in) :: max_iter
78 class(
pc_t),
optional,
intent(in),
target :: M
79 integer,
optional,
intent(in) :: lgmres
80 real(kind=
rp),
optional,
intent(in) :: rel_tol
81 real(kind=
rp),
optional,
intent(in) :: abs_tol
82 logical,
optional,
intent(in) :: monitor
84 if (
present(lgmres))
then
101 allocate(this%wk1(n))
103 allocate(this%c(this%lgmres))
104 allocate(this%s(this%lgmres))
105 allocate(this%gam(this%lgmres + 1))
107 allocate(this%z(n,this%lgmres))
108 allocate(this%v(n,this%lgmres))
110 allocate(this%h(this%lgmres,this%lgmres))
113 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
114 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
115 else if (
present(rel_tol) .and.
present(abs_tol))
then
116 call this%ksp_init(max_iter, rel_tol, abs_tol)
117 else if (
present(monitor) .and.
present(abs_tol))
then
118 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
119 else if (
present(rel_tol) .and.
present(monitor))
then
120 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
121 else if (
present(rel_tol))
then
122 call this%ksp_init(max_iter, rel_tol = rel_tol)
123 else if (
present(abs_tol))
then
124 call this%ksp_init(max_iter, abs_tol = abs_tol)
125 else if (
present(monitor))
then
126 call this%ksp_init(max_iter, monitor = monitor)
128 call this%ksp_init(max_iter)
139 if (
allocated(this%w))
then
143 if (
allocated(this%c))
then
147 if (
allocated(this%r))
then
151 if (
allocated(this%z))
then
155 if (
allocated(this%h))
then
159 if (
allocated(this%ml))
then
163 if (
allocated(this%v))
then
167 if (
allocated(this%s))
then
171 if (
allocated(this%mu))
then
175 if (
allocated(this%gam))
then
179 if (
allocated(this%wk1))
then
188 function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
190 class(
ax_t),
intent(in) :: ax
191 type(
field_t),
intent(inout) :: x
192 integer,
intent(in) :: n
193 real(kind=
rp),
dimension(n),
intent(in) :: f
194 type(
coef_t),
intent(inout) :: coef
196 type(
gs_t),
intent(inout) :: gs_h
198 integer,
optional,
intent(in) :: niter
199 integer :: iter, max_iter, glb_n
200 integer :: i, j, k, ierr
201 real(kind=
rp),
parameter :: one = 1.0
202 real(kind=
rp) :: rnorm
203 real(kind=
rp) :: alpha, temp, l
204 real(kind=
rp) :: ratio, div0, norm_fac
210 glb_n = n / x%msh%nelv * x%msh%glb_nelv
212 if (
present(niter))
then
215 max_iter = this%max_iter
218 call rone(this%ml, n)
219 call rone(this%mu, n)
220 norm_fac = one / sqrt(coef%volume)
222 call rzero(this%gam, this%lgmres + 1)
223 call rone(this%s, this%lgmres)
224 call rone(this%c, this%lgmres)
225 call rzero(this%h, this%lgmres * this%lgmres)
227 call this%monitor_start(
'GMRES')
228 do while (.not. conv .and. iter .lt. max_iter)
232 call col3(this%r,this%ml,f,n)
235 call copy (this%r,f,n)
236 call ax%compute(this%w, x%x, coef, x%msh, x%Xh)
237 call gs_h%op(this%w, n, gs_op_add)
239 call add2s2(this%r,this%w,-one,n)
240 call col2(this%r,this%ml,n)
242 this%gam(1) = sqrt(
glsc3(this%r, this%r, coef%mult, n))
244 div0 = this%gam(1) * norm_fac
245 ksp_results%res_start = div0
248 if (
abscmp(this%gam(1), 0.0_rp))
return
251 temp = one / this%gam(1)
252 call cmult2(this%v(1,1), this%r, temp, n)
253 do j = 1, this%lgmres
255 call col3(this%w, this%mu, this%v(1,j), n)
258 call this%M%solve(this%z(1,j), this%w, n)
260 call ax%compute(this%w, this%z(1,j), coef, x%msh, x%Xh)
261 call gs_h%op(this%w, n, gs_op_add)
263 call col2(this%w, this%ml, n)
268 this%h(i,j) = this%h(i,j) + &
269 this%w(k) * this%v(k,i) * coef%mult(k,1,1,1)
274 call mpi_allreduce(this%h(1,j), this%wk1, j, &
276 call copy(this%h(1,j), this%wk1, j)
280 this%w(k) = this%w(k) - this%h(i,j) * this%v(k,i)
287 this%h(i ,j) = this%c(i)*temp + this%s(i)*this%h(i+1,j)
288 this%h(i+1,j) = -this%s(i)*temp + this%c(i)*this%h(i+1,j)
291 alpha = sqrt(
glsc3(this%w, this%w, coef%mult, n))
293 if(
abscmp(alpha, 0.0_rp))
then
297 l = sqrt(this%h(j,j) * this%h(j,j) + alpha**2)
299 this%c(j) = this%h(j,j) * temp
300 this%s(j) = alpha * temp
302 this%gam(j+1) = -this%s(j) * this%gam(j)
303 this%gam(j) = this%c(j) * this%gam(j)
305 rnorm = abs(this%gam(j+1)) * norm_fac
306 call this%monitor_iter(iter, rnorm)
308 if (rnorm .lt. this%abs_tol)
then
313 if (iter + 1 .gt. max_iter)
exit
315 if( j .lt. this%lgmres)
then
317 call cmult2(this%v(1,j+1), this%w, temp, n)
320 j = min(j, this%lgmres)
325 temp = temp - this%h(k,i) * this%c(i)
327 this%c(k) = temp / this%h(k,k)
332 x%x(k,1,1,1) = x%x(k,1,1,1) + this%c(i) * this%z(k,i)
336 call this%monitor_stop()
337 ksp_results%res_final = rnorm
338 ksp_results%iter = iter
343 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
345 class(
ax_t),
intent(in) :: ax
346 type(
field_t),
intent(inout) :: x
347 type(
field_t),
intent(inout) :: y
348 type(
field_t),
intent(inout) :: z
349 integer,
intent(in) :: n
350 real(kind=
rp),
dimension(n),
intent(in) :: fx
351 real(kind=
rp),
dimension(n),
intent(in) :: fy
352 real(kind=
rp),
dimension(n),
intent(in) :: fz
353 type(
coef_t),
intent(inout) :: coef
357 type(
gs_t),
intent(inout) :: gs_h
359 integer,
optional,
intent(in) :: niter
361 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
362 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
363 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.
subroutine sx_gmres_free(this)
Deallocate a standard GMRES solver.
subroutine sx_gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol, monitor)
Initialise a standard GMRES solver.
type(ksp_monitor_t) function, dimension(3) sx_gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard GMRES coupled solve.
type(ksp_monitor_t) function sx_gmres_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.
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 col2(a, b, n)
Vector multiplication .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
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)
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 (SX version)
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.