← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/pjac.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/pjac.f
 
1
*DECK PJAC
 
2
      SUBROUTINE PJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F,
 
3
     +   JAC, RPAR, IPAR)
 
4
C***BEGIN PROLOGUE  PJAC
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to DEBDF
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      SINGLE PRECISION (PJAC-S, DPJAC-D)
 
9
C***AUTHOR  Watts, H. A., (SNLA)
 
10
C***DESCRIPTION
 
11
C
 
12
C   PJAC sets up the iteration matrix (involving the Jacobian) for the
 
13
C   integration package DEBDF.
 
14
C
 
15
C***SEE ALSO  DEBDF
 
16
C***ROUTINES CALLED  SGBFA, SGEFA, VNWRMS
 
17
C***COMMON BLOCKS    DEBDF1
 
18
C***REVISION HISTORY  (YYMMDD)
 
19
C   800901  DATE WRITTEN
 
20
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
21
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
22
C   900328  Added TYPE section.  (WRB)
 
23
C   910722  Updated AUTHOR section.  (ALS)
 
24
C   920422  Changed DIMENSION statement.  (WRB)
 
25
C***END PROLOGUE  PJAC
 
26
C
 
27
CLLL. OPTIMIZE
 
28
      INTEGER NEQ, NYH, IWM, I, I1, I2, IER, II, IOWND, IOWNS, J, J1,
 
29
     1   JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND, MEB1, MEBAND,
 
30
     2   METH, MITER, ML, ML3, MU, N, NFE, NJE, NQ, NQU, NST
 
31
      EXTERNAL F, JAC
 
32
      REAL Y, YH, EWT, FTEM, SAVF, WM,
 
33
     1   ROWND, ROWNS, EL0, H, HMIN, HMXI, HU, TN, UROUND,
 
34
     2   CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, VNWRMS
 
35
      DIMENSION         Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
 
36
     1   WM(*), IWM(*), RPAR(*), IPAR(*)
 
37
      COMMON /DEBDF1/ ROWND, ROWNS(210),
 
38
     1   EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(14), IOWNS(6),
 
39
     2   IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE,
 
40
     3   NJE, NQU
 
41
C-----------------------------------------------------------------------
 
42
C PJAC IS CALLED BY STOD  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 JAC 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 SGEFA IF MITER = 1 OR 2, AND BY SGBFA IF MITER = 4 OR 5.
 
51
C
 
52
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
 
53
C WITH PJAC USES THE FOLLOWING..
 
54
C Y    = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
 
55
C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STOD ).
 
56
C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
 
57
C WM   = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
 
58
C        INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
 
59
C        OF P IF MITER IS 1, 2 , 4, OR 5.
 
60
C        STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
 
61
C        WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
 
62
C        WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
 
63
C        WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
 
64
C IWM  = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
 
65
C        IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS THE
 
66
C        BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
 
67
C EL0  = EL(1) (INPUT).
 
68
C IER  = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .NE. 0 IF
 
69
C        P MATRIX FOUND TO BE SINGULAR.
 
70
C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
 
71
C MITER, N, NFE, AND NJE.
 
72
C-----------------------------------------------------------------------
 
73
C***FIRST EXECUTABLE STATEMENT  PJAC
 
74
      NJE = NJE + 1
 
75
      HL0 = H*EL0
 
76
      GO TO (100, 200, 300, 400, 500), MITER
 
77
C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 
78
 100  LENP = N*N
 
79
      DO 110 I = 1,LENP
 
80
 110    WM(I+2) = 0.0E0
 
81
      CALL JAC (TN, Y, WM(3), N, RPAR, IPAR)
 
82
      CON = -HL0
 
83
      DO 120 I = 1,LENP
 
84
 120    WM(I+2) = WM(I+2)*CON
 
85
      GO TO 240
 
86
C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
 
87
 200  FAC = VNWRMS (N, SAVF, EWT)
 
88
      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
 
89
      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
 
90
      SRUR = WM(1)
 
91
      J1 = 2
 
92
      DO 230 J = 1,N
 
93
        YJ = Y(J)
 
94
        R = MAX(SRUR*ABS(YJ),R0*EWT(J))
 
