db/d54/ddaini_8f_source.html

       SUBROUTINE ddaini (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,

      +   IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)

 C***BEGIN PROLOGUE  DDAINI

 C***SUBSIDIARY

 C***PURPOSE  Initialization routine for DDASSL.

 C***LIBRARY   SLATEC (DASSL)

 C***TYPE      DOUBLE PRECISION (SDAINI-S, DDAINI-D)

 C***AUTHOR  PETZOLD, LINDA R., (LLNL)

 C***DESCRIPTION

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

 C     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER

 C     WITH THE BACKWARD EULER METHOD, TO

 C     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE

 C     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO

 C     SOLVE THE CORRECTOR ITERATION.

 C

 C     THE INITIAL GUESS FOR YPRIME IS USED IN THE

 C     PREDICTION, AND IN FORMING THE ITERATION

 C     MATRIX, BUT IS NOT INVOLVED IN THE

 C     ERROR TEST. THIS MAY HAVE TROUBLE

 C     CONVERGING IF THE INITIAL GUESS IS NO

 C     GOOD, OR IF G(X,Y,YPRIME) DEPENDS

 C     NONLINEARLY ON YPRIME.

 C

 C     THE PARAMETERS REPRESENT:

 C     X --         INDEPENDENT VARIABLE

 C     Y --         SOLUTION VECTOR AT X

 C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR

 C     NEQ --       NUMBER OF EQUATIONS

 C     H --         STEPSIZE. IMDER MAY USE A STEPSIZE

 C                  SMALLER THAN H.

 C     WT --        VECTOR OF WEIGHTS FOR ERROR

 C                  CRITERION

 C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS

 C                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY

 C                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME

 C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS

 C                  THAT ARE NOT ALTERED BY DDAINI

 C     PHI --       WORK SPACE FOR DDAINI

 C     DELTA,E --   WORK SPACE FOR DDAINI

 C     WM,IWM --    REAL AND INTEGER ARRAYS STORING

 C                  MATRIX INFORMATION

 C

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

 C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV

 C***REVISION HISTORY  (YYMMDD)

 C   830315  DATE WRITTEN

 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)

 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.

 C   901026  Added explicit declarations for all variables and minor

 C           cosmetic changes to prologue.  (FNF)

 C   901030  Minor corrections to declarations.  (FNF)

 C***END PROLOGUE  DDAINI

 C

       INTEGER  NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP

       double precision

      *   x, y(*), yprime(*), h, wt(*), rpar(*), phi(neq,*), delta(*),

      *   e(*), wm(*), hmin, uround

       EXTERNAL  res, jac

 C

       EXTERNAL  ddajac, ddanrm, ddaslv

       DOUBLE PRECISION  DDANRM

 C

       INTEGER  I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF,

      *   nef, nsf

       double precision

      *   cj, damp, delnrm, err, oldnrm, r, rate, s, xold, ynorm

       LOGICAL  CONVGD

 C

       parameter(lnre=12)

       parameter(lnje=13)

 C

       DATA maxit/10/,mjac/5/

       DATA damp/0.75d0/

 C

 C

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

 C     BLOCK 1.

 C     INITIALIZATIONS.

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

 C

 C***FIRST EXECUTABLE STATEMENT  DDAINI

       idid=1

       nef=0

       ncf=0

       nsf=0

       xold=x

       ynorm=ddanrm(neq,y,wt,rpar,ipar)

 C

 C     SAVE Y AND YPRIME IN PHI

       DO 100 i=1,neq

          phi(i,1)=y(i)

 100      phi(i,2)=yprime(i)

 C

 C

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

 C     BLOCK 2.

 C     DO ONE BACKWARD EULER STEP.

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

 C

 C     SET UP FOR START OF CORRECTOR ITERATION

 200   cj=1.0d0/h

       x=x+h

 C

 C     PREDICT SOLUTION AND DERIVATIVE

       DO 250 i=1,neq

 250     y(i)=y(i)+h*yprime(i)

 C

       jcalc=-1

       m=0

       convgd=.true.

 C

 C

 C     CORRECTOR LOOP.

 300   iwm(lnre)=iwm(lnre)+1

       ires=0

 C

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

       IF (ires.LT.0) GO TO 430

 C

 C

 C     EVALUATE THE ITERATION MATRIX

       IF (jcalc.NE.-1) GO TO 310

       iwm(lnje)=iwm(lnje)+1

       jcalc=0

       CALL ddajac(neq,x,y,yprime,delta,cj,h,

      *   ier,wt,e,wm,iwm,res,ires,

      *   uround,jac,rpar,ipar,ntemp)

 C

       s=1000000.d0

       IF (ires.LT.0) GO TO 430

       IF (ier.NE.0) GO TO 430

       nsf=0

 C

 C

 C

 C     MULTIPLY RESIDUAL BY DAMPING FACTOR

 310   CONTINUE

       DO 320 i=1,neq

 320      delta(i)=delta(i)*damp

 C

 C     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)

 C     STORE THE CORRECTION IN DELTA

 C

       CALL ddaslv(neq,delta,wm,iwm)

 C

 C     UPDATE Y AND YPRIME

       DO 330 i=1,neq

          y(i)=y(i)-delta(i)

 330      yprime(i)=yprime(i)-cj*delta(i)

 C

 C     TEST FOR CONVERGENCE OF THE ITERATION.

 C

       delnrm=ddanrm(neq,delta,wt,rpar,ipar)

       IF (delnrm.LE.100.d0*uround*ynorm)

      *   GO TO 400

 C

       IF (m.GT.0) GO TO 340

          oldnrm=delnrm

          GO TO 350

 C

 340   rate=(delnrm/oldnrm)**(1.0d0/m)

       IF (rate.GT.0.90d0) GO TO 430

       s=rate/(1.0d0-rate)

 C

 350   IF (s*delnrm .LE. 0.33d0) GO TO 400

 C

 C

 C     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE

 C     M AND AND TEST WHETHER THE MAXIMUM

 C     NUMBER OF ITERATIONS HAVE BEEN TRIED.

 C     EVERY MJAC ITERATIONS, GET A NEW

 C     ITERATION MATRIX.

 C

       m=m+1

       IF (m.GE.maxit) GO TO 430

 C

       IF ((m/mjac)*mjac.EQ.m) jcalc=-1

       GO TO 300

 C

 C

 C     THE ITERATION HAS CONVERGED.

 C     CHECK NONNEGATIVITY CONSTRAINTS

 400   IF (nonneg.EQ.0) GO TO 450

       DO 410 i=1,neq

 410      delta(i)=min(y(i),0.0d0)

 C

       delnrm=ddanrm(neq,delta,wt,rpar,ipar)

       IF (delnrm.GT.0.33d0) GO TO 430

 C

       DO 420 i=1,neq

          y(i)=y(i)-delta(i)

 420      yprime(i)=yprime(i)-cj*delta(i)

       GO TO 450

 C

 C

 C     EXITS FROM CORRECTOR LOOP.

 430   convgd=.false.

 450   IF (.NOT.convgd) GO TO 600

 C

 C

 C

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

 C     BLOCK 3.

 C     THE CORRECTOR ITERATION CONVERGED.

 C     DO ERROR TEST.

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

 C

       DO 510 i=1,neq

 510      e(i)=y(i)-phi(i,1)

       err=ddanrm(neq,e,wt,rpar,ipar)

 C

       IF (err.LE.1.0d0) RETURN

 C

 C

 C

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

 C     BLOCK 4.

 C     THE BACKWARD EULER STEP FAILED. RESTORE X, Y

 C     AND YPRIME TO THEIR ORIGINAL VALUES.

 C     REDUCE STEPSIZE AND TRY AGAIN, IF

 C     POSSIBLE.

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

 C

 600   CONTINUE

       x = xold

       DO 610 i=1,neq

          y(i)=phi(i,1)

 610      yprime(i)=phi(i,2)

 C

       IF (convgd) GO TO 640

       IF (ier.EQ.0) GO TO 620

          nsf=nsf+1

          h=h*0.25d0

          IF (nsf.LT.3.AND.abs(h).GE.hmin) GO TO 690

          idid=-12

          RETURN

 620   IF (ires.GT.-2) GO TO 630

          idid=-12

          RETURN

 630   ncf=ncf+1

       h=h*0.25d0

       IF (ncf.LT.10.AND.abs(h).GE.hmin) GO TO 690

          idid=-12

          RETURN

 C

 640   nef=nef+1

       r=0.90d0/(2.0d0*err+0.0001d0)

       r=max(0.1d0,min(0.5d0,r))

       h=h*r

       IF (abs(h).GE.hmin.AND.nef.LT.10) GO TO 690

          idid=-12

          RETURN

 690      GO TO 200

 C

 C-------------END OF SUBROUTINE DDAINI----------------------

       END

max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230

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

ddaini
subroutine ddaini(X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
Definition: ddaini.f:3

ddajac
subroutine ddajac(NEQ, X, Y, YPRIME, DELTA, CJ, H, IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR, IPAR, NTEMP)
Definition: ddajac.f:4

ddanrm
double precision function ddanrm(NEQ, V, WT, RPAR, IPAR)
Definition: ddanrm.f:2

ddaslv
subroutine ddaslv(NEQ, DELTA, WM, IWM)
Definition: ddaslv.f:2

abs
static T abs(T x)
Definition: pr-output.cc:1678