dspigm Subroutine

subroutine dspigm(neq, t, y, ydot, savr, r, wght, maxl, kmp, epslin, cj, res, ires, nres, psol, npsol, z, v, hes, q, lgmr, rwp, iwp, wk, dl, rhok, iflag, irst, nrsts, rpar, ipar)

Uses

  • proc~~dspigm~~UsesGraph proc~dspigm dspigm module~blas_interfaces blas_interfaces proc~dspigm->module~blas_interfaces module~daskr daskr proc~dspigm->module~daskr module~daskr_kinds daskr_kinds proc~dspigm->module~daskr_kinds module~blas_interfaces->module~daskr_kinds module~daskr->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

This routine solves the linear system using a scaled preconditioned version of the generalized minimum residual method. An initial guess of is assumed.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: neq

Problem size.

real(kind=rk), intent(in) :: t

Current independent variable.

real(kind=rk), intent(in) :: y(neq)

Current dependent variables.

real(kind=rk), intent(in) :: ydot(neq)

Current derivatives of dependent variables.

real(kind=rk), intent(in) :: savr(neq)

Current residual evaluated at (t, y, ydot).

real(kind=rk), intent(inout) :: r(neq)

On entry, the right hand side vector. Also used as workspace when computing the final approximation and will therefore be destroyed. r is the same as v(:,maxl+1) in the call to this routine.

real(kind=rk), intent(in) :: wght(neq)

Nonzero elements of the diagonal scaling matrix.

integer, intent(in) :: maxl

Maximum allowable order of the matrix hes.

integer, intent(in) :: kmp

Number of previous vectors the new vector vnew must be made orthogonal to (kmp <= maxl).

real(kind=rk), intent(in) :: epslin

Tolerance on residuals in weighted rms norm.

real(kind=rk), intent(in) :: cj

Scalar used in forming the system Jacobian.

procedure(res_t) :: res

User-defined residuals routine.

integer, intent(out) :: ires

Error flag from res.

integer, intent(inout) :: nres

Number of calls to res.

procedure(psol_t) :: psol

User-defined preconditioner routine.

integer, intent(inout) :: npsol

Number of calls to psol.

real(kind=rk), intent(out) :: z(neq)

Final computed approximation to the solution of the system .

real(kind=rk), intent(out) :: v(neq,*)

Matrix of shape (neq,lgmr+1) containing the lgmr orthogonal vectors v(:,1) to v(:,lgmr).

real(kind=rk), intent(out) :: hes(maxl+1,maxl)

Upper triangular factor of the QR decomposition of the upper Hessenberg matrix whose entries are the scaled inner-products of A*v(:,i) and v(:,k).

real(kind=rk), intent(out) :: q(2*maxl)

Components of the Givens rotations used in the QR decomposition of hes. It is loaded in dheqr and used in dhels.

integer, intent(out) :: lgmr

The number of iterations performed and the current order of the upper Hessenberg matrix hes.

real(kind=rk), intent(inout) :: rwp(*)

Real work array used by preconditioner psol.

integer, intent(inout) :: iwp(*)

Integer work array used by preconditioner psol.

real(kind=rk), intent(inout) :: wk(*)

Real work array used by datv and psol.

real(kind=rk), intent(inout) :: dl(*)

Real work array used for calculation of the residual norm rhok when the method is incomplete (kmp < maxl) and/or when using restarting.

real(kind=rk), intent(out) :: rhok

Weighted norm of the final preconditioned residual.

integer, intent(out) :: iflag

Error flag. 0: convergence in lgmr iterations, lgmr <= maxl. 1: convergence test did not pass in maxl iterations, but the new residual norm (rho) is less than the old residual norm (rnrm), and so z is computed. 2: convergence test did not pass in maxl iterations, new residual norm (rho) is greater than or equal to old residual norm (rnrm), and the initial guess, z = 0, is returned. 3: recoverable error in psol caused by the preconditioner being out of date. -1: unrecoverable error in psol.

