Branch data 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 : : real(dp) :: rcond
503 : 36 : real(dp), target :: w(n), z(n)
504 : 36 : real(dp) :: h_(m+1,m), y(m+1), p(m+1), work(4*m+1)
505 : 18 : integer :: imx, piv(m), rank, i
506 : : real(dp), save :: beta
507 : : integer, save :: j
508 : : logical :: done
509 : : integer :: ierr
510 : : character(len=500) :: msg
511 : :
512 [ - + ]: 18 : if(info==2) then
513 : 0 : call hookstep(j,h,m,beta,del, y)
514 [ # # # # : 0 : z = matmul(v(:,1:j),y(1:j))
# # ]
515 : 0 : call psolve(n, z)
516 [ # # ]: 0 : x = z
517 : 0 : info = 0
518 : 0 : return
519 : : end if
520 : :
521 : 18 : tol = res
522 : 18 : imx = its
523 : 18 : its = 0
524 [ + + + + ]: 2551896 : v = 0d0
525 [ + + ]: 378 : piv = 0
526 : 18 : rcond = 1.0d-14
527 : :
528 : : 1 continue
529 : 18 : res_ = 1d99
530 : 18 : stgn = 1d0 - 1d-14
531 : :
532 : 18 : beta = dsqrt(dotprd(n,x,x))
533 [ + - + + ]: 121518 : if(beta==0d0) w = 0d0
534 [ - + ]: 18 : if(beta/=0d0) call matvec(n,x, w)
535 [ + + ]: 121518 : w = b - w
536 : 18 : beta = dsqrt(dotprd(n,w,w))
537 [ + + ]: 121518 : v(:,1) = w / beta
538 : :
539 [ + + + + ]: 7938 : h = 0d0
540 [ + - ]: 90 : do j = 1, m
541 : 90 : its = its + 1
542 [ + + ]: 607590 : z = v(:,j)
543 : 90 : call psolve(n, z)
544 : 90 : call matvec(n, z, w)
545 [ + + ]: 370 : do i = 1, j
546 : 280 : h(i,j) = dotprd(n,w,v(1,i))
547 [ + + ]: 1890370 : w = w - h(i,j)*v(:,i)
548 : : end do
549 : 90 : h(j+1,j) = dsqrt(dotprd(n,w,w))
550 [ + + ]: 607590 : v(:,j+1) = w / h(j+1,j)
551 : :
552 : 90 : p(1) = beta
553 [ + + ]: 370 : p(2:j+1) = 0d0
554 [ + + + + ]: 1760 : h_(1:j+1,1:j) = h(1:j+1,1:j)
555 : : !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)
556 : 90 : call dgelsy(j+1,j,1,h_,m+1,p,m+1,piv,rcond,rank,work,4*m+1,i)
557 [ - + ]: 90 : if(i/=0) stop 'gmresm: dgelsy'
558 [ + + ]: 1980 : y = p
559 : :
560 [ + - + - : 2680 : p(1:j+1) = - matmul(h(1:j+1,1:j),y(1:j))
+ - - + -
+ + + + +
+ + + + +
- ]
561 : 90 : p(1) = p(1) + beta
562 [ + + ]: 460 : res = dsqrt(dot_product(p(1:j+1),p(1:j+1)))
563 : : ! MPI aware: broadcast the 'res' value of master to avoid desynchronization.
564 : 90 : call xmpi_bcast(res, 0, xmpi_world, ierr)
565 : : !if(info==1) print*, 'gmresm: it=', its,' res=', real(res)
566 : 90 : write(msg, *)'gmresm: it=', its,' res=', real(res)
567 [ - + ]: 90 : if(info==1) call wrtout(std_out, msg)
568 : :
569 [ + + + - : 90 : done = (res<=tol .or. its==imx .or. res>res_)
+ - ]
570 [ - + ]: 72 : if(done .or. j==m) then
571 [ - + ]: 18 : if(del>0d0) call hookstep(j,h,m,beta,del, y)
572 [ + + + + : 729108 : z = matmul(v(:,1:j),y(1:j))
+ + ]
573 : 18 : call psolve(n, z)
574 [ + + ]: 121518 : x = x + z
575 [ - + ]: 18 : if(its==imx) info = 2
576 [ - + ]: 18 : if(res>res_) info = 1
577 [ + - ]: 18 : if(res<=tol) info = 0
578 [ - + ]: 18 : if(done) return
579 : : !if(del>0d0) print*, 'gmres: warning! restart affects hookstep'
580 [ # # ]: 0 : if(del>0d0) call wrtout(std_out, 'gmres: warning! restart affects hookstep')
581 : : goto 1 ! (j==m) restart
582 : : end if
583 : 72 : res_ = res*stgn
584 : :
585 : : end do
586 : :
587 : : end subroutine gmresm
588 : :
589 : :
590 : : !-----------------------------------------------------------------
591 : : ! replace y with a vector that generates a hookstep
592 : : ! c.f. Viswanath (2008) arXiv:0809.1498
593 : : !-----------------------------------------------------------------
594 : 0 : subroutine hookstep(j,h,m,beta,del, y)
595 : : implicit none
596 : : integer, intent(in) :: j, m
597 : : real(dp), intent(in) :: h(m+1,j), beta
598 : : real(dp), intent(inout) :: del
599 : : real(dp), intent(out) :: y(j)
600 : 0 : real(dp) :: a(j+1,j), s(j), u(j+1,j+1), vt(j,j), work(5*(j+1))
601 : 0 : real(dp) :: p(j+1), q(j), mu, qn
602 : : integer :: info
603 : :
604 [ # # # # ]: 0 : a = h(1:j+1,1:j)
605 : :
606 : 0 : call dgesvd('A','A',j+1,j,a,j+1,s,u,j+1,vt,j,work,5*(j+1),info)
607 [ # # ]: 0 : if(info/=0) stop 'hookstep: dgesvd'
608 : :
609 [ # # ]: 0 : p(1:j) = beta * u(1,1:j)
610 : :
611 : 0 : mu = max(s(j)*s(j)*1d-6,1d-99)
612 : 0 : qn = 1d99
613 [ # # ]: 0 : do while(qn>del)
614 : 0 : mu = mu * 1.1d0
615 [ # # ]: 0 : q = p(1:j)*s/(mu+s*s)
616 [ # # ]: 0 : qn = dsqrt(dot_product(q,q))
617 : : end do
618 : :
619 [ # # # # : 0 : y = matmul(q,vt)
# # # # ]
620 : :
621 [ # # # # : 0 : p = - matmul(h(1:j+1,1:j),y(1:j))
# # # # #
# # # # #
# # # # #
# ]
622 : 0 : p(1) = p(1) + beta
623 [ # # ]: 0 : del = dsqrt(dot_product(p,p))
624 : :
625 : 0 : end subroutine hookstep
626 : :
627 [ - + - - ]: 90 : end module m_iterative_solvers
628 : : !!***
|
