home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / odepack / prepj.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  178 lines
  1.       SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
  2.      1   F, JAC, IERR)
  3. CLLL. OPTIMIZE
  4.       EXTERNAL F, JAC
  5.       INTEGER NEQ, NYH, IWM
  6.       INTEGER IOWND, IOWNS,
  7.      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  8.      2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  9.       INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
  10.      1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
  11.       DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM 
  12.       DOUBLE PRECISION ROWNS, 
  13.      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  14.       DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
  15.      1   VNORM
  16.       DIMENSION NEQ(1), Y(1), YH(NYH,1), EWT(1), FTEM(1), SAVF(1),
  17.      1   WM(1), IWM(1)
  18.       COMMON /LS0001/ ROWNS(209),
  19.      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  20.      3   IOWND(14), IOWNS(6), 
  21.      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  22.      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  23. C-----------------------------------------------------------------------
  24. C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
  25. C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
  26. C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
  27. C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
  28. C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
  29. C J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
  30. C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
  31. C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
  32. C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
  33. C
  34. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
  35. C WITH PREPJ USES THE FOLLOWING..
  36. C Y     = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
  37. C FTEM  = WORK ARRAY OF LENGTH N (ACOR IN STODE). 
  38. C SAVF  = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
  39. C WM    = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
  40. C         INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
  41. C         OF P IF MITER IS 1, 2 , 4, OR 5.
  42. C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
  43. C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
  44. C         WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
  45. C         WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
  46. C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
  47. C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
  48. C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
  49. C EL0   = EL(1) (INPUT).
  50. C IERPJ = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .GT. 0 IF
  51. C         P MATRIX FOUND TO BE SINGULAR.
  52. C JCUR  = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
  53. C         (OR APPROXIMATION) IS NOW CURRENT.
  54. C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
  55. C MITER, N, NFE, AND NJE.
  56. C-----------------------------------------------------------------------
  57.       NJE = NJE + 1 
  58.       IERPJ = 0
  59.       JCUR = 1
  60.       HL0 = H*EL0
  61.       GO TO (100, 200, 300, 400, 500), MITER
  62. C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
  63.  100  LENP = N*N
  64.       DO 110 I = 1,LENP
  65.  110    WM(I+2) = 0.0D0
  66.       CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
  67.       CON = -HL0
  68.       DO 120 I = 1,LENP
  69.  120    WM(I+2) = WM(I+2)*CON 
  70.       GO TO 240
  71. C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
  72.  200  FAC = VNORM (N, SAVF, EWT)
  73.       R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
  74.       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
  75.       SRUR = WM(1)
  76.       J1 = 2
  77.       DO 230 J = 1,N
  78.         YJ = Y(J)
  79.         R = DMAX1(SRUR*DABS(YJ),R0/EWT(J))
  80.         Y(J) = Y(J) + R
  81.         FAC = -HL0/R
  82.         IERR = 0
  83.         CALL F (NEQ, TN, Y, FTEM, IERR)
  84.         IF (IERR .LT. 0) RETURN
  85.         DO 220 I = 1,N
  86.  220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
  87.         Y(J) = YJ
  88.         J1 = J1 + N 
  89.  230    CONTINUE
  90.       NFE = NFE + N 
  91. C ADD IDENTITY MATRIX. -------------------------------------------------
  92.  240  J = 3
  93.       NP1 = N + 1
  94.       DO 250 I = 1,N
  95.         WM(J) = WM(J) + 1.0D0 
  96.  250    J = J + NP1 
  97. C DO LU DECOMPOSITION ON P. --------------------------------------------
  98.       CALL DGEFA (WM(3), N, N, IWM(21), IER)
  99.       IF (IER .NE. 0) IERPJ = 1
  100.       RETURN
  101. C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
  102.  300  WM(2) = HL0
  103.       R = EL0*0.1D0 
  104.       DO 310 I = 1,N
  105.  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
  106.       IERR = 0
  107.       CALL F (NEQ, TN, Y, WM(3), IERR)
  108.       IF (IERR .LT. 0) RETURN
  109.       NFE = NFE + 1 
  110.       DO 320 I = 1,N
  111.         R0 = H*SAVF(I) - YH(I,2)
  112.         DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
  113.         WM(I+2) = 1.0D0
  114.         IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
  115.         IF (DABS(DI) .EQ. 0.0D0) GO TO 330
  116.         WM(I+2) = 0.1D0*R0/DI 
  117.  320    CONTINUE
  118.       RETURN
  119.  330  IERPJ = 1
  120.       RETURN
  121. C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
  122.  400  ML = IWM(1)
  123.       MU = IWM(2)
  124.       ML3 = ML + 3
  125.       MBAND = ML + MU + 1
  126.       MEBAND = MBAND + ML
  127.       LENP = MEBAND*N
  128.       DO 410 I = 1,LENP
  129.  410    WM(I+2) = 0.0D0
  130.       CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
  131.       CON = -HL0
  132.       DO 420 I = 1,LENP
  133.  420    WM(I+2) = WM(I+2)*CON 
  134.       GO TO 570
  135. C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
  136.  500  ML = IWM(1)
  137.       MU = IWM(2)
  138.       MBAND = ML + MU + 1
  139.       MBA = MIN0(MBAND,N)
  140.       MEBAND = MBAND + ML
  141.       MEB1 = MEBAND - 1
  142.       SRUR = WM(1)
  143.       FAC = VNORM (N, SAVF, EWT)
  144.       R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
  145.       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
  146.       DO 560 J = 1,MBA
  147.         DO 530 I = J,N,MBAND
  148.           YI = Y(I) 
  149.           R = DMAX1(SRUR*DABS(YI),R0/EWT(I))
  150.  530      Y(I) = Y(I) + R
  151.         IERR = 0
  152.         CALL F (NEQ, TN, Y, FTEM, IERR)
  153.         IF (IERR .LT. 0) RETURN
  154.         DO 550 JJ = J,N,MBAND 
  155.           Y(JJ) = YH(JJ,1)
  156.           YJJ = Y(JJ)
  157.           R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ))
  158.           FAC = -HL0/R
  159.           I1 = MAX0(JJ-MU,1)
  160.           I2 = MIN0(JJ+ML,N)
  161.           II = JJ*MEB1 - ML + 2
  162.           DO 540 I = I1,I2
  163.  540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
  164.  550      CONTINUE
  165.  560    CONTINUE
  166.       NFE = NFE + MBA
  167. C ADD IDENTITY MATRIX. -------------------------------------------------
  168.  570  II = MBAND + 2
  169.       DO 580 I = 1,N
  170.         WM(II) = WM(II) + 1.0D0
  171.  580    II = II + MEBAND
  172. C DO LU DECOMPOSITION OF P. --------------------------------------------
  173.       CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
  174.       IF (IER .NE. 0) IERPJ = 1
  175.       RETURN
  176. C----------------------- END OF SUBROUTINE PREPJ -----------------------
  177.       END 
  178.