Mercurial > hg > octave-nkf
annotate liboctave/cruft/odepack/intdy.f @ 18174:b72bcf5f78cc stable release-3-8-0
Version 3.8.0 released.
* configure.ac (OCTAVE_VERSION): Now 3.8.0.
(OCTAVE_RELEASE_DATE): Now 2013-12-27.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 27 Dec 2013 16:59:04 -0500 |
| parents | 648dabbb4c6b |
| children | 446c46af4b42 |
| rev | line source |
|---|---|
| 2329 | 1 SUBROUTINE INTDY (T, K, YH, NYH, DKY, IFLAG) |
| 2 CLLL. OPTIMIZE | |
| 3 INTEGER K, NYH, IFLAG | |
| 4 INTEGER IOWND, IOWNS, | |
| 5 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, | |
| 6 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU | |
| 7 INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1 | |
| 8 DOUBLE PRECISION T, YH, DKY | |
| 9 DOUBLE PRECISION ROWNS, | |
| 10 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND | |
| 11 DOUBLE PRECISION C, R, S, TP | |
|
7575
d20a2f261306
use * instead of 1 for assumed-size fortran arrays
John W. Eaton <jwe@octave.org>
parents:
4040
diff
changeset
|
12 DIMENSION YH(NYH,*), DKY(*) |
| 2329 | 13 COMMON /LS0001/ ROWNS(209), |
| 14 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, | |
| 15 3 IOWND(14), IOWNS(6), | |
| 16 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, | |
| 17 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU | |
| 18 C----------------------------------------------------------------------- | |
| 19 C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE | |
| 20 C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY. THIS ROUTINE | |
| 21 C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY | |
| 22 C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER. | |
| 23 C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.) | |
| 24 C----------------------------------------------------------------------- | |
| 25 C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE | |
| 26 C NORDSIECK HISTORY ARRAY YH. THIS ARRAY CORRESPONDS UNIQUELY TO A | |
| 27 C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET | |
| 28 C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T. | |
| 29 C THE FORMULA FOR DKY IS.. | |
| 30 C Q | |
| 31 C DKY(I) = SUM C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1) | |
| 32 C J=K | |
| 33 C WHERE C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR. | |
| 34 C THE QUANTITIES NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE | |
| 35 C COMMUNICATED BY COMMON. THE ABOVE SUM IS DONE IN REVERSE ORDER. | |
| 36 C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS. | |
| 37 C----------------------------------------------------------------------- | |
| 38 IFLAG = 0 | |
| 39 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 | |
| 40 TP = TN - HU - 100.0D0*UROUND*(TN + HU) | |
| 41 IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90 | |
| 42 C | |
| 43 S = (T - TN)/H | |
| 44 IC = 1 | |
| 45 IF (K .EQ. 0) GO TO 15 | |
| 46 JJ1 = L - K | |
| 47 DO 10 JJ = JJ1,NQ | |
| 48 10 IC = IC*JJ | |
| 49 15 C = DBLE(IC) | |
| 50 DO 20 I = 1,N | |
| 51 20 DKY(I) = C*YH(I,L) | |
| 52 IF (K .EQ. NQ) GO TO 55 | |
| 53 JB2 = NQ - K | |
| 54 DO 50 JB = 1,JB2 | |
| 55 J = NQ - JB | |
| 56 JP1 = J + 1 | |
| 57 IC = 1 | |
| 58 IF (K .EQ. 0) GO TO 35 | |
| 59 JJ1 = JP1 - K | |
| 60 DO 30 JJ = JJ1,J | |
| 61 30 IC = IC*JJ | |
| 62 35 C = DBLE(IC) | |
| 63 DO 40 I = 1,N | |
| 64 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) | |
| 65 50 CONTINUE | |
| 66 IF (K .EQ. 0) RETURN | |
| 67 55 R = H**(-K) | |
| 68 DO 60 I = 1,N | |
| 69 60 DKY(I) = R*DKY(I) | |
| 70 RETURN | |
| 71 C | |
| 4040 | 72 80 CALL XERRWD('INTDY-- K (=I1) ILLEGAL ', |
| 2329 | 73 1 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0) |
| 74 IFLAG = -1 | |
| 75 RETURN | |
| 4040 | 76 90 CALL XERRWD('INTDY-- T (=R1) ILLEGAL ', |
| 2329 | 77 1 30, 52, 0, 0, 0, 0, 1, T, 0.0D0) |
| 4040 | 78 CALL XERRWD( |
| 79 1 ' T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2) ', | |
| 2329 | 80 1 60, 52, 0, 0, 0, 0, 2, TP, TN) |
| 81 IFLAG = -2 | |
| 82 RETURN | |
| 83 C----------------------- END OF SUBROUTINE INTDY ----------------------- | |
| 84 END |