d7/d62/ddasik_8f_source.html

C Work perfored under the auspices of the U.S. Department of Energy

C by Lawrence Livermore National Laboratory under contract number

C W-7405-Eng-48.

C


      SUBROUTINE ddasik(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,WT,

     *   JSKIP,RPAR,IPAR,SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,

     *   EPLI,SQRTN,RSQRTN,EPCON,RATEMX,STPTOL,JFLG,

     *   ICNFLG,ICNSTR,IERNLS)

C

C***BEGIN PROLOGUE  DDASIK

C***REFER TO  DDASPK

C***DATE WRITTEN   941026   (YYMMDD)

C***REVISION DATE  950808   (YYMMDD)

C***REVISION DATE  951110   Removed unreachable block 390.

C

C

C-----------------------------------------------------------------------

C***DESCRIPTION

C

C

C     DDASIK solves a nonlinear system of algebraic equations of the

C     form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in

C     the initial conditions.

C

C     An initial value for Y and initial guess for YPRIME are input.

C

C     The method used is a Newton scheme with Krylov iteration and a

C     linesearch algorithm.

C

C     The parameters represent

C

C     X         -- Independent variable.

C     Y         -- Solution vector at x.

C     YPRIME    -- Derivative of solution vector.

C     NEQ       -- Number of equations to be integrated.

C     ICOPT     -- Initial condition option chosen (1 or 2).

C     ID        -- Array of dimension NEQ, which must be initialized

C                  if ICOPT = 1.  See DDASIC.

C     RES       -- External user-supplied subroutine

C                  to evaluate the residual.  See RES description

C                  in DDASPK prologue.

C     JACK     --  External user-supplied routine to update

C                  the preconditioner.  (This is optional).

C                  See JAC description for the case

C                  INFO(12) = 1 in the DDASPK prologue.

C     PSOL      -- External user-supplied routine to solve

C                  a linear system using preconditioning.

C                  (This is optional).  See explanation inside DDASPK.

C     H         -- Scaling factor for this initial condition calc.

C     WT        -- Vector of weights for error criterion.

C     JSKIP     -- input flag to signal if initial JAC call is to be

C                  skipped.  1 => skip the call, 0 => do not skip call.

C     RPAR,IPAR -- Real and integer arrays used for communication

C                  between the calling program and external user

C                  routines.  They are not altered within DASPK.

C     SAVR      -- Work vector for DDASIK of length NEQ.

C     DELTA     -- Work vector for DDASIK of length NEQ.

C     R         -- Work vector for DDASIK of length NEQ.

C     YIC,YPIC  -- Work vectors for DDASIK, each of length NEQ.

C     PWK       -- Work vector for DDASIK of length NEQ.

C     WM,IWM    -- Real and integer arrays storing

C                  matrix information for linear system

C                  solvers, and various other information.

C     CJ        -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).

C     UROUND    -- Unit roundoff.

C     EPLI      -- convergence test constant.

C                  See DDASPK prologue for more details.

C     SQRTN     -- Square root of NEQ.

C     RSQRTN    -- reciprical of square root of NEQ.

C     EPCON     -- Tolerance to test for convergence of the Newton

C                  iteration.

C     RATEMX    -- Maximum convergence rate for which Newton iteration

C                  is considered converging.

C     JFLG      -- Flag showing whether a Jacobian routine is supplied.

C     ICNFLG    -- Integer scalar.  If nonzero, then constraint

C                  violations in the proposed new approximate solution

C                  will be checked for, and the maximum step length

C                  will be adjusted accordingly.

C     ICNSTR    -- Integer array of length NEQ containing flags for

C                  checking constraints.

C     IERNLS    -- Error flag for nonlinear solver.

C                   0   ==> nonlinear solver converged.

C                   1,2 ==> recoverable error inside nonlinear solver.

C                           1 => retry with current Y, YPRIME

C                           2 => retry with original Y, YPRIME

C                  -1   ==> unrecoverable error in nonlinear solver.

