ddajac.f

Go to the documentation of this file.
00001       SUBROUTINE DDAJAC (NEQ, X, Y, YPRIME, DELTA, CJ, H,
00002      +   IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR,
00003      +   IPAR, NTEMP)
00004 C***BEGIN PROLOGUE  DDAJAC
00005 C***SUBSIDIARY
00006 C***PURPOSE  Compute the iteration matrix for DDASSL and form the
00007 C            LU-decomposition.
00008 C***LIBRARY   SLATEC (DASSL)
00009 C***TYPE      DOUBLE PRECISION (SDAJAC-S, DDAJAC-D)
00010 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
00011 C***DESCRIPTION
00012 C-----------------------------------------------------------------------
00013 C     THIS ROUTINE COMPUTES THE ITERATION MATRIX
00014 C     PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0).
00015 C     HERE PD IS COMPUTED BY THE USER-SUPPLIED
00016 C     ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND
00017 C     IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING
00018 C     IF IWM(MTYPE)IS 2 OR 5
00019 C     THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
00020 C     Y        = ARRAY CONTAINING PREDICTED VALUES
00021 C     YPRIME   = ARRAY CONTAINING PREDICTED DERIVATIVES
00022 C     DELTA    = RESIDUAL EVALUATED AT (X,Y,YPRIME)
00023 C                (USED ONLY IF IWM(MTYPE)=2 OR 5)
00024 C     CJ       = SCALAR PARAMETER DEFINING ITERATION MATRIX
00025 C     H        = CURRENT STEPSIZE IN INTEGRATION
00026 C     IER      = VARIABLE WHICH IS .NE. 0
00027 C                IF ITERATION MATRIX IS SINGULAR,
00028 C                AND 0 OTHERWISE.
00029 C     WT       = VECTOR OF WEIGHTS FOR COMPUTING NORMS
00030 C     E        = WORK SPACE (TEMPORARY) OF LENGTH NEQ
00031 C     WM       = REAL WORK SPACE FOR MATRICES. ON
00032 C                OUTPUT IT CONTAINS THE LU DECOMPOSITION
00033 C                OF THE ITERATION MATRIX.
00034 C     IWM      = INTEGER WORK SPACE CONTAINING
00035 C                MATRIX INFORMATION
00036 C     RES      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
00037 C                TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
00038 C     IRES     = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
00039 C                IN RES, AND LESS THAN ZERO OTHERWISE.  (IF IRES
00040 C                IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
00041 C                IN THIS CASE (IF IRES .LT. 0), THEN IER = 0.
00042 C     UROUND   = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED.
00043 C     JAC      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
00044 C                TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE
00045 C                IS ONLY USED IF IWM(MTYPE) IS 1 OR 4)
00046 C-----------------------------------------------------------------------
00047 C***ROUTINES CALLED  DGBTRF, DGETRF
00048 C***REVISION HISTORY  (YYMMDD)
00049 C   830315  DATE WRITTEN
00050 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
00051 C   901010  Modified three MAX calls to be all on one line.  (FNF)
00052 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
00053 C   901026  Added explicit declarations for all variables and minor
00054 C           cosmetic changes to prologue.  (FNF)
00055 C   901101  Corrected PURPOSE.  (FNF)
00056 C   020204  Convert to use LAPACK
00057 C***END PROLOGUE  DDAJAC
00058 C
00059       INTEGER  NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP
00060       DOUBLE PRECISION
00061      *   X, Y(*), YPRIME(*), DELTA(*), CJ, H, WT(*), E(*), WM(*),
00062      *   UROUND, RPAR(*)
00063       EXTERNAL  RES, JAC
00064 C
00065       EXTERNAL  DGBTRF, DGETRF
00066 C
00067       INTEGER  I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT,
00068      *   LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N,
00069      *   NPD, NPDM1, NROW
00070       DOUBLE PRECISION  DEL, DELINV, SQUR, YPSAVE, YSAVE
00071 C
00072       PARAMETER (NPD=1)
00073       PARAMETER (LML=1)
00074       PARAMETER (LMU=2)
00075       PARAMETER (LMTYPE=4)
00076       PARAMETER (LIPVT=22)
00077 C
00078 C***FIRST EXECUTABLE STATEMENT  DDAJAC
00079       IER = 0
00080       NPDM1=NPD-1
00081       MTYPE=IWM(LMTYPE)
00082       GO TO (100,200,300,400,500),MTYPE
00083 C
00084 C
00085 C     DENSE USER-SUPPLIED MATRIX
00086 100   LENPD=NEQ*NEQ
00087       DO 110 I=1,LENPD
00088 110      WM(NPDM1+I)=0.0D0
00089       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
00090       GO TO 230
00091 C
00092 C
00093 C     DENSE FINITE-DIFFERENCE-GENERATED MATRIX
00094 200   IRES=0
00095       NROW=NPDM1
00096       SQUR = SQRT(UROUND)
00097       DO 210 I=1,NEQ
00098          DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),ABS(WT(I)))
00099          DEL=SIGN(DEL,H*YPRIME(I))
00100          DEL=(Y(I)+DEL)-Y(I)
00101          YSAVE=Y(I)
00102          YPSAVE=YPRIME(I)
00103          Y(I)=Y(I)+DEL
00104          YPRIME(I)=YPRIME(I)+CJ*DEL
00105          CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
00106          IF (IRES .LT. 0) RETURN
00107          DELINV=1.0D0/DEL
00108          DO 220 L=1,NEQ
00109 220      WM(NROW+L)=(E(L)-DELTA(L))*DELINV
00110       NROW=NROW+NEQ
00111       Y(I)=YSAVE
00112       YPRIME(I)=YPSAVE
00113 210   CONTINUE
00114 C
00115 C
00116 C     DO DENSE-MATRIX LU DECOMPOSITION ON PD
00117 230      CALL DGETRF( NEQ, NEQ, WM(NPD), NEQ, IWM(LIPVT), IER)
00118       RETURN
00119 C
00120 C
00121 C     DUMMY SECTION FOR IWM(MTYPE)=3
00122 300   RETURN
00123 C
00124 C
00125 C     BANDED USER-SUPPLIED MATRIX
00126 400   LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ
00127       DO 410 I=1,LENPD
00128 410      WM(NPDM1+I)=0.0D0
00129       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
00130       MEBAND=2*IWM(LML)+IWM(LMU)+1
00131       GO TO 550
00132 C
00133 C
00134 C     BANDED FINITE-DIFFERENCE-GENERATED MATRIX
00135 500   MBAND=IWM(LML)+IWM(LMU)+1
00136       MBA=MIN(MBAND,NEQ)
00137       MEBAND=MBAND+IWM(LML)
00138       MEB1=MEBAND-1
00139       MSAVE=(NEQ/MBAND)+1
00140       ISAVE=NTEMP-1
00141       IPSAVE=ISAVE+MSAVE
00142       IRES=0
00143       SQUR=SQRT(UROUND)
00144       DO 540 J=1,MBA
00145          DO 510 N=J,NEQ,MBAND
00146           K= (N-J)/MBAND + 1
00147           WM(ISAVE+K)=Y(N)
00148           WM(IPSAVE+K)=YPRIME(N)
00149           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
00150           DEL=SIGN(DEL,H*YPRIME(N))
00151           DEL=(Y(N)+DEL)-Y(N)
00152           Y(N)=Y(N)+DEL
00153 510       YPRIME(N)=YPRIME(N)+CJ*DEL
00154       CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
00155       IF (IRES .LT. 0) RETURN
00156       DO 530 N=J,NEQ,MBAND
00157           K= (N-J)/MBAND + 1
00158           Y(N)=WM(ISAVE+K)
00159           YPRIME(N)=WM(IPSAVE+K)
00160           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
00161           DEL=SIGN(DEL,H*YPRIME(N))
00162           DEL=(Y(N)+DEL)-Y(N)
00163           DELINV=1.0D0/DEL
00164           I1=MAX(1,(N-IWM(LMU)))
00165           I2=MIN(NEQ,(N+IWM(LML)))
00166           II=N*MEB1-IWM(LML)+NPDM1
00167           DO 520 I=I1,I2
00168 520         WM(II+I)=(E(I)-DELTA(I))*DELINV
00169 530      CONTINUE
00170 540   CONTINUE
00171 C
00172 C
00173 C     DO LU DECOMPOSITION OF BANDED PD
00174 550   CALL DGBTRF(NEQ, NEQ, IWM(LML), IWM(LMU), WM(NPD), MEBAND,
00175      *     IWM(LIPVT), IER)
00176       RETURN
00177 C------END OF SUBROUTINE DDAJAC------
00178       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines