← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dpjac.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/dpjac.f
 
1
*DECK DPJAC
 
2
      SUBROUTINE DPJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, DF,
 
3
     +   DJAC, RPAR, IPAR)
 
4
C***BEGIN PROLOGUE  DPJAC
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to DDEBDF
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      DOUBLE PRECISION (PJAC-S, DPJAC-D)
 
9
C***AUTHOR  Watts, H. A., (SNLA)
 
10
C***DESCRIPTION
 
11
C
 
12
C   DPJAC sets up the iteration matrix (involving the Jacobian) for the
 
13
C   integration package DDEBDF.
 
14
C
 
15
C***SEE ALSO  DDEBDF
 
16
C***ROUTINES CALLED  DGBFA, DGEFA, DVNRMS
 
17
C***COMMON BLOCKS    DDEBD1
 
18
C***REVISION HISTORY  (YYMMDD)
 
19
C   820301  DATE WRITTEN
 
20
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
21
C   890911  Removed unnecessary intrinsics.  (WRB)
 
22
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
23
C   900328  Added TYPE section.  (WRB)
 
24
C   910722  Updated AUTHOR section.  (ALS)
 
25
C   920422  Changed DIMENSION statement.  (WRB)
 
26
C***END PROLOGUE  DPJAC
 
27
C
 
28
      INTEGER I, I1, I2, IER, II, IOWND, IOWNS, IPAR, IWM, J, J1,
 
29
     1      JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND,
 
30
     2      MEB1, MEBAND, METH, MITER, ML, ML3, MU, N, NEQ,
 
31
     3      NFE, NJE, NQ, NQU, NST, NYH
 
32
      DOUBLE PRECISION CON, DI, DVNRMS, EL0, EWT,
 
33
     1      FAC, FTEM, H, HL0, HMIN, HMXI, HU, R, R0, ROWND, ROWNS,
 
34
     2      RPAR, SAVF, SRUR, TN, UROUND, WM, Y, YH, YI, YJ, YJJ
 
35
      EXTERNAL DF, DJAC
 
36
      DIMENSION Y(*),YH(NYH,*),EWT(*),FTEM(*),SAVF(*),WM(*),IWM(*),
 
37
     1          RPAR(*),IPAR(*)
 
38
      COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND,
 
39
     1                IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER,
 
40
     2                MAXORD,N,NQ,NST,NFE,NJE,NQU
 
41
C     ------------------------------------------------------------------
 
42
C      DPJAC IS CALLED BY DSTOD  TO COMPUTE AND PROCESS THE MATRIX
 
43
C      P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
 
44
C      HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE DJAC IF
 
45
C      MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
 
46
C      IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
 
47
C      J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
 
48
C      SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
 
49
C      OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
 
50
C      BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
 
51
C
 
52
C      IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
 
53
C      WITH DPJAC USES THE FOLLOWING..
 
54
C      Y    = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
 
55
C      FTEM = WORK ARRAY OF LENGTH N (ACOR IN DSTOD ).
 
56
C      SAVF = ARRAY CONTAINING DF EVALUATED AT PREDICTED Y.
 
57
C      WM   = DOUBLE PRECISION WORK SPACE FOR MATRICES.  ON OUTPUT IT
 
58
C      CONTAINS THE
 
59
C             INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU
 
60
C             DECOMPOSITION OF P IF MITER IS 1, 2 , 4, OR 5.
 
61
C             STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
 
62
C             WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
 
63
C             WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN
 
64
C             INCREMENTS.  WM(2) = H*EL0, SAVED FOR LATER USE IF MITER =
 
65
C             3.
 
66
C      IWM  = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
 
67
C             AT IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS
 
68
C             THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER
 
69
C             IS 4 OR 5.
 
70
C      EL0  = EL(1) (INPUT).
 
71
C      IER  = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .NE. 0 IF
 
72
C             P MATRIX FOUND TO BE SINGULAR.
 
73
C      THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
 
74
C      MITER, N, NFE, AND NJE.
 
75
C-----------------------------------------------------------------------
 
76
C     BEGIN BLOCK PERMITTING ...EXITS TO 240
 
77
C        BEGIN BLOCK PERMITTING ...EXITS TO 220
 
78
C           BEGIN BLOCK PERMITTING ...EXITS TO 130
 
79
C              BEGIN BLOCK PERMITTING ...EXITS TO 70
 
80
C***FIRST EXECUTABLE STATEMENT  DPJAC
 
81
                  NJE = NJE + 1
 
82
                  HL0 = H*EL0
 
83
                  GO TO (10,40,90,140,170), MITER
 
84
C                 IF MITER = 1, CALL DJAC AND MULTIPLY BY SCALAR.
 
85
C                 -----------------------
 
86
   10             CONTINUE
 
87
                  LENP = N*N
 
88
                  DO 20 I = 1, LENP
 
89
                     WM(I+2) = 0.0D0
 
90
   20             CONTINUE
 
91
                  CALL DJAC(TN,Y,WM(3),N,RPAR,IPAR)
 
92
                  CON = -HL0
 
93
                  DO 30 I = 1, LENP
 
94
                     WM(I+2) = WM(I+2)*CON
 
95
   30             CONTINUE
 
96
C              ...EXIT
 
97
                  GO TO 70
 
98
C                 IF MITER = 2, MAKE N CALLS TO DF TO APPROXIMATE J.
 
99
C                 --------------------
 
100
   40             CONTINUE
 
101
                  FAC = DVNRMS(N,SAVF,EWT)
 
102
                  R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
 
103
                  IF (R0 .EQ. 0.0D0) R0 = 1.0D0
 
104
                  SRUR = WM(1)
 
