Line data Source code
1 : !!****m* ABINIT/m_iterative_solvers
2 : !! NAME
3 : !! m_iterative_solvers
4 : !!
5 : !! FUNCTION
6 : !! This module contains basic (matrix-free) MPI-aware iterative solvers (GMRES, CG).
7 : !!
8 : !! COPYRIGHT
9 : !! TODO
10 : !!
11 : !! SOURCE
12 :
13 : #if defined HAVE_CONFIG_H
14 : #include "config.h"
15 : #endif
16 :
17 : #include "abi_common.h"
18 :
19 : module m_iterative_solvers
20 :
21 : use m_errors
22 : use defs_basis
23 : use m_xmpi
24 : use m_abicore
25 : use m_specialmsg
26 :
27 : implicit none
28 : private
29 : public :: cg_linear_solver, gmres_linear_solver, cg_eigen_solver_treshold
30 :
31 : contains
32 :
33 : !-------------------------------------------------------------------------------------------
34 : ! Eigensolvers
35 : !-------------------------------------------------------------------------------------------
36 :
37 : !****f* m_iterative_solvers/cg_eigen_solver_treshold
38 : !! NAME
39 : !! cg_eigen_solver_treshold
40 : !!
41 : !! FUNCTION
42 : !! Compute the smallest eigenvalues and corresponding eigenvectors of a matrix using
43 : !! the Conjugate Gradient method to minimize the Rayleigh quotient.
44 : !! The stopping criterion is : at least one of the computed eigenvalues is above 'treshold'.
45 : !!
46 : !! INPUTS
47 : !! n = Size of the matrix.
48 : !! matvec = Subroutine that performs matrix-vector multiplication.
49 : !! x0 = Initial guess for the eigenvector.
50 : !! tol = Convergence tolerance for the residual norm.
51 : !! max_iter = Maximum number of iterations for the Conjugate Gradient method.
52 : !! max_neig = Maximum number of eigenvalues to compute.
53 : !! eigenvalue_threshold = Threshold below which to stop computing further eigenvalues.
54 : !!
55 : !! OUTPUTS
56 : !! eigenvalues = Array of computed smallest eigenvalues.
57 : !! eigenvectors = Matrix of corresponding eigenvectors.
58 : !! n_eig = Number of computed eigenvalues.
59 : !!
60 : !! SOURCE
61 0 : subroutine cg_eigen_solver_treshold(n, matvec, x0, tol, max_iter, max_neig, eigenvalue_threshold, eigenvalues, eigenvectors, n_eig)
62 :
63 : ! Input parameters
64 : integer, intent(in) :: n, max_iter, max_neig
65 : real(dp), intent(in) :: tol, eigenvalue_threshold
66 : real(dp), intent(in) :: x0(n)
67 : interface
68 : subroutine matvec(n_, x, y)
69 : integer, intent(in) :: n_
70 : double precision, intent(inout), target :: x(n_), y(n_)
71 : end subroutine matvec
72 : end interface
73 :
74 : ! Output parameters
75 : real(dp), intent(out) :: eigenvalues(max_neig)
76 : real(dp), intent(out) :: eigenvectors(n, max_neig)
77 : integer, intent(out) :: n_eig
78 :
79 : ! Local variables
80 0 : real(dp) :: x(n), Ax(n), p(n), g(n), beta
81 : real(dp) :: residual_norm, rayleigh_quotient
82 : integer :: i, j, iter
83 : logical :: do_exit
84 : integer :: ierr
85 :
86 : ! *************************************************************************
87 :
88 : ! Initialize variables
89 0 : x = x0
90 0 : call matvec(n, x, Ax) ! Compute initial matrix-vector product
91 0 : rayleigh_quotient = dot_product(x, Ax) / dot_product(x, x)
92 0 : g = 2/dot_product(x, x) * (Ax - rayleigh_quotient * x) ! Gradient of Rayleigh quotient
93 0 : p = 0
94 0 : beta = 0
95 :
96 : ! Conjugate Gradient iterations to minimize Rayleigh quotient
97 0 : do iter = 1, max_iter
98 0 : call cg_update(n, matvec, x, Ax, p, g, beta, rayleigh_quotient, residual_norm)
99 0 : do_exit = (residual_norm < tol)
100 0 : call xmpi_bcast(do_exit, 0, xmpi_world, ierr)
101 0 : if (do_exit) exit
102 : end do
103 :
104 : ! Store the computed eigenvalue and eigenvector
105 0 : eigenvalues(1) = rayleigh_quotient
106 0 : eigenvectors(:, 1) = x/norm2(x)
107 0 : n_eig = 1
108 :
109 : ! If more eigenvalues are required, use deflation to compute subsequent eigenvalues
110 0 : do j = 2, max_neig
111 :
112 : ! Check if the eigenvalue is below the threshold
113 0 : do_exit = (eigenvalues(j-1) < eigenvalue_threshold)
114 0 : call xmpi_bcast(do_exit, 0, xmpi_world, ierr)
115 0 : if (do_exit) exit
116 :
117 : ! If it is not, we need to compute more eigenvalues :
118 : ! Orthogonalize the initial guess against previously computed eigenvectors
119 0 : x = x0
120 0 : do i = 1, j - 1
121 0 : x = x - dot_product(p, eigenvectors(:, i))/dot_product(eigenvectors(:, i), eigenvectors(:, i)) * eigenvectors(:, i)
122 : end do
123 :
124 : ! Initialization
125 0 : call matvec(n, x, Ax) ! Compute initial matrix-vector product
126 0 : rayleigh_quotient = dot_product(x, Ax) / dot_product(x, x)
127 0 : g = 2/dot_product(x, x) * (Ax - rayleigh_quotient * x) ! Gradient of Rayleigh quotient
128 0 : p = 0
129 0 : beta = 0
130 :
131 : ! Perform conjugate gradient iterations
132 0 : do iter = 1, max_iter
133 0 : call cg_update(n, matvec, x, Ax, p, g, beta, rayleigh_quotient, residual_norm, eigenvectors(:, 1:j-1))
134 0 : do_exit = (residual_norm < tol)
135 0 : call xmpi_bcast(do_exit, 0, xmpi_world, ierr) ! MPI aware: avoid desynchronization.
136 0 : if (do_exit) exit
137 : end do
138 :
139 : ! Store the computed eigenvalue and eigenvector
140 0 : eigenvalues(j) = rayleigh_quotient
141 0 : eigenvectors(:, j) = x/norm2(x)
142 0 : n_eig = n_eig + 1
143 :
144 : end do
145 :
146 0 : end subroutine cg_eigen_solver_treshold
147 :
148 : !****f* m_iterative_solvers/cg_update
149 : !! NAME
150 : !! cg_update
151 : !!
152 : !! FUNCTION
153 : !! Perform a single Conjugate Gradient update step to minimize the Rayleigh quotient
154 : !! of a self-adjoint operator given by matvec. This method is used to approximate
155 : !! the eigenvalue and eigenvector of the operator.
156 : !!
157 : !! INPUTS
158 : !! n = Integer, size of the matrix (number of rows/columns).
159 : !! matvec = Subroutine, performs matrix-vector multiplication (A * v).
160 : !!
161 : !! INPUT/OUTPUTS
162 : !! x = Current solution vector (eigenvector approximation).
163 : !! Ax = A * x, where A is the operator defined by matvec.
164 : !! p = Previous search direction vector.
165 : !! g = Gradient of the Rayleigh-quotient at x.
166 : !! beta = Scalar, used to update the search direction.
167 : !!
168 : !! OUTPUTS
169 : !! rayleigh_quotient = Current approximation of the eigenvalue (lambda).
170 : !! residual_norm = Norm of the residual.
171 : !!
172 : !! NOTES
173 : !! - The input vectors (x, r, p, Ap) must be properly initialized before calling
174 : !! this subroutine.
175 : !! - The subroutine assumes that the operator is self-adjoint (Hermitian).
176 : !!
177 : !! SOURCE
178 0 : subroutine cg_update(n, matvec, x, Ax, p, g, beta, rayleigh_quotient, residual_norm, eigenvectors)
179 : ! Arguments
180 : integer, intent(in) :: n
181 : real(dp), optional, intent(in) :: eigenvectors(:,:)
182 : interface
183 : subroutine matvec(n_, x, y)
184 : integer, intent(in) :: n_
185 : double precision, intent(inout), target :: x(n_), y(n_)
186 : end subroutine matvec
187 : end interface
188 : real(dp), intent(inout) :: x(n), Ax(n), g(n), p(n), beta
189 : real(dp), intent(inout) :: rayleigh_quotient, residual_norm
190 :
191 : ! Local variables
192 0 : real(dp) :: Ap(n), a, b, c, d, e, f, alpha, alpha_(2)
193 : integer :: n_eig, i
194 :
195 : ! *************************************************************************
196 0 : if (present(eigenvectors)) then
197 0 : n_eig = size(eigenvectors, 2) ! Number of eigenvectors already computed
198 : else
199 : n_eig = 0
200 : end if
201 :
202 : ! Update search direction
203 0 : p = -g + beta*p
204 : ! Orthogonalize the search direction against previously computed eigenvectors
205 0 : do i = 1, n_eig
206 0 : p = p - dot_product(p, eigenvectors(:, i))/dot_product(eigenvectors(:, i), eigenvectors(:, i)) * eigenvectors(:, i)
207 : end do
208 :
209 : ! Compute A*p
210 0 : call matvec(n, p, Ap)
211 :
212 : ! Compute alpha (step size = minimizer of R(x+alpha*p) that is the solution (+) of a quadratic problem)
213 0 : a = dot_product(p, Ap)
214 0 : b = 2*dot_product(x, Ap)
215 0 : c = dot_product(x, Ax)
216 0 : d = dot_product(p, p)
217 0 : e = 2*dot_product(x, p)
218 0 : f = dot_product(x, x)
219 0 : alpha_ = quadratic_roots(a*e-b*d, 2*(f*a-d*c), b*f-c*e)
220 0 : alpha = alpha_(1)
221 :
222 : ! Update solution vector
223 0 : x = x + alpha * p
224 0 : Ax = Ax + alpha * Ap
225 :
226 : ! Update the Rayleigh-quotient
227 0 : rayleigh_quotient = dot_product(x, Ax)/dot_product(x, x)
228 :
229 : ! Update the (x-normalized) gradient and beta = dot(g, g)/dot(g_old, g_old)
230 0 : beta = 1/dot_product(g, g)
231 0 : g = 2/norm2(x) * (Ax - rayleigh_quotient * x)
232 0 : beta = beta * dot_product(g, g)
233 :
234 : ! Update the residual norm (= norm of the gradient)
235 0 : residual_norm = norm2(g)
236 :
237 0 : end subroutine cg_update
238 :
239 0 : function quadratic_roots(a, b, c) result(r)
240 : ! Arguments
241 : real(dp) :: a, b, c
242 : real(dp) :: r(2)
243 :
244 : ! Local variables
245 : real(dp) :: discriminant, sqrt_discriminant
246 :
247 : ! *************************************************************************
248 :
249 0 : if (a == 0.0_dp) then
250 0 : r(1) = 0.0_dp
251 0 : r(2) = 0.0_dp
252 0 : return
253 : end if
254 :
255 0 : discriminant = b**2 - 4.0_dp * a * c
256 0 : if (discriminant < 0.0_dp) then
257 0 : r(1) = 0.0_dp
258 0 : r(2) = 0.0_dp
259 0 : return
260 : end if
261 :
262 0 : sqrt_discriminant = sqrt(discriminant)
263 0 : r(1) = (-b + sqrt_discriminant) / (2.0_dp * a)
264 0 : r(2) = (-b - sqrt_discriminant) / (2.0_dp * a)
265 :
266 : end function quadratic_roots
267 :
268 : !-------------------------------------------------------------------------------------------
269 : ! Linear solvers (CG and GMRES)
270 : !-------------------------------------------------------------------------------------------
271 :
272 0 : subroutine cg_linear_solver(n, matvec, rhs, est, cg_maxiter, cg_rtol, verbose)
273 :
274 : !Arguments ------------------------------------
275 : integer, intent(in) :: n, cg_maxiter
276 : real(dp), intent(in) :: cg_rtol
277 : real(dp), intent(in) :: rhs(:)
278 : logical, intent(in) :: verbose
279 : real(dp),intent(inout) :: est(:)
280 : character(len=500) :: msg
281 : interface
282 : subroutine matvec(n_, x, y)
283 : integer, intent(in) :: n_
284 : double precision, intent(inout), target :: x(n_), y(n_)
285 : end subroutine matvec
286 : end interface
287 :
288 : !Local variables-------------------------------
289 0 : real(dp) :: r(n), p(n), Ap(n)
290 : real(dp) :: alpha, beta, rsold, rsnew
291 : integer :: iter
292 : logical :: do_exit
293 : integer :: ierr
294 :
295 : ! *************************************************************************
296 :
297 : ! Initialize
298 0 : call matvec(n, est, Ap)
299 0 : r = rhs - Ap
300 0 : p = r
301 0 : rsold = dot_product(r, r)
302 :
303 : ! Conjugate Gradient iterations
304 0 : do iter = 1, cg_maxiter
305 0 : call matvec(n, p, Ap)
306 0 : alpha = rsold / dot_product(p, Ap)
307 0 : est = est + alpha * p
308 0 : r = r - alpha * Ap
309 0 : rsnew = dot_product(r, r)
310 :
311 : ! Check for convergence
312 0 : do_exit = (sqrt(rsnew) < cg_rtol)
313 0 : call xmpi_bcast(do_exit, 0, xmpi_world, ierr) ! MPI aware: avoid desynchronization.
314 0 : if (do_exit) exit
315 :
316 0 : beta = rsnew / rsold
317 0 : p = r + beta * p
318 0 : write(msg, *)'cg: it=', iter,' res=', sqrt(rsnew)
319 0 : if (verbose) call wrtout(std_out, msg)
320 : !if (verbose) write(std_out,*) 'cg: it=', iter,' res=', sqrt(rsnew)
321 0 : rsold = rsnew
322 : end do
323 :
324 0 : end subroutine cg_linear_solver
325 :
326 : !****f* m_iterative_solvers/call_gmresm
327 : !! NAME
328 : !! call_gmresm
329 : !!
330 : !! FUNCTION
331 : !! Call the gmresm (Willis, A. (2017) SoftwareX 6, 124-127, code at the end of this file) routine to solve a linear system.
332 : !!
333 : !! INPUTS
334 : !! n = Size of the matrix.
335 : !! matvec = Subroutine that performs matrix-vector multiplication.
336 : !! rhs = Right-hand side vector of the linear system.
337 : !! gmres_maxiter = Maximum number of iterations for the GMRES algorithm.
338 : !! gmres_rtol = Relative tolerance for convergence.
339 : !! verbose = Logical, if true, print residuals at each iteration.
340 : !!
341 : !! INPUT/OUTPUTS
342 : !! est = Initial guess for the solution vector, updated with the computed solution.
343 : !!
344 : !! SOURCE
345 18 : subroutine call_gmresm(n, matvec, est, rhs, gmres_maxiter, gmres_rtol, verbose)
346 : !Arguments ------------------------------------
347 : integer, intent(in) :: n, gmres_maxiter
348 : real(dp), intent(in) :: gmres_rtol
349 : real(dp), intent(in) :: rhs(n)
350 : logical, intent(in) :: verbose
351 : real(dp), intent(inout) :: est(n)
352 : interface
353 : subroutine matvec(n_, x, y)
354 : integer, intent(in) :: n_
355 : double precision, intent(inout), target :: x(n_), y(n_)
356 : end subroutine matvec
357 : end interface
358 : !Local variables-------------------------------
359 : integer :: its, info, m
360 : real(dp) :: res, del
361 18 : real(dp), allocatable :: h(:, :), v(:, :)
362 :
363 : ! *************************************************************************
364 :
365 18 : m = gmres_maxiter
366 72 : ABI_MALLOC(h, (m+1, m))
367 72 : ABI_MALLOC(v, (n, m+1))
368 121518 : res = gmres_rtol * norm2(rhs)
369 18 : del = 0
370 18 : its = gmres_maxiter ! No restart
371 18 : info = 0
372 18 : if (verbose) then
373 0 : info = 1
374 : end if
375 18 : call gmresm(m, n, est, rhs, matvec, psolve, dotprd, h, v, res, del, its, info)
376 18 : ABI_FREE(h)
377 18 : ABI_FREE(v)
378 :
379 : contains
380 : ! Dummy : No preconditioning
381 108 : subroutine psolve(n_, x)
382 : integer, intent(in) :: n_
383 : real(dp), intent(inout) :: x(n_)
384 : ! ***********************
385 : ! We do nothing here but don't wan't to be flashed by abirule.
386 : if (.false.) then
387 : x = zero
388 : end if
389 108 : end subroutine psolve
390 : ! Dot product
391 406 : function dotprd(n_, a, b) result(c)
392 : integer, intent(in) :: n_
393 : real(dp), intent(inout) :: a(n_), b(n_)
394 : real(dp) :: c
395 : ! ***********************
396 2740906 : c = dot_product(a, b)
397 406 : end function dotprd
398 :
399 : end subroutine call_gmresm
400 :
401 : !****f* m_iterative_solvers/gmres_linear_solver
402 : !! NAME
403 : !! gmres_linear_solver
404 : !!
405 : !! FUNCTION
406 : !! Solve a linear system using GMRES. Depending on the availability of MKL,
407 : !! it either calls the MKL FGMRES routine or the gmresm routine.
408 : !!
409 : !! INPUTS
410 : !! n = Size of the matrix.
411 : !! matvec = Subroutine that performs matrix-vector multiplication.
412 : !! rhs = Right-hand side vector of the linear system.
413 : !! gmres_maxiter = Maximum number of iterations for the GMRES algorithm.
414 : !! gmres_rtol = Relative tolerance for convergence.
415 : !!
416 : !! INPUT/OUTPUTS
417 : !! est = Initial guess for the solution vector, updated with the computed solution.
418 : !!
419 : !! SOURCE
420 18 : subroutine gmres_linear_solver(n, matvec, rhs, est, gmres_maxiter, gmres_rtol, verbose)
421 : !Arguments ------------------------------------
422 : integer, intent(in) :: n, gmres_maxiter
423 : real(dp), intent(in) :: gmres_rtol
424 : real(dp), intent(in) :: rhs(n)
425 : logical :: verbose
426 : real(dp), intent(inout) :: est(n)
427 : interface
428 : subroutine matvec(n_, x, y)
429 : integer, intent(in) :: n_
430 : double precision, intent(inout), target :: x(n_), y(n_)
431 : end subroutine matvec
432 : end interface
433 :
434 : ! *************************************************************************
435 :
436 18 : call call_gmresm(n, matvec, est, rhs, gmres_maxiter, gmres_rtol, verbose)
437 :
438 18 : end subroutine gmres_linear_solver
439 :
440 : !-------------------------------------------------------------------------------------------
441 :
442 : !----------------------------------------------------------------------
443 : ! Openpipeflow.org. If used in your work, please cite
444 : ! Willis, A. (2017) SoftwareX 6, 124-127.
445 : ! https://doi.org/10.1016/j.softx.2017.05.003 (open access)
446 : ! Thanks in advance! Ashley 2019.
447 : !----------------------------------------------------------------------
448 : ! solve A x = b for x ;
449 : ! minimise |Ax-b| subject to constraint |x| < delta .
450 : ! requires lapack routines dgelsy, dgesvd.
451 : !----------------------------------------------------------------------
452 : ! m gmres dimension
453 : ! n dimension of x
454 : ! x on input: guess for x, can be 0
455 : ! on exit: solution x, subject to constraint if del>0
456 : ! b input b
457 : ! matvec performs y := A x, call matvec(N,x, y)
458 : ! psolve preconditioner, solve M x_out = x_in, call psolve(N,x)
459 : ! dotprd dot product, d = dotprd(n,a,b)
460 : ! h Hessian matrix, size (m+1)*m
461 : ! v Krylov subspace, size n*(m+1)
462 : ! res on input: |Ax-b|/|b|<res;
463 : ! on exit: residual reached
464 : ! del on input: if(del>0) then the x returned is the hookstep
465 : ! on exit: norm of next b predicted by hook
466 : ! its on input: max num its;
467 : ! on exit: number of its taken
468 : ! info on input: if(info==1) print* residuals
469 : ! if(info==2) recalc hookstep with new del
470 : ! on exit: 0 sucessful, 1 method breakdown, 2 max its
471 : ! A.P.Willis 2008
472 : !----------------------------------------------------------------------
473 :
474 18 : subroutine gmresm(m,n,x,b,matvec,psolve,dotprd,h,v,res,del,its,info)
475 : implicit none
476 : integer, intent(in) :: m
477 : integer, intent(in) :: n
478 : real(dp), intent(inout), target :: x(n)
479 : real(dp), intent(in) :: b(n)
480 : interface
481 : subroutine matvec(n_, x, y)
482 : integer, intent(in) :: n_
483 : double precision, intent(inout), target :: x(n_), y(n_)
484 : end subroutine matvec
485 : subroutine psolve(n_, x)
486 : integer, intent(in) :: n_
487 : double precision, intent(inout) :: x(n_)
488 : end subroutine psolve
489 : function dotprd(n_, a, b) result(c)
490 : integer, intent(in) :: n_
491 : double precision, intent(inout) :: a(n_), b(n_)
492 : double precision :: c
493 : end function dotprd
494 : end interface
495 : real(dp), intent(inout) :: h(m+1,m)
496 : real(dp), intent(inout) :: v(n,m+1)
497 : real(dp), intent(inout) :: res
498 : real(dp), intent(inout) :: del
499 : integer, intent(inout) :: its
500 : integer, intent(inout) :: info
501 : real(dp) :: tol,res_,stgn
502 36 : real(dp), target :: w(n), z(n)
503 36 : real(dp) :: h_(m+1,m), y(m+1), p(m+1), work(4*m+1)
504 18 : integer :: imx, piv(m), rank, i
505 : real(dp), save :: beta
506 : integer, save :: j
507 : logical :: done
508 : integer :: ierr
509 : character(len=500) :: msg
510 :
511 18 : if(info==2) then
512 0 : call hookstep(j,h,m,beta,del, y)
513 0 : z = matmul(v(:,1:j),y(1:j))
514 0 : call psolve(n, z)
515 0 : x = z
516 0 : info = 0
517 0 : return
518 : end if
519 :
520 18 : tol = res
521 18 : imx = its
522 18 : its = 0
523 2551896 : v = 0d0
524 :
525 : 1 continue
526 18 : res_ = 1d99
527 18 : stgn = 1d0 - 1d-14
528 :
529 18 : beta = dsqrt(dotprd(n,x,x))
530 121518 : if(beta==0d0) w = 0d0
531 18 : if(beta/=0d0) call matvec(n,x, w)
532 121518 : w = b - w
533 18 : beta = dsqrt(dotprd(n,w,w))
534 121518 : v(:,1) = w / beta
535 :
536 7938 : h = 0d0
537 90 : do j = 1, m
538 90 : its = its + 1
539 607590 : z = v(:,j)
540 90 : call psolve(n, z)
541 90 : call matvec(n, z, w)
542 370 : do i = 1, j
543 280 : h(i,j) = dotprd(n,w,v(1,i))
544 1890370 : w = w - h(i,j)*v(:,i)
545 : end do
546 90 : h(j+1,j) = dsqrt(dotprd(n,w,w))
547 607590 : v(:,j+1) = w / h(j+1,j)
548 :
549 90 : p(1) = beta
550 370 : p(2:j+1) = 0d0
551 1760 : h_(1:j+1,1:j) = h(1:j+1,1:j)
552 : !call dgelsy(j+1,j,1,h_(1:m+1, 1:j),m+1,p,m+1,piv,m,rank,work,4*m+1,i)
553 90 : call dgelsy(j+1,j,1,h_,m+1,p,m+1,piv,m,rank,work,4*m+1,i)
554 90 : if(i/=0) stop 'gmresm: dgelsy'
555 1980 : y = p
556 :
557 2680 : p(1:j+1) = - matmul(h(1:j+1,1:j),y(1:j))
558 90 : p(1) = p(1) + beta
559 460 : res = dsqrt(dot_product(p(1:j+1),p(1:j+1)))
560 : ! MPI aware: broadcast the 'res' value of master to avoid desynchronization.
561 90 : call xmpi_bcast(res, 0, xmpi_world, ierr)
562 : !if(info==1) print*, 'gmresm: it=', its,' res=', real(res)
563 90 : write(msg, *)'gmresm: it=', its,' res=', real(res)
564 90 : if(info==1) call wrtout(std_out, msg)
565 :
566 90 : done = (res<=tol .or. its==imx .or. res>res_)
567 72 : if(done .or. j==m) then
568 18 : if(del>0d0) call hookstep(j,h,m,beta,del, y)
569 729108 : z = matmul(v(:,1:j),y(1:j))
570 18 : call psolve(n, z)
571 121518 : x = x + z
572 18 : if(its==imx) info = 2
573 18 : if(res>res_) info = 1
574 18 : if(res<=tol) info = 0
575 18 : if(done) return
576 : !if(del>0d0) print*, 'gmres: warning! restart affects hookstep'
577 0 : if(del>0d0) call wrtout(std_out, 'gmres: warning! restart affects hookstep')
578 : goto 1 ! (j==m) restart
579 : end if
580 72 : res_ = res*stgn
581 :
582 : end do
583 :
584 : end subroutine gmresm
585 :
586 :
587 : !-----------------------------------------------------------------
588 : ! replace y with a vector that generates a hookstep
589 : ! c.f. Viswanath (2008) arXiv:0809.1498
590 : !-----------------------------------------------------------------
591 0 : subroutine hookstep(j,h,m,beta,del, y)
592 : implicit none
593 : integer, intent(in) :: j, m
594 : real(dp), intent(in) :: h(m+1,j), beta
595 : real(dp), intent(inout) :: del
596 : real(dp), intent(out) :: y(j)
597 0 : real(dp) :: a(j+1,j), s(j), u(j+1,j+1), vt(j,j), work(5*(j+1))
598 0 : real(dp) :: p(j+1), q(j), mu, qn
599 : integer :: info
600 :
601 0 : a = h(1:j+1,1:j)
602 :
603 0 : call dgesvd('A','A',j+1,j,a,j+1,s,u,j+1,vt,j,work,5*(j+1),info)
604 0 : if(info/=0) stop 'hookstep: dgesvd'
605 :
606 0 : p(1:j) = beta * u(1,1:j)
607 :
608 0 : mu = max(s(j)*s(j)*1d-6,1d-99)
609 0 : qn = 1d99
610 0 : do while(qn>del)
611 0 : mu = mu * 1.1d0
612 0 : q = p(1:j)*s/(mu+s*s)
613 0 : qn = dsqrt(dot_product(q,q))
614 : end do
615 :
616 0 : y = matmul(q,vt)
617 :
618 0 : p = - matmul(h(1:j+1,1:j),y(1:j))
619 0 : p(1) = p(1) + beta
620 0 : del = dsqrt(dot_product(p,p))
621 :
622 0 : end subroutine hookstep
623 :
624 90 : end module m_iterative_solvers
625 : !!***
|
