ddajac.f

Go to the documentation of this file.

      SUBROUTINE ddajac (NEQ, X, Y, YPRIME, DELTA, CJ, H,
     +   ier, wt, e, wm, iwm, res, ires, uround, jac, rpar,
     +   ipar, ntemp)
C***BEGIN PROLOGUE  DDAJAC
C***SUBSIDIARY
C***PURPOSE  Compute the iteration matrix for DDASSL and form the
C            LU-decomposition.
C***LIBRARY   SLATEC (DASSL)
C***TYPE      DOUBLE PRECISION (SDAJAC-S, DDAJAC-D)
C***AUTHOR  PETZOLD, LINDA R., (LLNL)
C***DESCRIPTION
C-----------------------------------------------------------------------
C     THIS ROUTINE COMPUTES THE ITERATION MATRIX
C     PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0).
C     HERE PD IS COMPUTED BY THE USER-SUPPLIED
C     ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND
C     IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING
C     IF IWM(MTYPE)IS 2 OR 5
C     THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
C     Y        = ARRAY CONTAINING PREDICTED VALUES
C     YPRIME   = ARRAY CONTAINING PREDICTED DERIVATIVES
C     DELTA    = RESIDUAL EVALUATED AT (X,Y,YPRIME)
C                (USED ONLY IF IWM(MTYPE)=2 OR 5)
C     CJ       = SCALAR PARAMETER DEFINING ITERATION MATRIX
C     H        = CURRENT STEPSIZE IN INTEGRATION
C     IER      = VARIABLE WHICH IS .NE. 0
C                IF ITERATION MATRIX IS SINGULAR,
C                AND 0 OTHERWISE.
C     WT       = VECTOR OF WEIGHTS FOR COMPUTING NORMS
C     E        = WORK SPACE (TEMPORARY) OF LENGTH NEQ
C     WM       = REAL WORK SPACE FOR MATRICES. ON
C                OUTPUT IT CONTAINS THE LU DECOMPOSITION
C                OF THE ITERATION MATRIX.
C     IWM      = INTEGER WORK SPACE CONTAINING
C                MATRIX INFORMATION
C     RES      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
C                TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
C     IRES     = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
C                IN RES, AND LESS THAN ZERO OTHERWISE.  (IF IRES
C                IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
C                IN THIS CASE (IF IRES .LT. 0), THEN IER = 0.
C     UROUND   = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED.
C     JAC      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
C                TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE
C                IS ONLY USED IF IWM(MTYPE) IS 1 OR 4)
C-----------------------------------------------------------------------
C***ROUTINES CALLED  DGBTRF, DGETRF
C***REVISION HISTORY  (YYMMDD)
C   830315  DATE WRITTEN
C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
C   901010  Modified three MAX calls to be all on one line.  (FNF)
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   901101  Corrected PURPOSE.  (FNF)
C   020204  Convert to use LAPACK
C***END PROLOGUE  DDAJAC
C
      INTEGER  neq, ier, iwm(*), ires, ipar(*), ntemp
      double precision
     *   x, y(*), yprime(*), delta(*), cj, h, wt(*), e(*), wm(*),
     *   uround, rpar(*)
      EXTERNAL  res, jac
C
      EXTERNAL  dgbtrf, dgetrf
C
      INTEGER  i, i1, i2, ii, ipsave, isave, j, k, l, lenpd, lipvt,
     *   lml, lmtype, lmu, mba, mband, meb1, meband, msave, mtype, n,
     *   npd, npdm1, nrow
      DOUBLE PRECISION  del, delinv, squr, ypsave, ysave
C
      parameter(npd=1)
      parameter(lml=1)
      parameter(lmu=2)
      parameter(lmtype=4)
      parameter(lipvt=22)
C
C***FIRST EXECUTABLE STATEMENT  DDAJAC
      ier = 0
      npdm1=npd-1
      mtype=iwm(lmtype)
      go to(100,200,300,400,500),mtype
