72 rel_tol, abs_tol, monitor)
73 class(
gmres_t),
intent(inout) :: this
74 integer,
intent(in) :: n
75 integer,
intent(in) :: max_iter
76 class(
pc_t),
optional,
intent(in),
target :: M
77 integer,
optional,
intent(in) :: lgmres
78 real(kind=
rp),
optional,
intent(in) :: rel_tol
79 real(kind=
rp),
optional,
intent(in) :: abs_tol
80 logical,
optional,
intent(in) :: monitor
82 if (
present(lgmres))
then
98 allocate(this%c(this%lgmres))
99 allocate(this%s(this%lgmres))
100 allocate(this%gam(this%lgmres + 1))
102 allocate(this%z(n, this%lgmres))
103 allocate(this%v(n, this%lgmres))
105 allocate(this%h(this%lgmres, this%lgmres))
107 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
108 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
109 else if (
present(rel_tol) .and.
present(abs_tol))
then
110 call this%ksp_init(max_iter, rel_tol, abs_tol)
111 else if (
present(monitor) .and.
present(abs_tol))
then
112 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
113 else if (
present(rel_tol) .and.
present(monitor))
then
114 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
115 else if (
present(rel_tol))
then
116 call this%ksp_init(max_iter, rel_tol = rel_tol)
117 else if (
present(abs_tol))
then
118 call this%ksp_init(max_iter, abs_tol = abs_tol)
119 else if (
present(monitor))
then
120 call this%ksp_init(max_iter, monitor = monitor)
122 call this%ksp_init(max_iter)
171 function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
173 class(
gmres_t),
intent(inout) :: this
174 class(
ax_t),
intent(in) :: ax
175 type(
field_t),
intent(inout) :: x
176 integer,
intent(in) :: n
177 real(kind=
rp),
dimension(n),
intent(in) :: f
178 type(
coef_t),
intent(inout) :: coef
180 type(
gs_t),
intent(inout) :: gs_h
182 integer,
optional,
intent(in) :: niter
183 integer :: iter, max_iter
184 integer :: i, j, k, l, ierr
186 real(kind=
xp) :: alpha, lr, alpha2, norm_fac
187 real(kind=
rp) :: temp, rnorm
193 if (
present(niter))
then
196 max_iter = this%max_iter
199 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
200 v => this%v, s => this%s, gam => this%gam)
202 norm_fac = 1.0_rp / sqrt(coef%volume)
208 call this%monitor_start(
'GMRES')
209 do while (.not. conv .and. iter .lt. max_iter)
211 if (iter .eq. 0)
then
215 call ax%compute(w, x%x, coef, x%msh, x%Xh)
216 call gs_h%op(w, n, gs_op_add)
217 call blst%apply(w, n)
221 gam(1) = sqrt(
glsc3(r, r, coef%mult, n))
222 if (iter .eq. 0)
then
223 ksp_results%res_start = gam(1) * norm_fac
226 if (
abscmp(gam(1), 0.0_xp))
return
229 temp = 1.0_rp / gam(1)
230 call cmult2(v(1,1), r, temp, n)
231 do j = 1, this%lgmres
234 call this%M%solve(z(1,j), v(1,j), n)
236 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
237 call gs_h%op(w, n, gs_op_add)
238 call blst%apply(w, n)
249 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
256 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
262 call mpi_allreduce(mpi_in_place, h(1,j), j, &
273 w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
277 w(i+k) = w(i+k) + w_plus(k)
278 alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
284 w_plus(1) = w_plus(1) - h(l,j) * v(i+k,l)
286 w(i+k) = w(i+k) + w_plus(1)
287 alpha2 = alpha2 + (w(i+k)**2) * coef%mult(i+k,1,1,1)
292 call mpi_allreduce(mpi_in_place,alpha2, 1, &
297 h(i,j) = c(i)*temp + s(i) * h(i+1,j)
298 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
302 if (
abscmp(alpha, 0.0_xp))
then
307 lr = sqrt(h(j,j) * h(j,j) + alpha2)
312 gam(j+1) = -s(j) * gam(j)
313 gam(j) = c(j) * gam(j)
314 rnorm = abs(gam(j+1)) * norm_fac
315 call this%monitor_iter(iter, rnorm)
316 if (rnorm .lt. this%abs_tol)
then
321 if (iter + 1 .gt. max_iter)
exit
323 if (j .lt. this%lgmres)
then
324 temp = 1.0_rp / alpha
325 call cmult2(v(1,j+1), w, temp, n)
330 j = min(j, this%lgmres)
334 temp = temp - h(k,i) * c(i)
346 x_plus(k) = x_plus(k) + c(l) * z(i+k,l)
350 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
356 x_plus(1) = x_plus(1) + c(l) * z(i+k,l)
358 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
365 call this%monitor_stop()
366 ksp_results%res_final = rnorm
367 ksp_results%iter = iter
368 ksp_results%converged = this%is_converged(iter, rnorm)
374 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
375 class(
gmres_t),
intent(inout) :: this
376 class(ax_t),
intent(in) :: ax
377 type(field_t),
intent(inout) :: x
378 type(field_t),
intent(inout) :: y
379 type(field_t),
intent(inout) :: z
380 integer,
intent(in) :: n
381 real(kind=rp),
dimension(n),
intent(in) :: fx
382 real(kind=rp),
dimension(n),
intent(in) :: fy
383 real(kind=rp),
dimension(n),
intent(in) :: fz
384 type(coef_t),
intent(inout) :: coef
385 type(bc_list_t),
intent(inout) :: blstx
386 type(bc_list_t),
intent(inout) :: blsty
387 type(bc_list_t),
intent(inout) :: blstz
388 type(gs_t),
intent(inout) :: gs_h
389 type(ksp_monitor_t),
dimension(3) :: ksp_results
390 integer,
optional,
intent(in) :: niter
392 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
393 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
394 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
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.