A wraper for the nitsol implementation of the generalised minimal residual method (GMRES) for iteratively solving linear systems. This provides a more general interface not specifically intended for use with Newton iteration. It also uses Fortran 90 features to provide a more convenient call signature.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=r8), | intent(inout), | dimension(:) | :: | solution | On input, an initial guess of the solution to the linear system. On output, the iteratively determined solution. |
|
| procedure(mat_mult) | :: | lhs | The linear operator on the left hand side of the linear system. |
|||
| real(kind=r8), | intent(in), | dimension(:) | :: | rhs | The right hand side of the linear system being solved |
|
| real(kind=r8), | intent(out) | :: | resid_norm | GMRES residual norm on return |
||
| integer, | intent(out) | :: | flag | Termination flag. Values have the following meanings:
|
||
| integer, | intent(out), | optional | :: | nlhs | Number of evaluations of the left hand side of the system |
|
| integer, | intent(out), | optional | :: | nrpre | Number of evaluations of the right-preconditioner |
|
| integer, | intent(out), | optional | :: | nli | Number of iterations |
|
| real(kind=r8), | intent(in), | optional | :: | tol | The tolerance for the solution. Default is |
|
| procedure(mat_mult), | optional | :: | precond | A right-preconditioner which may be used to improve convergence of the solution. |
||
| real(kind=r8), | intent(inout), | optional | dimension(*) | :: | rpar | Parameter/work array passed to the |
| integer, | intent(inout), | optional | dimension(*) | :: | ipar | Parameter/work array passed to the |
| integer, | intent(in), | optional | :: | resid_update | Residual update flag. On GMRES restarts, the residual can
be updated using a linear combination ( |
|
| integer, | intent(in), | optional | :: | iter_max | Maximum allowable number of GMRES iterations. Default is 1000. |
|
| integer, | intent(in), | optional | :: | krylov_dim | Maximum Krylov subspace dimension; default 10. |
|
| procedure(dinpr_intr), | optional | :: | inner_prod | Inner-product routine, either user-supplied or BLAS |
||
| procedure(dnorm_intr), | optional | :: | norm | Norm routine, either user supplied or BLAS dnrm2. |
subroutine gmres_solve(solution, lhs, rhs, resid_norm, flag, nlhs, nrpre, &
nli, tol, precond, rpar, ipar, resid_update, &
iter_max, krylov_dim, inner_prod, norm)
!* Author: Chris MacMackin
! Date: March 2017
!
! A wraper for the
! [nitsol](http://users.wpi.edu/~walker/Papers/nitsol,SISC_19,1998,302-318.pdf)
! implementation of the generalised minimal residual method
! ([GMRES](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method))
! for iteratively solving linear systems. This provides a more
! general interface not specifically intended for use with Newton
! iteration. It also uses Fortran 90 features to provide a more
! convenient call signature.
!
real(r8), dimension(:), intent(inout) :: solution
!! On input, an initial guess of the solution to the linear
!! system. On output, the iteratively determined solution.
procedure(mat_mult) :: lhs
!! The linear operator on the left hand side of the linear
!! system.
real(r8), dimension(:), intent(in) :: rhs
!! The right hand side of the linear system being solved
real(r8), intent(out) :: resid_norm
!! GMRES residual norm on return
integer, intent(out) :: flag
!! Termination flag. Values have the following meanings:
!!
!!0
!!: normal termination: acceptable solution found
!!
!!1
!!: \(J\vec{v}\) failure
!!
!!2
!!: \(P^{-1}\vec{v}\) failure
!!
!!3
!!: Acceptable solution not found in `iter_max` GMRES iterations
!!
!!4
!!: Insignificant residual norm reduction of a cycle of `kdmax`
!! steps (stagnation) before an acceptable solution has been
!! found.
!!
!!5
!!: Dangerous ill-conditioning detected before an acceptable
!! solution has been found.
real(r8), intent(in), optional :: tol
!! The tolerance for the solution. Default is `size(solution) * 1e-8`.
integer, intent(out), optional :: nlhs
!! Number of evaluations of the left hand side of the system
integer, intent(out), optional :: nrpre
!! Number of evaluations of the right-preconditioner
integer, intent(out), optional :: nli
!! Number of iterations
procedure(mat_mult), optional :: precond
!! A right-preconditioner which may be used to improve
!! convergence of the solution.
real(r8), dimension(*), intent(inout), optional :: rpar
!! Parameter/work array passed to the `lhs` and `precond` routines.
integer, dimension(*), intent(inout), optional :: ipar
!! Parameter/work array passed to the `lhs` and `precond` routines
integer, intent(in), optional :: resid_update
!! Residual update flag. On GMRES restarts, the residual can
!! be updated using a linear combination (`iresup == 0`) or by
!! direct evaluation (`iresup == 1`). The first is cheap (one
!! n-vector saxpy) but may lose accuracy with extreme residual
!! reduction; the second retains accuracy better but costs one
!! \(J\vec{v}\) product. Default is 0.
integer, intent(in), optional :: iter_max
!! Maximum allowable number of GMRES iterations. Default is
!! 1000.
integer, intent(in), optional :: krylov_dim
!! Maximum Krylov subspace dimension; default 10.
procedure(dinpr_intr), optional :: inner_prod
!! Inner-product routine, either user-supplied or BLAS `ddot`.
procedure(dnorm_intr), optional :: norm
!! Norm routine, either user supplied or BLAS dnrm2.
integer :: npoints, preflag, resup, itmax, kdim, nfe, lnlhs, lnrpre, lnli
real(r8) :: eta
procedure(dinpr_intr), pointer :: dinpr
procedure(dnorm_intr), pointer :: dnorm
real(r8), dimension(:), allocatable, save :: xcur, svbig, svsml, w, rwork
real(r8), dimension(:,:), allocatable, save :: vv, rr
lnlhs = 0
lnrpre = 0
lnli = 0
npoints = size(solution)
if (present(precond)) then
preflag = 1
else
preflag = 0
end if
if (present(tol)) then
eta = tol
else
eta = 1.e-8_r8 * npoints
end if
if (present(resid_update)) then
resup = resid_update
else
resup = 0
end if
if (present(iter_max)) then
itmax = iter_max
else
itmax = 1000
end if
if (present(krylov_dim)) then
kdim = krylov_dim
else
kdim = 10
end if
if (present(inner_prod)) then
dinpr => inner_prod
else
dinpr => ddot
end if
if (present(norm)) then
dnorm => norm
else
dnorm => dnrm2
end if
xcur = solution
if (allocated(vv)) then
if (any(shape(vv) /= [npoints, kdim+1])) then
deallocate(vv)
allocate(vv(npoints, kdim+1))
end if
if (size(svbig) /= kdim) then
deallocate(rr)
deallocate(svbig)
deallocate(svsml)
deallocate(w)
allocate(rr(kdim, kdim))
allocate(svbig(kdim))
allocate(svsml(kdim))
allocate(w(kdim))
end if
if (size(rwork) /= npoints) then
deallocate(rwork)
allocate(rwork(npoints))
end if
else
allocate(vv(npoints, kdim+1))
allocate(rr(kdim, kdim))
allocate(svbig(kdim))
allocate(svsml(kdim))
allocate(w(kdim))
allocate(rwork(npoints))
end if
call nitgm2(npoints, xcur, -rhs, solution, eta, dummy_f, jacv, &
rpar, ipar, 1, preflag, itmax, resup, 1, nfe, &
lnlhs, lnrpre, lnli, kdim, kdim+1, vv, rr, svbig, &
svsml, w, rwork, resid_norm, dinpr, dnorm, flag)
if (present(nlhs)) nlhs = lnlhs
if (present(nrpre)) nrpre = lnrpre
if (present(nli)) nli = lnli
contains
subroutine jacv(n, xcur, fcur, ijob, v, z, rpar, ipar, itrmjv)
!! A wrapper on the user-provided routines for the linear
!! operator and preconditioner, to put them in the form NITSOL
!! expects.
!!
integer, intent(in) :: n
! Dimension of the problem
real(r8), dimension(n), intent(in) :: xcur
! Array of length `n` containing the current \(x\) value
real(r8), dimension(n), intent(in) :: fcur
! Array of length `n` containing the current \(f(x)\) value
integer, intent(in) :: ijob
! Integer flag indicating which product is desired. 0
! indicates \(z = J\vec{v}\). 1 indicates \(z = P^{-1}\vec{v}\).
real(r8), dimension(n), intent(in) :: v
! An array of length `n` to be multiplied
real(r8), dimension(n), intent(out) :: z
! An array of length n containing the desired product on
! output.
real(r8), dimension(*), intent(inout) :: rpar
! Parameter/work array
integer, dimension(*), intent(inout) :: ipar
! Parameter/work array
integer, intent(out) :: itrmjv
! Termination flag. 0 indcates normal termination, 1
! indicatesfailure to prodce \(J\vec{v}\), and 2 indicates
! failure to produce \(P^{-1}\vec{v}\)
logical :: success
itrmjv = 0
if (ijob == 0) then
z = lhs(v, xcur, fcur, rpar, ipar, success)
if (.not. success) itrmjv = 1
else if (ijob == 1) then
z = precond(v, xcur, fcur, rpar, ipar, success)
if (.not. success) itrmjv = 2
else
error stop ('`ijob` not equal to 0 or 1.')
end if
end subroutine jacv
end subroutine gmres_solve