← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/sdntl.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/sdntl.f
 
1
*DECK SDNTL
 
2
      SUBROUTINE SDNTL (EPS, F, FA, HMAX, HOLD, IMPL, JTASK, MATDIM,
 
3
     8   MAXORD, MINT, MITER, ML, MU, N, NDE, SAVE1, T, UROUND, USERS,
 
4
     8   Y, YWT, H, MNTOLD, MTROLD, NFE, RC, YH, A, CONVRG, EL, FAC,
 
5
     8   IER, IPVT, NQ, NWAIT, RH, RMAX, SAVE2, TQ, TREND, ISWFLG,
 
6
     8   JSTATE)
 
7
C***BEGIN PROLOGUE  SDNTL
 
8
C***SUBSIDIARY
 
9
C***PURPOSE  Subroutine SDNTL is called to set parameters on the first
 
10
C            call to SDSTP, on an internal restart, or when the user has
 
11
C            altered MINT, MITER, and/or H.
 
12
C***LIBRARY   SLATEC (SDRIVE)
 
13
C***TYPE      SINGLE PRECISION (SDNTL-S, DDNTL-D, CDNTL-C)
 
14
C***AUTHOR  Kahaner, D. K., (NIST)
 
15
C             National Institute of Standards and Technology
 
16
C             Gaithersburg, MD  20899
 
17
C           Sutherland, C. D., (LANL)
 
18
C             Mail Stop D466
 
19
C             Los Alamos National Laboratory
 
20
C             Los Alamos, NM  87545
 
21
C***DESCRIPTION
 
22
C
 
23
C  On the first call, the order is set to 1 and the initial derivatives
 
24
C  are calculated.  RMAX is the maximum ratio by which H can be
 
25
C  increased in one step.  It is initially RMINIT to compensate
 
26
C  for the small initial H, but then is normally equal to RMNORM.
 
27
C  If a failure occurs (in corrector convergence or error test), RMAX
 
28
C  is set at RMFAIL for the next increase.
 
29
C  If the caller has changed MINT, or if JTASK = 0, SDCST is called
 
30
C  to set the coefficients of the method.  If the caller has changed H,
 
31
C  YH must be rescaled.  If H or MINT has been changed, NWAIT is
 
32
C  reset to NQ + 2 to prevent further increases in H for that many
 
33
C  steps.  Also, RC is reset.  RC is the ratio of new to old values of
 
34
C  the coefficient L(0)*H.  If the caller has changed MITER, RC is
 
35
C  set to 0 to force the partials to be updated, if partials are used.
 
36
C***ROUTINES CALLED  SDCST, SDSCL, SGBFA, SGBSL, SGEFA, SGESL, SNRM2
 
37
C***REVISION HISTORY  (YYMMDD)
 
38
C   790601  DATE WRITTEN
 
39
C   900329  Initial submission to SLATEC.
 
40
C***END PROLOGUE  SDNTL
 
41
      INTEGER I, IFLAG, IMPL, INFO, ISWFLG, JSTATE, JTASK, MATDIM,
 
42
     8        MAXORD, MINT, MITER, ML, MNTOLD, MTROLD, MU, N, NDE, NFE,
 
43
     8        NQ, NWAIT
 
44
      REAL A(MATDIM,*), EL(13,12), EPS, FAC(*), H, HMAX,
 
45
     8     HOLD, OLDL0, RC, RH, RMAX, RMINIT, SAVE1(*), SAVE2(*), SNRM2,
 
46
     8     SUM, T, TQ(3,12), TREND, UROUND, Y(*), YH(N,*), YWT(*)
 
47
      INTEGER IPVT(*)
 
48
      LOGICAL CONVRG, IER
 
49
      PARAMETER(RMINIT = 10000.E0)
 
50
C***FIRST EXECUTABLE STATEMENT  SDNTL
 
51
      IER = .FALSE.
 
52
      IF (JTASK .GE. 0) THEN
 
53
        IF (JTASK .EQ. 0) THEN
 
54
          CALL SDCST (MAXORD, MINT, ISWFLG,  EL, TQ)
 
55
          RMAX = RMINIT
 
56
        END IF
 
57
        RC = 0.E0
 
58
        CONVRG = .FALSE.
 
59
        TREND = 1.E0
 
60
        NQ = 1
 
61
        NWAIT = 3
 
62
        CALL F (N, T, Y, SAVE2)
 
63
        IF (N .EQ. 0) THEN
 
64
          JSTATE = 6
 
65
          RETURN
 
66
        END IF
 
67
        NFE = NFE + 1
 
68
        IF (IMPL .NE. 0) THEN
 
69
          IF (MITER .EQ. 3) THEN
 
70
            IFLAG = 0
 
71
            CALL USERS (Y, YH, YWT, SAVE1, SAVE2, T, H, EL, IMPL, N,
 
72
     8                  NDE, IFLAG)
 
