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
72 subroutine sx_gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol)
74 integer,
intent(in) :: n
75 integer,
intent(in) :: max_iter
76 class(
pc_t),
optional,
intent(inout),
target :: M
77 integer,
optional,
intent(inout) :: lgmres
78 real(kind=
rp),
optional,
intent(inout) :: rel_tol
79 real(kind=
rp),
optional,
intent(inout) :: abs_tol
81 if (
present(lgmres))
then
100 allocate(this%c(this%lgmres))
101 allocate(this%s(this%lgmres))
102 allocate(this%gam(this%lgmres + 1))
104 allocate(this%z(n,this%lgmres))
105 allocate(this%v(n,this%lgmres))
107 allocate(this%h(this%lgmres,this%lgmres))
110 if (
present(rel_tol) .and.
present(abs_tol))
then
111 call this%ksp_init(max_iter, rel_tol, abs_tol)
112 else if (
present(rel_tol))
then
113 call this%ksp_init(max_iter, rel_tol=rel_tol)
114 else if (
present(abs_tol))
then
115 call this%ksp_init(max_iter, abs_tol=abs_tol)
117 call this%ksp_init(max_iter)
128 if (
allocated(this%w))
then
132 if (
allocated(this%c))
then
136 if (
allocated(this%r))
then
140 if (
allocated(this%z))
then
144 if (
allocated(this%h))
then
148 if (
allocated(this%ml))
then
152 if (
allocated(this%v))
then
156 if (
allocated(this%s))
then
160 if (
allocated(this%mu))
then
164 if (
allocated(this%gam))
then
168 if (
allocated(this%wk1))
then
177 function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
179 class(
ax_t),
intent(inout) :: ax
180 type(
field_t),
intent(inout) :: x
181 integer,
intent(in) :: n
182 real(kind=
rp),
dimension(n),
intent(inout) :: f
183 type(
coef_t),
intent(inout) :: coef
185 type(
gs_t),
intent(inout) :: gs_h
187 integer,
optional,
intent(in) :: niter
188 integer :: iter, max_iter, glb_n
189 integer :: i, j, k, ierr
190 real(kind=
rp),
parameter :: one = 1.0
191 real(kind=
rp) :: rnorm
192 real(kind=
rp) :: alpha, temp, l
193 real(kind=
rp) :: ratio, div0, norm_fac
199 glb_n = n / x%msh%nelv * x%msh%glb_nelv
201 if (
present(niter))
then
204 max_iter = this%max_iter
207 call rone(this%ml, n)
208 call rone(this%mu, n)
209 norm_fac = one / sqrt(coef%volume)
211 call rzero(this%gam, this%lgmres + 1)
212 call rone(this%s, this%lgmres)
213 call rone(this%c, this%lgmres)
214 call rzero(this%h, this%lgmres * this%lgmres)
216 do while (.not. conv .and. iter .lt. max_iter)
220 call col3(this%r,this%ml,f,n)
223 call copy (this%r,f,n)
224 call ax%compute(this%w, x%x, coef, x%msh, x%Xh)
225 call gs_h%op(this%w, n, gs_op_add)
227 call add2s2(this%r,this%w,-one,n)
228 call col2(this%r,this%ml,n)
230 this%gam(1) = sqrt(
glsc3(this%r, this%r, coef%mult, n))
232 div0 = this%gam(1) * norm_fac
233 ksp_results%res_start = div0
236 if (
abscmp(this%gam(1), 0.0_rp))
return
239 temp = one / this%gam(1)
240 call cmult2(this%v(1,1), this%r, temp, n)
241 do j = 1, this%lgmres
243 call col3(this%w, this%mu, this%v(1,j), n)
246 call this%M%solve(this%z(1,j), this%w, n)
248 call ax%compute(this%w, this%z(1,j), coef, x%msh, x%Xh)
249 call gs_h%op(this%w, n, gs_op_add)
251 call col2(this%w, this%ml, n)
256 this%h(i,j) = this%h(i,j) + &
257 this%w(k) * this%v(k,i) * coef%mult(k,1,1,1)
262 call mpi_allreduce(this%h(1,j), this%wk1, j, &
264 call copy(this%h(1,j), this%wk1, j)
268 this%w(k) = this%w(k) - this%h(i,j) * this%v(k,i)
275 this%h(i ,j) = this%c(i)*temp + this%s(i)*this%h(i+1,j)
276 this%h(i+1,j) = -this%s(i)*temp + this%c(i)*this%h(i+1,j)
279 alpha = sqrt(
glsc3(this%w, this%w, coef%mult, n))
281 if(
abscmp(alpha, 0.0_rp))
then
285 l = sqrt(this%h(j,j) * this%h(j,j) + alpha**2)
287 this%c(j) = this%h(j,j) * temp
288 this%s(j) = alpha * temp
290 this%gam(j+1) = -this%s(j) * this%gam(j)
291 this%gam(j) = this%c(j) * this%gam(j)
293 rnorm = abs(this%gam(j+1)) * norm_fac
295 if (rnorm .lt. this%abs_tol)
then
300 if (iter + 1 .gt. max_iter)
exit
302 if( j .lt. this%lgmres)
then
304 call cmult2(this%v(1,j+1), this%w, temp, n)
307 j = min(j, this%lgmres)
312 temp = temp - this%h(k,i) * this%c(i)
314 this%c(k) = temp / this%h(k,k)
319 x%x(k,1,1,1) = x%x(k,1,1,1) + this%c(i) * this%z(k,i)
324 ksp_results%res_final = rnorm
325 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 sx_gmres_free(this)
Deallocate a standard GMRES solver.
subroutine sx_gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol)
Initialise a standard GMRES solver.
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.