dvnlsd Subroutine

private subroutine dvnlsd(me, y, yh, ldyh, vsav, savf, ewt, acor, iwm, wm, nflag)

subroutine dvnlsd is a nonlinear system solver, which uses functional iteration or a chord (modified newton) method. for the chord method direct linear algebraic system solvers are used. subroutine dvnlsd then handles the corrector phase of this integration package.

Type Bound

dvode_t

Arguments

Type IntentOptional Attributes Name
class(dvode_t), intent(inout) :: me
real(kind=wp), intent(inout) :: y(*)

the dependent variable, a vector of length n
real(kind=wp), intent(inout) :: yh(ldyh,*)

the nordsieck (taylor) array, ldyh by lmax, input and output. on input, it contains predicted values.
integer, intent(in) :: ldyh

a constant >= n, the first dimension of yh
real(kind=wp) :: vsav(*)

unused work array.
real(kind=wp) :: savf(*)

a work array of length n.
real(kind=wp), intent(in) :: ewt(*)

an error weight vector of length n
real(kind=wp), intent(inout) :: acor(*)

a work array of length n, used for the accumulated corrections to the predicted y vector.
integer, intent(inout) :: iwm(*)

integer work array associated with matrix operations in chord iteration (miter /= 0).
real(kind=wp), intent(inout) :: wm(*)

real work array associated with matrix operations in chord iteration (miter /= 0).
integer, intent(inout) :: nflag

input/output flag, with values and meanings as follows:

  • input:

    • 0 first call for this time step.
    • -1 convergence failure in previous call to dvnlsd.
    • -2 error test failure in dvstep.

  • output:

    • 0 successful completion of nonlinear solver.
    • -1 convergence failure or singular matrix.
    • -2 unrecoverable error in matrix preprocessing (cannot occur here).
    • -3 unrecoverable error in solution (cannot occur here).

Calls

proc~~dvnlsd~~CallsGraph proc~dvnlsd dvode_module::dvode_t%dvnlsd proc~daxpy dvode_blas_module::daxpy proc~dvnlsd->proc~daxpy proc~dcopy dvode_blas_module::dcopy proc~dvnlsd->proc~dcopy proc~dscal dvode_blas_module::dscal proc~dvnlsd->proc~dscal proc~dvjac dvode_module::dvode_t%dvjac proc~dvnlsd->proc~dvjac proc~dvsol dvode_module::dvode_t%dvsol proc~dvnlsd->proc~dvsol proc~dvjac->proc~dcopy proc~dvjac->proc~dscal proc~dacopy dvode_module::dacopy proc~dvjac->proc~dacopy proc~dgbfa dvode_linpack_module::dgbfa proc~dvjac->proc~dgbfa proc~dgefa dvode_linpack_module::dgefa proc~dvjac->proc~dgefa proc~dgbsl dvode_linpack_module::dgbsl proc~dvsol->proc~dgbsl proc~dgesl dvode_linpack_module::dgesl proc~dvsol->proc~dgesl proc~dacopy->proc~dcopy proc~dgbfa->proc~daxpy proc~dgbfa->proc~dscal proc~idamax dvode_blas_module::idamax proc~dgbfa->proc~idamax proc~dgbsl->proc~daxpy proc~ddot dvode_blas_module::ddot proc~dgbsl->proc~ddot proc~dgefa->proc~daxpy proc~dgefa->proc~dscal proc~dgefa->proc~idamax proc~dgesl->proc~daxpy proc~dgesl->proc~ddot

Called by

proc~~dvnlsd~~CalledByGraph proc~dvnlsd dvode_module::dvode_t%dvnlsd proc~dvstep dvode_module::dvode_t%dvstep proc~dvstep->proc~dvnlsd proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvstep

Source Code

   subroutine dvnlsd(me,y,yh,ldyh,vsav,savf,ewt,acor,iwm,wm,nflag)

      class(dvode_t),intent(inout) :: me
      real(wp),intent(inout) :: y(*) !! the dependent variable, a vector of length n
      integer,intent(in) :: ldyh !! a constant >= n, the first dimension of yh
      real(wp),intent(inout) :: yh(ldyh,*) !! the nordsieck (taylor) array, ldyh by lmax, input
                                           !! and output.  on input, it contains predicted values.
      real(wp) :: vsav(*) !! unused work array.
      real(wp) :: savf(*) !! a work array of length n.
      real(wp),intent(in) :: ewt(*) !! an error weight vector of length n
      real(wp),intent(inout) :: acor(*) !! a work array of length n, used for the accumulated
                                        !! corrections to the predicted y vector.
      real(wp),intent(inout) :: wm(*) !! real work array associated with matrix
                                      !! operations in chord iteration (miter /= 0).
      integer,intent(inout) :: iwm(*) !! integer work array associated with matrix
                                      !! operations in chord iteration (miter /= 0).
      integer,intent(inout) :: nflag !! input/output flag, with values and meanings as follows:
                                     !!
                                     !! * **input:**
                                     !!     *  0 first call for this time step.
                                     !!     * -1 convergence failure in previous call to dvnlsd.
                                     !!     * -2 error test failure in dvstep.
                                     !!
                                     !! * **output:**
                                     !!     *  0 successful completion of nonlinear solver.
                                     !!     * -1 convergence failure or singular matrix.
                                     !!     * -2 unrecoverable error in matrix preprocessing (cannot occur here).
                                     !!     * -3 unrecoverable error in solution (cannot occur here).

      real(wp) :: cscale , dcon , del , delp
      integer :: i , ierpj , iersl , m

      integer,parameter :: maxcor = 3
      integer,parameter :: msbp = 20
      real(wp),parameter :: ccmax = 0.3_wp
      real(wp),parameter :: crdown = 0.3_wp
      real(wp),parameter :: rdiv = 2.0_wp
      real(wp),parameter :: one = 1.0_wp
      real(wp),parameter :: two = 2.0_wp

      !-----------------------------------------------------------------------
      ! on the first step, on a change of method order, or after a
      ! nonlinear convergence failure with nflag = -2, set ipup = miter
      ! to force a jacobian update when miter /= 0.
      !-----------------------------------------------------------------------
      if ( me%dat%jstart==0 ) me%dat%nslp = 0
      if ( nflag==0 ) me%dat%icf = 0
      if ( nflag==-2 ) me%dat%ipup = me%dat%miter
      if ( (me%dat%jstart==0) .or. (me%dat%jstart==-1) ) me%dat%ipup = me%dat%miter
      ! if this is functional iteration, set crate == 1 and drop to 220
      if ( me%dat%miter==0 ) then
         me%dat%crate = one
      else
         !-----------------------------------------------------------------------
         ! rc is the ratio of new to old values of the coefficient h/el(2)=h/l1.
         ! when rc differs from 1 by more than ccmax, ipup is set to miter
         ! to force dvjac to be called, if a jacobian is involved.
         ! in any case, dvjac is called at least every msbp steps.
         !-----------------------------------------------------------------------
         me%dat%drc = abs(me%dat%rc-one)
         if ( me%dat%drc>ccmax .or. me%dat%nst>=me%dat%nslp+msbp ) me%dat%ipup = me%dat%miter
      end if

      corrector : do
         !-----------------------------------------------------------------------
         ! up to maxcor corrector iterations are taken.  a convergence test is
         ! made on the r.m.s. norm of each correction, weighted by the error
         ! weight vector ewt.  the sum of the corrections is accumulated in the
         ! vector acor(i).  the yh array is not altered in the corrector loop.
         !-----------------------------------------------------------------------
         m = 0
         delp = zero
         call dcopy(me%dat%n,yh(1,1),1,y,1)
         call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),savf(1:me%dat%n))
         me%dat%nfe = me%dat%nfe + 1
         if ( me%dat%ipup>0 ) then
            !-----------------------------------------------------------------------
            ! if indicated, the matrix p = i - h*rl1*j is reevaluated and
            ! preprocessed before starting the corrector iteration.  ipup is set
            ! to 0 as an indicator that this has been done.
            !-----------------------------------------------------------------------
            call me%dvjac(y,yh,ldyh,ewt,acor,savf,wm,iwm,ierpj)
            me%dat%ipup = 0
            me%dat%rc = one
            me%dat%drc = zero
            me%dat%crate = one
            me%dat%nslp = me%dat%nst
            ! if matrix is singular, take error return to force cut in step size. --
            if ( ierpj/=0 ) exit corrector
         endif
         do i = 1 , me%dat%n
            acor(i) = zero
         enddo

         do
            ! this is a looping point for the corrector iteration. -----------------
            if ( me%dat%miter/=0 ) then
               !-----------------------------------------------------------------------
               ! in the case of the chord method, compute the corrector error,
               ! and solve the linear system with that as right-hand side and
               ! p as coefficient matrix.  the correction is scaled by the factor
               ! 2/(1+rc) to account for changes in h*rl1 since the last dvjac call.
               !-----------------------------------------------------------------------
               do i = 1 , me%dat%n
                  y(i) = (me%dat%rl1*me%dat%h)*savf(i) - (me%dat%rl1*yh(i,2)+acor(i))
               enddo
               call me%dvsol(wm,iwm,y,iersl)
               me%dat%nni = me%dat%nni + 1
               if ( iersl>0 ) exit
               if ( me%dat%meth==2 .and. me%dat%rc/=one ) then
                  cscale = two/(one+me%dat%rc)
                  call dscal(me%dat%n,cscale,y,1)
               endif
               del = me%dvnorm(me%dat%n,y(1:me%dat%n),ewt(1:me%dat%n))
               call daxpy(me%dat%n,one,y,1,acor,1)
               do i = 1 , me%dat%n
                  y(i) = yh(i,1) + acor(i)
               enddo
            else
               !-----------------------------------------------------------------------
               ! in the case of functional iteration, update y directly from
               ! the result of the last function evaluation.
               !-----------------------------------------------------------------------
               do i = 1 , me%dat%n
                  savf(i) = me%dat%rl1*(me%dat%h*savf(i)-yh(i,2))
               enddo
               do i = 1 , me%dat%n
                  y(i) = savf(i) - acor(i)
               enddo
               del = me%dvnorm(me%dat%n,y(1:me%dat%n),ewt(1:me%dat%n))
               do i = 1 , me%dat%n
                  y(i) = yh(i,1) + savf(i)
               enddo
               call dcopy(me%dat%n,savf,1,acor,1)
            endif
            !-----------------------------------------------------------------------
            ! test for convergence.  if m > 0, an estimate of the convergence
            ! rate constant is stored in crate, and this is used in the test.
            !-----------------------------------------------------------------------
            if ( m/=0 ) me%dat%crate = max(crdown*me%dat%crate,del/delp)
            dcon = del*min(one,me%dat%crate)/me%dat%tq(4)
            if ( dcon<=one ) then
               ! return for successful step. ------------------------------------------
               nflag = 0
               me%dat%jcur = 0
               me%dat%icf = 0
               if ( m==0 ) me%dat%acnrm = del
               if ( m>0 ) me%dat%acnrm = me%dvnorm(me%dat%n,acor(1:me%dat%n),ewt(1:me%dat%n))
               return
            else
               m = m + 1
               if ( m/=maxcor ) then
                  if ( m<2 .or. del<=rdiv*delp ) then
                     delp = del
                     call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),savf(1:me%dat%n))
                     me%dat%nfe = me%dat%nfe + 1
                     cycle
                  endif
               endif
               exit
            endif
         end do

         if ( me%dat%miter/=0 .and. me%dat%jcur/=1 ) then
            me%dat%icf = 1
            me%dat%ipup = me%dat%miter
         else
            exit corrector
         endif

      end do corrector

      nflag = -1
      me%dat%icf = 2
      me%dat%ipup = me%dat%miter

   end subroutine dvnlsd