integer, intent(in) :: irst

Restarting flag. If irst > 0, then restarting is being performed.

integer, intent(in) :: nrsts

Number of restarts on the current call to dspigm. If nrsts > 0, then the residual r is already scaled, and so scaling of r is not necessary.

real(kind=rk), intent(inout) :: rpar(*)

User real workspace.

integer, intent(inout) :: ipar(*)

User integer workspace.


Calls

proc~~dspigm~~CallsGraph proc~dspigm dspigm datv datv proc~dspigm->datv dhels dhels proc~dspigm->dhels dheqr dheqr proc~dspigm->dheqr dorth dorth proc~dspigm->dorth interface~daxpy daxpy proc~dspigm->interface~daxpy interface~dcopy dcopy proc~dspigm->interface~dcopy interface~dnrm2 dnrm2 proc~dspigm->interface~dnrm2 interface~dscal dscal proc~dspigm->interface~dscal

Variables

Type Visibility Attributes Name Initial
integer, public :: i
integer, public :: ierr
integer, public :: info
integer, public :: ip1
integer, public :: i2
integer, public :: k
integer, public :: ll
integer, public :: llp1
integer, public :: maxlm1
integer, public :: maxlp1
real(kind=rk), public :: c
real(kind=rk), public :: dlnorm
real(kind=rk), public :: prod
real(kind=rk), public :: rho
real(kind=rk), public :: rnorm
real(kind=rk), public :: s
real(kind=rk), public :: snormw
real(kind=rk), public :: tem

Source Code

subroutine dspigm( &
   neq, t, y, ydot, savr, r, wght, maxl, &
   kmp, epslin, cj, res, ires, nres, psol, npsol, z, v, &
   hes, q, lgmr, rwp, iwp, wk, dl, rhok, iflag, irst, nrsts, &
   rpar, ipar)