73
            IF (IFLAG .EQ. -1) THEN
 
74
              IER = .TRUE.
 
75
              RETURN
 
76
            END IF
 
77
            IF (N .EQ. 0) THEN
 
78
              JSTATE = 10
 
79
              RETURN
 
80
            END IF
 
81
          ELSE IF (IMPL .EQ. 1) THEN
 
82
            IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
 
83
              CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
 
84
              IF (N .EQ. 0) THEN
 
85
                JSTATE = 9
 
86
                RETURN
 
87
              END IF
 
88
              CALL SGEFA (A, MATDIM, N, IPVT, INFO)
 
89
              IF (INFO .NE. 0) THEN
 
90
                IER = .TRUE.
 
91
                RETURN
 
92
              END IF
 
93
              CALL SGESL (A, MATDIM, N, IPVT, SAVE2, 0)
 
94
            ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
 
95
              CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
 
96
              IF (N .EQ. 0) THEN
 
97
                JSTATE = 9
 
98
                RETURN
 
99
              END IF
 
100
              CALL SGBFA (A, MATDIM, N, ML, MU, IPVT, INFO)
 
101
              IF (INFO .NE. 0) THEN
 
102
                IER = .TRUE.
 
103
                RETURN
 
104
              END IF
 
105
              CALL SGBSL (A, MATDIM, N, ML, MU, IPVT, SAVE2, 0)
 
106
            END IF
 
107
          ELSE IF (IMPL .EQ. 2) THEN
 
108
            CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
 
109
            IF (N .EQ. 0) THEN
 
110
              JSTATE = 9
 
111
              RETURN
 
112
            END IF
 
113
            DO 150 I = 1,NDE
 
114
              IF (A(I,1) .EQ. 0.E0) THEN
 
115
                IER = .TRUE.
 
116
                RETURN
 
117
              ELSE
 
118
                SAVE2(I) = SAVE2(I)/A(I,1)
 
119
              END IF
 
120
 150          CONTINUE
 
121
            DO 155 I = NDE+1,N
 
122
 155          A(I,1) = 0.E0
 
123
          ELSE IF (IMPL .EQ. 3) THEN
 
124
            IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
 
125
              CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
 
126
              IF (N .EQ. 0) THEN
 
127
                JSTATE = 9
 
128
                RETURN
 
129
              END IF
 
130
              CALL SGEFA (A, MATDIM, NDE, IPVT, INFO)
 
131
              IF (INFO .NE. 0) THEN
 
132
                IER = .TRUE.
 
133
                RETURN
 
134
              END IF
 
135
              CALL SGESL (A, MATDIM, NDE, IPVT, SAVE2, 0)
 
136
            ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
 
137
              CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
 
138
              IF (N .EQ. 0) THEN
 
139
                JSTATE = 9
 
140
                RETURN
 
141
              END IF
 
142
              CALL SGBFA (A, MATDIM, NDE, ML, MU, IPVT, INFO)
 
143
              IF (INFO .NE. 0) THEN
 
144
                IER = .TRUE.
 
145
                RETURN
 
146
              END IF
 
147
              CALL SGBSL (A, MATDIM, NDE, ML, MU, IPVT, SAVE2, 0)
 
148
            END IF
 
149
          END IF
 
150
        END IF
 
151
        DO 170 I = 1,NDE
 
152
 170      SAVE1(I) = SAVE2(I)/MAX(1.E0, YWT(I))
 
153
        SUM = SNRM2(NDE, SAVE1, 1)/SQRT(REAL(NDE))
 
154
        IF (SUM .GT. EPS/ABS(H)) H = SIGN(EPS/SUM, H)
 
155
        DO 180 I = 1,N
 
156
 180      YH(I,2) = H*SAVE2(I)
 
157
        IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. ISWFLG .EQ. 3) THEN
 
158
          DO 20 I = 1,N
 
159
 20         FAC(I) = SQRT(UROUND)
 
160
        END IF
 
161
      ELSE
 
162
        IF (MITER .NE. MTROLD) THEN
 
163
          MTROLD = MITER
 
164
          RC = 0.E0
 
165
          CONVRG = .FALSE.
 
166
        END IF
 
167
        IF (MINT .NE. MNTOLD) THEN
 
168
          MNTOLD = MINT
 
169
          OLDL0 = EL(1,NQ)
 
170
          CALL SDCST (MAXORD, MINT, ISWFLG,  EL, TQ)
 
171
          RC = RC*EL(1,NQ)/OLDL0
 
172
          NWAIT = NQ + 2
 
173
        END IF
 
174
        IF (H .NE. HOLD) THEN
 
175
          NWAIT = NQ + 2
 
176
          RH = H/HOLD
 
177
          CALL SDSCL (HMAX, N, NQ, RMAX,  HOLD, RC, RH, YH)
 
178
        END IF
 
179
      END IF
 
180
      RETURN
 
181
      END