A wraper for the nitsol implementation of the transpose-free quasi-minimal residual method (TFQMR) 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 | TFQMR 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 TFQMR 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 tfqmr_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 transpose-free quasi-minimal residual method
! ([TFQMR](https://epubs.siam.org/doi/pdf/10.1137/0914029))
! 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
!! TFQMR residual norm on return
integer, intent(out) :: flag
!! Termination flag. Values have the following meanings:
!!
!!0
!!: normal termination: acceptable step found
!!
!!1
!!: \(J\vec{v}\) failure in `nitjv`
!!
!!2
!!: \(P^{-1}\vec{v}\) failure in `nitjv`
!!
!!3
!!: Acceptable step not found in `iksmax` TFQMR iterations
!!
!!4
!!: TFQMR breakdown
!!
!!5
!!: Floating point error (the underlying CGS iteration
!! has probably blown up)
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. Meaningless in this routine, but kept
!! for consistent interface.
integer, intent(in), optional :: iter_max
!! Maximum allowable number of TFQMR 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
real(r8), dimension(:,:), allocatable, save :: workers
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(workers)) then
if (size(workers,1) /= npoints) then
deallocate(workers)
allocate(workers(npoints, 11))
end if
else
allocate(workers(npoints, 11))
end if
call nittfq2(npoints, xcur, -rhs, solution, eta, dummy_f, jacv, &
rpar, ipar, 1, preflag, itmax, 1, nfe, lnlhs, lnrpre, &
lnli, workers(:,1), workers(:,2), workers(:,3), &
workers(:,4), workers(:,5), workers(:,6), &
workers(:,6), workers(:,8), workers(:,9), &
workers(:,10), workers(:,11), 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 tfqmr_solve