105
                  J1 = 2
 
106
                  DO 60 J = 1, N
 
107
                     YJ = Y(J)
 
108
                     R = MAX(SRUR*ABS(YJ),R0*EWT(J))
 
109
                     Y(J) = Y(J) + R
 
110
                     FAC = -HL0/R
 
111
                     CALL DF(TN,Y,FTEM,RPAR,IPAR)
 
112
                     DO 50 I = 1, N
 
113
                        WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
 
114
   50                CONTINUE
 
115
                     Y(J) = YJ
 
116
                     J1 = J1 + N
 
117
   60             CONTINUE
 
118
                  NFE = NFE + N
 
119
   70          CONTINUE
 
120
C              ADD IDENTITY MATRIX.
 
121
C              -------------------------------------------------
 
122
               J = 3
 
123
               DO 80 I = 1, N
 
124
                  WM(J) = WM(J) + 1.0D0
 
125
                  J = J + (N + 1)
 
126
   80          CONTINUE
 
127
C              DO LU DECOMPOSITION ON P.
 
128
C              --------------------------------------------
 
129
               CALL DGEFA(WM(3),N,N,IWM(21),IER)
 
130
C     .........EXIT
 
131
               GO TO 240
 
132
C              IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND
 
133
C              P. ---------
 
134
   90          CONTINUE
 
135
               WM(2) = HL0
 
136
               IER = 0
 
137
               R = EL0*0.1D0
 
138
               DO 100 I = 1, N
 
139
                  Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
 
140
  100          CONTINUE
 
141
               CALL DF(TN,Y,WM(3),RPAR,IPAR)
 
142
               NFE = NFE + 1
 
143
               DO 120 I = 1, N
 
144
                  R0 = H*SAVF(I) - YH(I,2)
 
145
                  DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
 
146
                  WM(I+2) = 1.0D0
 
147
                  IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 110
 
148
C           .........EXIT
 
149
                     IF (ABS(DI) .EQ. 0.0D0) GO TO 130
 
150
                     WM(I+2) = 0.1D0*R0/DI
 
151
  110             CONTINUE
 
152
  120          CONTINUE
 
153
C     .........EXIT
 
154
               GO TO 240
 
155
  130       CONTINUE
 
156
            IER = -1
 
157
C     ......EXIT
 
158
            GO TO 240
 
159
C           IF MITER = 4, CALL DJAC AND MULTIPLY BY SCALAR.
 
160
C           -----------------------
 
161
  140       CONTINUE
 
162
            ML = IWM(1)
 
163
            MU = IWM(2)
 
164
            ML3 = 3
 
165
            MBAND = ML + MU + 1
 
166
            MEBAND = MBAND + ML
 
167
            LENP = MEBAND*N
 
168
            DO 150 I = 1, LENP
 
169
               WM(I+2) = 0.0D0
 
170
  150       CONTINUE
 
171
            CALL DJAC(TN,Y,WM(ML3),MEBAND,RPAR,IPAR)
 
172
            CON = -HL0
 
173
            DO 160 I = 1, LENP
 
174
               WM(I+2) = WM(I+2)*CON
 
175
  160       CONTINUE
 
176
C        ...EXIT
 
177
            GO TO 220
 
178
C           IF MITER = 5, MAKE MBAND CALLS TO DF TO APPROXIMATE J.
 
179
C           ----------------
 
180
  170       CONTINUE
 
181
            ML = IWM(1)
 
182
            MU = IWM(2)
 
183
            MBAND = ML + MU + 1
 
184
            MBA = MIN(MBAND,N)
 
185
            MEBAND = MBAND + ML
 
186
            MEB1 = MEBAND - 1
 
187
            SRUR = WM(1)
 
188
            FAC = DVNRMS(N,SAVF,EWT)
 
189
            R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
 
190
            IF (R0 .EQ. 0.0D0) R0 = 1.0D0
 
191
            DO 210 J = 1, MBA
 
192
               DO 180 I = J, N, MBAND
 
193
                  YI = Y(I)
 
194
                  R = MAX(SRUR*ABS(YI),R0*EWT(I))
 
195
                  Y(I) = Y(I) + R
 
196
  180          CONTINUE
 
197
               CALL DF(TN,Y,FTEM,RPAR,IPAR)
 
198
               DO 200 JJ = J, N, MBAND
 
199
                  Y(JJ) = YH(JJ,1)
 
200
                  YJJ = Y(JJ)
 
201
                  R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ))
 
202
                  FAC = -HL0/R
 
203
                  I1 = MAX(JJ-MU,1)
 
204
                  I2 = MIN(JJ+ML,N)
 
205
                  II = JJ*MEB1 - ML + 2
 
206
                  DO 190 I = I1, I2
 
207
                     WM(II+I) = (FTEM(I) - SAVF(I))*FAC
 
208
  190             CONTINUE
 
209
  200          CONTINUE
 
210
  210       CONTINUE
 
211
            NFE = NFE + MBA
 
212
  220    CONTINUE
 
213
C        ADD IDENTITY MATRIX.
 
214
C        -------------------------------------------------
 
215
         II = MBAND + 2
 
216
         DO 230 I = 1, N
 
217
            WM(II) = WM(II) + 1.0D0
 
218
            II = II + MEBAND
 
219
  230    CONTINUE
 
220
C        DO LU DECOMPOSITION OF P.
 
221
C        --------------------------------------------
 
222
         CALL DGBFA(WM(3),MEBAND,N,ML,MU,IWM(21),IER)
 
223
  240 CONTINUE
 
224
      RETURN
 
225
C     ----------------------- END OF SUBROUTINE DPJAC
 
226
C     -----------------------
 
227
      END