139 integer,
intent(in) :: n
140 integer,
intent(in) :: max_iter
141 class(pc_t),
optional,
intent(in),
target :: M
142 real(kind=rp),
optional,
intent(in) :: rel_tol
143 real(kind=rp),
optional,
intent(in) :: abs_tol
144 logical,
optional,
intent(in) :: monitor
145 type(device_ident_t),
target :: M_ident
147 integer(c_size_t) :: z_size
160 call device_map(this%w, this%w_d, n)
161 call device_map(this%r, this%r_d, n)
163 allocate(this%c(this%m_restart))
164 allocate(this%s(this%m_restart))
165 allocate(this%gam(this%m_restart + 1))
166 call device_map(this%c, this%c_d, this%m_restart)
167 call device_map(this%s, this%s_d, this%m_restart)
168 call device_map(this%gam, this%gam_d, this%m_restart+1)
170 allocate(this%z(n, this%m_restart))
171 allocate(this%v(n, this%m_restart))
172 allocate(this%h(this%m_restart, this%m_restart))
173 allocate(this%z_d(this%m_restart))
174 allocate(this%v_d(this%m_restart))
175 allocate(this%h_d(this%m_restart))
176 do i = 1, this%m_restart
177 this%z_d(i) = c_null_ptr
178 call device_map(this%z(:,i), this%z_d(i), n)
180 this%v_d(i) = c_null_ptr
181 call device_map(this%v(:,i), this%v_d(i), n)
183 this%h_d(i) = c_null_ptr
184 call device_map(this%h(:,i), this%h_d(i), this%m_restart)
187 z_size = c_sizeof(c_null_ptr) * (this%m_restart)
188 call device_alloc(this%z_d_d, z_size)
189 call device_alloc(this%v_d_d, z_size)
190 call device_alloc(this%h_d_d, z_size)
191 ptr = c_loc(this%z_d)
192 call device_memcpy(ptr, this%z_d_d, z_size, &
193 host_to_device, sync = .false.)
194 ptr = c_loc(this%v_d)
195 call device_memcpy(ptr, this%v_d_d, z_size, &
196 host_to_device, sync = .false.)
197 ptr = c_loc(this%h_d)
198 call device_memcpy(ptr, this%h_d_d, z_size, &
199 host_to_device, sync = .false.)
202 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
203 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
204 else if (
present(rel_tol) .and.
present(abs_tol))
then
205 call this%ksp_init(max_iter, rel_tol, abs_tol)
206 else if (
present(monitor) .and.
present(abs_tol))
then
207 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
208 else if (
present(rel_tol) .and.
present(monitor))
then
209 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
210 else if (
present(rel_tol))
then
211 call this%ksp_init(max_iter, rel_tol = rel_tol)
212 else if (
present(abs_tol))
then
213 call this%ksp_init(max_iter, abs_tol = abs_tol)
214 else if (
present(monitor))
then
215 call this%ksp_init(max_iter, monitor = monitor)
217 call this%ksp_init(max_iter)
220 call device_event_create(this%gs_event, 2)
315 class(ax_t),
intent(in) :: ax
316 type(field_t),
intent(inout) :: x
317 integer,
intent(in) :: n
318 real(kind=rp),
dimension(n),
intent(in) :: f
319 type(coef_t),
intent(inout) :: coef
320 type(bc_list_t),
intent(inout) :: blst
321 type(gs_t),
intent(inout) :: gs_h
322 type(ksp_monitor_t) :: ksp_results
323 integer,
optional,
intent(in) :: niter
324 integer :: iter, max_iter
326 real(kind=rp) :: rnorm, alpha, temp, lr, alpha2, norm_fac
330 f_d = device_get_ptr(f)
336 if (
present(niter))
then
339 max_iter = this%max_iter
342 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
343 v => this%v, s => this%s, gam => this%gam, v_d => this%v_d, &
344 w_d => this%w_d, r_d => this%r_d, h_d => this%h_d, &
345 v_d_d => this%v_d_d, x_d => x%x_d, z_d_d => this%z_d_d, &
348 norm_fac = 1.0_rp / sqrt(coef%volume)
349 call rzero(gam, this%m_restart + 1)
350 call rone(s, this%m_restart)
351 call rone(c, this%m_restart)
352 call rzero(h, this%m_restart * this%m_restart)
353 call device_rzero(x%x_d, n)
354 call device_rzero(this%gam_d, this%m_restart + 1)
355 call device_rone(this%s_d, this%m_restart)
356 call device_rone(this%c_d, this%m_restart)
358 call rzero(this%h, this%m_restart**2)
363 call this%monitor_start(
'GMRES')
364 do while (.not. conv .and. iter .lt. max_iter)
366 if (iter .eq. 0)
then
367 call device_copy(r_d, f_d, n)
369 call device_copy(r_d, f_d, n)
370 call ax%compute(w, x%x, coef, x%msh, x%Xh)
371 call gs_h%op(w, n, gs_op_add, this%gs_event)
372 call device_event_sync(this%gs_event)
373 call blst%apply_scalar(w, n)
374 call device_sub2(r_d, w_d, n)
377 gam(1) = sqrt(device_glsc3(r_d, r_d, coef%mult_d, n))
378 if (iter .eq. 0)
then
379 ksp_results%res_start = gam(1) * norm_fac
382 if (abscmp(gam(1), 0.0_rp))
exit
385 temp = 1.0_rp / gam(1)
386 call device_cmult2(v_d(1), r_d, temp, n)
387 do j = 1, this%m_restart
390 call this%M%solve(z(1,j), v(1,j), n)
392 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
393 call gs_h%op(w, n, gs_op_add, this%gs_event)
394 call device_event_sync(this%gs_event)
395 call blst%apply_scalar(w, n)
397 if (neko_bcknd_opencl .eq. 1)
then
399 h(i,j) = device_glsc3(w_d, v_d(i), coef%mult_d, n)
401 call device_add2s2(w_d, v_d(i), -h(i,j), n)
403 alpha2 = device_glsc3(w_d, w_d, coef%mult_d, n)
406 call device_glsc3_many(h(1,j), w_d, v_d_d, coef%mult_d, j, n)
408 call device_memcpy(h(:,j), h_d(j), j, &
409 host_to_device, sync = .false.)
419 h(i,j) = c(i)*temp + s(i) * h(i+1,j)
420 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
424 if (abscmp(alpha, 0.0_rp))
then
429 lr = sqrt(h(j,j) * h(j,j) + alpha2)
434 call device_memcpy(h(:,j), h_d(j), j, &
435 host_to_device, sync = .false.)
436 gam(j+1) = -s(j) * gam(j)
437 gam(j) = c(j) * gam(j)
439 rnorm = abs(gam(j+1)) * norm_fac
440 call this%monitor_iter(iter, rnorm)
441 if (rnorm .lt. this%abs_tol)
then
446 if (iter + 1 .gt. max_iter)
exit
448 if (j .lt. this%m_restart)
then
449 temp = 1.0_rp / alpha
450 call device_cmult2(v_d(j+1), w_d, temp, n)
455 j = min(j, this%m_restart)
459 temp = temp - h(k,i) * c(i)
464 if (neko_bcknd_opencl .eq. 1)
then
466 call device_add2s2(x_d, this%z_d(i), c(i), n)
469 call device_memcpy(c, c_d, j, host_to_device, sync = .false.)
470 call device_add2s2_many(x_d, z_d_d, c_d, j, n)
475 call this%monitor_stop()
476 ksp_results%res_final = rnorm
477 ksp_results%iter = iter
478 ksp_results%converged = this%is_converged(iter, rnorm)
484 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
486 class(ax_t),
intent(in) :: ax
487 type(field_t),
intent(inout) :: x
488 type(field_t),
intent(inout) :: y
489 type(field_t),
intent(inout) :: z
490 integer,
intent(in) :: n
491 real(kind=rp),
dimension(n),
intent(in) :: fx
492 real(kind=rp),
dimension(n),
intent(in) :: fy
493 real(kind=rp),
dimension(n),
intent(in) :: fz
494 type(coef_t),
intent(inout) :: coef
495 type(bc_list_t),
intent(inout) :: blstx
496 type(bc_list_t),
intent(inout) :: blsty
497 type(bc_list_t),
intent(inout) :: blstz
498 type(gs_t),
intent(inout) :: gs_h
499 type(ksp_monitor_t),
dimension(3) :: ksp_results
500 integer,
optional,
intent(in) :: niter
502 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
503 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
504 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)