95
        Y(J) = Y(J) + R
 
96
        FAC = -HL0/R
 
97
        CALL F (TN, Y, FTEM, RPAR, IPAR)
 
98
        DO 220 I = 1,N
 
99
 220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
 
100
        Y(J) = YJ
 
101
        J1 = J1 + N
 
102
 230    CONTINUE
 
103
      NFE = NFE + N
 
104
C ADD IDENTITY MATRIX. -------------------------------------------------
 
105
 240  J = 3
 
106
      DO 250 I = 1,N
 
107
        WM(J) = WM(J) + 1.0E0
 
108
 250    J = J + (N + 1)
 
109
C DO LU DECOMPOSITION ON P. --------------------------------------------
 
110
      CALL SGEFA (WM(3), N, N, IWM(21), IER)
 
111
      RETURN
 
112
C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
 
113
 300  WM(2) = HL0
 
114
      IER = 0
 
115
      R = EL0*0.1E0
 
116
      DO 310 I = 1,N
 
117
 310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
 
118
      CALL F (TN, Y, WM(3), RPAR, IPAR)
 
119
      NFE = NFE + 1
 
120
      DO 320 I = 1,N
 
121
        R0 = H*SAVF(I) - YH(I,2)
 
122
        DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I))
 
123
        WM(I+2) = 1.0E0
 
124
        IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 320
 
125
        IF (ABS(DI) .EQ. 0.0E0) GO TO 330
 
126
        WM(I+2) = 0.1E0*R0/DI
 
127
 320    CONTINUE
 
128
      RETURN
 
129
 330  IER = -1
 
130
      RETURN
 
131
C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 
132
 400  ML = IWM(1)
 
133
      MU = IWM(2)
 
134
      ML3 =  3
 
135
      MBAND = ML + MU + 1
 
136
      MEBAND = MBAND + ML
 
137
      LENP = MEBAND*N
 
138
      DO 410 I = 1,LENP
 
139
 410    WM(I+2) = 0.0E0
 
140
      CALL JAC (TN, Y, WM(ML3), MEBAND, RPAR, IPAR)
 
141
      CON = -HL0
 
142
      DO 420 I = 1,LENP
 
143
 420    WM(I+2) = WM(I+2)*CON
 
144
      GO TO 570
 
145
C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
 
146
 500  ML = IWM(1)
 
147
      MU = IWM(2)
 
148
      MBAND = ML + MU + 1
 
149
      MBA = MIN(MBAND,N)
 
150
      MEBAND = MBAND + ML
 
151
      MEB1 = MEBAND - 1
 
152
      SRUR = WM(1)
 
153
      FAC = VNWRMS (N, SAVF, EWT)
 
154
      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
 
155
      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
 
156
      DO 560 J = 1,MBA
 
157
        DO 530 I = J,N,MBAND
 
158
          YI = Y(I)
 
159
          R = MAX(SRUR*ABS(YI),R0*EWT(I))
 
160
 530      Y(I) = Y(I) + R
 
161
        CALL F (TN, Y, FTEM, RPAR, IPAR)
 
162
        DO 550 JJ = J,N,MBAND
 
163
          Y(JJ) = YH(JJ,1)
 
164
          YJJ = Y(JJ)
 
165
          R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ))
 
166
          FAC = -HL0/R
 
167
          I1 = MAX(JJ-MU,1)
 
168
          I2 = MIN(JJ+ML,N)
 
169
          II = JJ*MEB1 - ML + 2
 
170
          DO 540 I = I1,I2
 
171
 540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
 
172
 550      CONTINUE
 
173
 560    CONTINUE
 
174
      NFE = NFE + MBA
 
175
C ADD IDENTITY MATRIX. -------------------------------------------------
 
176
 570  II = MBAND + 2
 
177
      DO 580 I = 1,N
 
178
        WM(II) = WM(II) + 1.0E0
 
179
 580    II = II + MEBAND
 
180
C DO LU DECOMPOSITION OF P. --------------------------------------------
 
181
      CALL SGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
 
182
      RETURN
 
183
C----------------------- END OF SUBROUTINE PJAC -----------------------
 
184
      END