!! This routine solves the linear system \(A z = r\) using a scaled preconditioned version
!! of the generalized minimum residual method. An initial guess of \(z=0\) is assumed.

   use daskr_kinds, only: rk, zero, one
   use daskr
   use blas_interfaces, only: daxpy, dcopy, dnrm2, dscal

   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: t
      !! Current independent variable.
   real(rk), intent(in) :: y(neq)
      !! Current dependent variables.
   real(rk), intent(in) :: ydot(neq)
      !! Current derivatives of dependent variables.
   real(rk), intent(in) :: savr(neq)
      !! Current residual evaluated at `(t, y, ydot)`.
   real(rk), intent(inout) :: r(neq)
      !! On entry, the right hand side vector. Also used as workspace when computing
      !! the final approximation and will therefore be destroyed. `r` is the same as
      !! `v(:,maxl+1)` in the call to this routine.
   real(rk), intent(in) :: wght(neq)
      !! Nonzero elements of the diagonal scaling matrix.
   integer, intent(in) :: maxl
      !! Maximum allowable order of the matrix `hes`.
   integer, intent(in) :: kmp
      !! Number of previous vectors the new vector `vnew` must be made orthogonal to
      !! (`kmp <= maxl`).
   real(rk), intent(in) :: epslin
      !! Tolerance on residuals \(Az - r\) in weighted rms norm.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   integer, intent(inout) :: nres
      !! Number of calls to `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   integer, intent(inout) :: npsol
      !! Number of calls to `psol`.
   real(rk), intent(out) :: z(neq)
      !! Final computed approximation to the solution of the system \(Az = r\).
   real(rk), intent(out) :: v(neq, *)
      !! Matrix of shape `(neq,lgmr+1)` containing the `lgmr` orthogonal vectors `v(:,1)`
      !! to `v(:,lgmr)`.
   real(rk), intent(out) :: hes(maxl + 1, maxl)
      !! Upper triangular factor of the QR decomposition of the upper Hessenberg matrix
      !! whose entries are the scaled inner-products of `A*v(:,i)` and `v(:,k)`.
   real(rk), intent(out) :: q(2*maxl)
      !! Components of the Givens rotations used in the QR decomposition of `hes`. It is
      !! loaded in [[dheqr]] and used in [[dhels]].
   integer, intent(out) :: lgmr
      !! The number of iterations performed and the current order of the upper
      !! Hessenberg matrix `hes`.
   real(rk), intent(inout) :: rwp(*)
      !! Real work array used by preconditioner `psol`.
   integer, intent(inout) :: iwp(*)
      !! Integer work array used by preconditioner `psol`.
   real(rk), intent(inout) :: wk(*)
      !! Real work array used by [[datv]] and `psol`.
   real(rk), intent(inout) :: dl(*)
      !! Real work array used for calculation of the residual norm `rhok` when the
      !! method is incomplete (`kmp < maxl`) and/or when using restarting.
   real(rk), intent(out) :: rhok
      !! Weighted norm of the final preconditioned residual.
   integer, intent(out) :: iflag
      !! Error flag.
      !! `0`: convergence in `lgmr` iterations, `lgmr <= maxl`.
      !! `1`: convergence test did not pass in `maxl` iterations, but the new
      !! residual norm (`rho`) is less than the old residual norm (`rnrm`), and so `z`
      !! is computed.
      !! `2`: convergence test did not pass in `maxl` iterations, new residual
      !! norm (`rho`) is greater than or equal to old residual norm (`rnrm`), and the
      !! initial guess, `z = 0`, is returned.
      !! `3`: recoverable error in `psol` caused by the preconditioner being out of date.
      !! `-1`: unrecoverable error in `psol`.
   integer, intent(in) :: irst
      !! Restarting flag. If `irst > 0`, then restarting is being performed.
   integer, intent(in) :: nrsts
      !! Number of restarts on the current call to [[dspigm]]. If `nrsts > 0`, then the
      !! residual `r` is already scaled, and so scaling of `r` is not necessary.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   integer :: i, ierr, info, ip1, i2, k, ll, llp1, maxlm1, maxlp1
   real(rk) :: c, dlnorm, prod, rho, rnorm, s, snormw, tem

   ierr = 0
   iflag = 0
   lgmr = 0
   npsol = 0
   nres = 0

   ! The initial guess for Z is 0. The initial residual is therefore
   ! the vector R. Initialize Z to 0.
   z = zero

   ! Apply inverse of left preconditioner to vector R if NRSTS == 0.
   ! Form V(:,1), the scaled preconditioned right hand side.
   if (nrsts == 0) then
      call psol(neq, t, y, ydot, savr, wk, cj, wght, rwp, iwp, r, epslin, ierr, rpar, ipar)
      npsol = 1
      if (ierr /= 0) goto 300
      v(:, 1) = r*wght
   else
      v(:, 1) = r
   end if

   ! Calculate norm of scaled vector V(:,1) and normalize it
   ! If, however, the norm of V(:,1) (i.e. the norm of the preconditioned
   ! residual) is <= EPSLIN, then return with Z=0.
   rnorm = dnrm2(neq, v, 1)
   if (rnorm <= epslin) then
      rhok = rnorm
      return
   end if
   tem = one/rnorm
   call dscal(neq, tem, v(1, 1), 1)

   ! Zero out the HES array.
   hes = zero

   ! Main loop to compute the vectors V(:,2) to V(:,MAXL).
   ! The running product PROD is needed for the convergence test.
   prod = one
   maxlp1 = maxl + 1
   do ll = 1, maxl
      lgmr = ll

      ! Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is
      ! the matrix A with scaling and inverse preconditioner factors applied.
      call datv(neq, y, t, ydot, savr, v(1, ll), wght, z, &
                res, ires, psol, v(1, ll + 1), wk, rwp, iwp, cj, epslin, &
                ierr, nres, npsol, rpar, ipar)
      if (ires < 0) return
      if (ierr /= 0) goto 300

      ! Call routine DORTH to orthogonalize the new vector VNEW = V(:,LL+1).
      call dorth(v(1, ll + 1), v, hes, neq, ll, maxlp1, kmp, snormw)
      hes(ll + 1, ll) = snormw

      ! Call routine DHEQR to update the factors of HES.
      call dheqr(hes, maxlp1, ll, q, info, ll)
      if (info == ll) goto 120

      ! Update RHO, the estimate of the norm of the residual R - A*ZL.
      ! If KMP < MAXL, then the vectors V(:,1),...,V(:,LL+1) are not
      ! necessarily orthogonal for LL > KMP.  The vector DL must then
      ! be computed, and its norm used in the calculation of RHO.
      prod = prod*q(2*ll)
      rho = abs(prod*rnorm)
      if ((ll > kmp) .and. (kmp < maxl)) then
         if (ll == kmp + 1) then
            call dcopy(neq, v(1, 1), 1, dl, 1)
            do i = 1, kmp
               ip1 = i + 1
               i2 = i*2
               s = q(i2)
               c = q(i2 - 1)
               do k = 1, neq
                  dl(k) = s*dl(k) + c*v(k, ip1)
               end do
            end do
         end if
         s = q(2*ll)
         c = q(2*ll - 1)/snormw
         llp1 = ll + 1
         do k = 1, neq
            dl(k) = s*dl(k) + c*v(k, llp1)
         end do
         dlnorm = dnrm2(neq, dl, 1)
         rho = rho*dlnorm
      end if

      ! Test for convergence. If passed, compute approximation ZL.
      ! If failed and LL < MAXL, then continue iterating.
      if (rho <= epslin) goto 200
      if (ll == maxl) goto 100

      ! Rescale so that the norm of V(1,LL+1) is one.
      tem = one/snormw
      call dscal(neq, tem, v(1, ll + 1), 1)

   end do

100 continue
   if (rho < rnorm) goto 150

120 continue
   iflag = 2
   z = zero
   return

150 continue
   ! The tolerance was not met, but the residual norm was reduced.
   ! If performing restarting (IRST > 0) calculate the residual vector
   ! RL and store it in the DL array.  If the incomplete version is
   ! being used (KMP < MAXL) then DL has already been calculated.
   iflag = 1
   if (irst > 0) then

      if (kmp == maxl) then
         !  Calculate DL from the V(I)'s.
         call dcopy(neq, v(1, 1), 1, dl, 1)
         maxlm1 = maxl - 1
         do i = 1, maxlm1
            ip1 = i + 1
            i2 = i*2
            s = q(i2)
            c = q(i2 - 1)
            do k = 1, neq
               dl(k) = s*dl(k) + c*v(k, ip1)
            end do
         end do
         s = q(2*maxl)
         c = q(2*maxl - 1)/snormw
         do k = 1, neq
            dl(k) = s*dl(k) + c*v(k, maxlp1)
         end do
      end if

      ! Scale DL by RNRM*PROD to obtain the residual RL.
      tem = rnorm*prod
      call dscal(neq, tem, dl, 1)
   end if

   ! Compute the approximation ZL to the solution.
   ! Since the vector Z was used as work space, and the initial guess
   ! of the Newton correction is zero, Z must be reset to zero.
200 continue
   ll = lgmr
   llp1 = ll + 1
   r(1:llp1) = zero
   r(1) = rnorm
   call dhels(hes, maxlp1, ll, q, r)
   z = zero
   do i = 1, ll
      call daxpy(neq, r(i), v(1, i), 1, z, 1)
   end do
   z = z/wght
   ! Load RHO into RHOK.
   rhok = rho
   return

   ! This block handles error returns forced by routine PSOL.
300 continue
   if (ierr < 0) iflag = -1
   if (ierr > 0) iflag = 3

end subroutine dspigm