dmatd.f

Go to the documentation of this file.

C Work performed 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 dmatd(NEQ,X,Y,YPRIME,DELTA,CJ,H,IER,EWT,E,
     *                 WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR)
C
C***BEGIN PROLOGUE  DMATD
C***REFER TO  DDASPK
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***REVISION DATE  940701   (YYMMDD) (new LIPVT)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     This routine computes the iteration matrix
C     J = dG/dY+CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
C     Here J is computed by:
C       the user-supplied routine JACD if IWM(MTYPE) is 1 or 4, or
C       by numerical difference quotients if IWM(MTYPE) is 2 or 5.
C
C     The parameters have the following meanings.
C     X        = Independent variable.
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 if iteration matrix
C                is singular, and 0 otherwise.
C     EWT      = Vector of error weights for computing norms.
C     E        = Work space (temporary) of length NEQ.
C     WM       = Real work space for matrices.  On output
C                it contains the LU decomposition
C                of the iteration matrix.
C     IWM      = Integer work space containing
C                matrix information.
C     RES      = External user-supplied subroutine
C                to evaluate the residual.  See RES description
C                in DDASPK prologue.
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     JACD     = 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                See JAC description for the case INFO(12) = 0
C                in DDASPK prologue.
C     RPAR,IPAR= Real and integer parameter arrays that
C                are used for communication between the
C                calling program and external user routines.
C                They are not altered by DMATD.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   JACD, RES, DGETRF, DGBTRF
C
C***END PROLOGUE  DMATD
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),delta(*),ewt(*),e(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*)
      EXTERNAL  res, jacd
C
      parameter(lml=1, lmu=2, lmtype=4, lnre=12, lnpd=22, llciwp=30)
C
      lipvt = iwm(llciwp)
      ier = 0
      mtype=iwm(lmtype)
      GO TO (100,200,300,400,500),mtype
C
C
C     Dense user-supplied matrix.
C
100   lenpd=iwm(lnpd)
      DO 110 i=1,lenpd
110      wm(i)=0.0d0
      CALL jacd(x,y,yprime,wm,cj,rpar,ipar)
      GO TO 230
C
C
C     Dense finite-difference-generated matrix.
C
200   ires=0
      nrow=0
      squr = sqrt(uround)
      DO 210 i=1,neq
         del=squr*max(abs(y(i)),abs(h*yprime(i)),
     *     abs(1.d0/ewt(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
         iwm(lnre)=iwm(lnre)+1
         CALL res(x,y,yprime,cj,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 J.
C
230      CALL dgetrf( neq, neq, wm, neq, iwm(lipvt), ier)
      RETURN
C
C
C     Dummy section for IWM(MTYPE)=3.
C
300   RETURN
C
C
C     Banded user-supplied matrix.
C
400   lenpd=iwm(lnpd)
      DO 410 i=1,lenpd
410      wm(i)=0.0d0
      CALL jacd(x,y,yprime,wm,cj,rpar,ipar)
      meband=2*iwm(lml)+iwm(lmu)+1
      GO TO 550
C
C
C     Banded finite-difference-generated matrix.
C
500   mband=iwm(lml)+iwm(lmu)+1
      mba=min0(mband,neq)
      meband=mband+iwm(lml)
      meb1=meband-1
      msave=(neq/mband)+1
      isave=iwm(lnpd)
      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(1.d0/ewt(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
        iwm(lnre)=iwm(lnre)+1
        CALL res(x,y,yprime,cj,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(1.d0/ewt(n)))
          del=sign(del,h*yprime(n))
          del=(y(n)+del)-y(n)
          delinv=1.0d0/del
          i1=max0(1,(n-iwm(lmu)))
          i2=min0(neq,(n+iwm(lml)))
          ii=n*meb1-iwm(lml)
          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 J.
C
550   CALL dgbtrf(neq, neq, iwm(lml), iwm(lmu), wm, meband,
     *     iwm(lipvt), ier)
      RETURN
C
C------END OF SUBROUTINE DMATD------------------------------------------
      END