C

C-----------------------------------------------------------------------

C

C***ROUTINES CALLED

C   RES, JACK, DNSIK, DCOPY

C

C***END PROLOGUE  DDASIK

C

C

      IMPLICIT DOUBLE PRECISION(a-h,o-z)

      dimension y(*),yprime(*),id(*),wt(*),icnstr(*)

      dimension savr(*),delta(*),r(*),yic(*),ypic(*),pwk(*)

      dimension wm(*),iwm(*), rpar(*),ipar(*)

      EXTERNAL res, jack, psol

C

      parameter(lnre=12, lnje=13, llocwp=29, llciwp=30)

      parameter(lmxnit=32, lmxnj=33)

C

C

C     Perform initializations.

C

      lwp = iwm(llocwp)

      liwp = iwm(llciwp)

      mxnit = iwm(lmxnit)

      mxnj = iwm(lmxnj)

      iernls = 0

      nj = 0

      eplin = epli*epcon

C

C     Call RES to initialize DELTA.

C

      ires = 0

      iwm(lnre) = iwm(lnre) + 1

      CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)

      IF (ires .LT. 0) GO TO 370

C

C     Looping point for updating the preconditioner.

C

 300  CONTINUE

C

C     Initialize all error flags to zero.

C

      ierpj = 0

      ires = 0

      iernew = 0

C

C     If a Jacobian routine was supplied, call it.

C

      IF (jflg .EQ. 1 .AND. jskip .EQ. 0) THEN

        nj = nj + 1

        iwm(lnje)=iwm(lnje)+1

        CALL jack (res, ires, neq, x, y, yprime, wt, delta, r, h, cj,

     *     wm(lwp), iwm(liwp), ierpj, rpar, ipar)

        IF (ires .LT. 0 .OR. ierpj .NE. 0) GO TO 370

        ENDIF

      jskip = 0

C

C     Call the nonlinear Newton solver for up to MXNIT iterations.

C

      CALL dnsik(x,y,yprime,neq,icopt,id,res,psol,wt,rpar,ipar,

     *   savr,delta,r,yic,ypic,pwk,wm,iwm,cj,sqrtn,rsqrtn,

     *   eplin,epcon,ratemx,mxnit,stptol,icnflg,icnstr,iernew)

C

      IF (iernew .EQ. 1 .AND. nj .LT. mxnj .AND. jflg .EQ. 1) THEN

C

C       Up to MXNIT iterations were done, the convergence rate is < 1,

C       a Jacobian routine is supplied, and the number of JACK calls

C       is less than MXNJ.

C       Copy the residual SAVR to DELTA, call JACK, and try again.

C

        CALL dcopy (neq,  savr, 1, delta, 1)

        GO TO 300

        ENDIF

C

      IF (iernew .NE. 0) GO TO 380

      RETURN

C

C

C     Unsuccessful exits from nonlinear solver.

C     Set IERNLS accordingly.

C

 370  iernls = 2

      IF (ires .LE. -2) iernls = -1

      RETURN

C

 380  iernls = min(iernew,2)

      RETURN

C

C----------------------- END OF SUBROUTINE DDASIK-----------------------


      END

min
charNDArray min(char d, const charNDArray &m)
Definition chNDArray.cc:207

ddasik
subroutine ddasik(x, y, yprime, neq, icopt, id, res, jack, psol, h, wt, jskip, rpar, ipar, savr, delta, r, yic, ypic, pwk, wm, iwm, cj, uround, epli, sqrtn, rsqrtn, epcon, ratemx, stptol, jflg, icnflg, icnstr, iernls)
Definition ddasik.f:9

dnsik
subroutine dnsik(x, y, yprime, neq, icopt, id, res, psol, wt, rpar, ipar, savr, delta, r, yic, ypic, pwk, wm, iwm, cj, sqrtn, rsqrtn, eplin, epcon, ratemx, maxit, stptol, icnflg, icnstr, iernew)
Definition dnsik.f:8