C
C
C     DENSE USER-SUPPLIED MATRIX
100   lenpd=neq*neq
      DO 110 i=1,lenpd
110      wm(npdm1+i)=0.0d0
      CALL jac(x,y,yprime,wm(npd),cj,rpar,ipar)
      go to 230
C
C
C     DENSE FINITE-DIFFERENCE-GENERATED MATRIX
200   ires=0
      nrow=npdm1
      squr = sqrt(uround)
      DO 210 i=1,neq
         del=squr*max(abs(y(i)),abs(h*yprime(i)),abs(wt(i)))
         del=sign(del,h*yprime(i))
         del=(y(i)+del)-y(i)
         ysave=y(i)
         ypsave=yprime(i)
         y(i)=y(i)+del
         yprime(i)=yprime(i)+cj*del
         CALL res(x,y,yprime,e,ires,rpar,ipar)
         IF (ires .LT. 0) RETURN
         delinv=1.0d0/del
         DO 220 l=1,neq
220      wm(nrow+l)=(e(l)-delta(l))*delinv
      nrow=nrow+neq
      y(i)=ysave
      yprime(i)=ypsave
210   CONTINUE
C
C
C     DO DENSE-MATRIX LU DECOMPOSITION ON PD
230      CALL dgetrf( neq, neq, wm(npd), neq, iwm(lipvt), ier)
      RETURN
C
C
C     DUMMY SECTION FOR IWM(MTYPE)=3
300   RETURN
C
C
C     BANDED USER-SUPPLIED MATRIX
400   lenpd=(2*iwm(lml)+iwm(lmu)+1)*neq
      DO 410 i=1,lenpd
410      wm(npdm1+i)=0.0d0
      CALL jac(x,y,yprime,wm(npd),cj,rpar,ipar)
      meband=2*iwm(lml)+iwm(lmu)+1
      go to 550
C
C
C     BANDED FINITE-DIFFERENCE-GENERATED MATRIX
500   mband=iwm(lml)+iwm(lmu)+1
      mba=min(mband,neq)
      meband=mband+iwm(lml)
      meb1=meband-1
      msave=(neq/mband)+1
      isave=ntemp-1
      ipsave=isave+msave
      ires=0
      squr=sqrt(uround)
      DO 540 j=1,mba
         DO 510 n=j,neq,mband
          k= (n-j)/mband + 1
          wm(isave+k)=y(n)
          wm(ipsave+k)=yprime(n)
          del=squr*max(abs(y(n)),abs(h*yprime(n)),abs(wt(n)))
          del=sign(del,h*yprime(n))
          del=(y(n)+del)-y(n)
          y(n)=y(n)+del
510       yprime(n)=yprime(n)+cj*del
      CALL res(x,y,yprime,e,ires,rpar,ipar)
      IF (ires .LT. 0) RETURN
      DO 530 n=j,neq,mband
          k= (n-j)/mband + 1
          y(n)=wm(isave+k)
          yprime(n)=wm(ipsave+k)
          del=squr*max(abs(y(n)),abs(h*yprime(n)),abs(wt(n)))
          del=sign(del,h*yprime(n))
          del=(y(n)+del)-y(n)
          delinv=1.0d0/del
          i1=max(1,(n-iwm(lmu)))
          i2=min(neq,(n+iwm(lml)))
          ii=n*meb1-iwm(lml)+npdm1
          DO 520 i=i1,i2
520         wm(ii+i)=(e(i)-delta(i))*delinv
530      CONTINUE
540   CONTINUE
C
C
C     DO LU DECOMPOSITION OF BANDED PD
550   CALL dgbtrf(neq, neq, iwm(lml), iwm(lmu), wm(npd), meband,
     *     iwm(lipvt), ier)
      RETURN
C------END OF SUBROUTINE DDAJAC------
      END