← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/integrate/odepack/stode.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/integrate/odepack/stode.f
1
 
      subroutine stode (neq, y, yh, nyh, yh1, ewt, savf, acor,
2
 
     1   wm, iwm, f, jac, pjac, slvs)
3
 
clll. optimize
4
 
      external f, jac, pjac, slvs
5
 
      integer neq, nyh, iwm
6
 
      integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
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, iredo, iret, j, jb, m, ncf, newq
10
 
      double precision y, yh, yh1, ewt, savf, acor, wm
11
 
      double precision conit, crate, el, elco, hold, rmax, tesco,
12
 
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
13
 
      double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
14
 
     1   r, rh, rhdn, rhsm, rhup, told, vnorm
15
 
      dimension neq(1), y(1), yh(nyh,*), yh1(1), ewt(1), savf(1),
16
 
     1   acor(1), wm(*), iwm(*)
17
 
      common /ls0001/ conit, crate, el(13), elco(13,12),
18
 
     1   hold, rmax, tesco(3,12),
19
 
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14),
20
 
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
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 stode performs one step of the integration of an initial value
25
 
c problem for a system of ordinary differential equations.
26
 
c note.. stode is independent of the value of the iteration method
27
 
c indicator miter, when this is .ne. 0, and hence is independent
28
 
c of the type of chord method used, or the jacobian structure.
29
 
c communication with stode is done with the following variables..
30
 
c
31
 
c neq    = integer array containing problem size in neq(1), and
32
 
c          passed as the neq argument in all calls to f and jac.
33
 
c y      = an array of length .ge. n used as the y argument in
34
 
c          all calls to f and jac.
35
 
c yh     = an nyh by lmax array containing the dependent variables
36
 
c          and their approximate scaled derivatives, where
37
 
c          lmax = maxord + 1.  yh(i,j+1) contains the approximate
38
 
c          j-th derivative of y(i), scaled by h**j/factorial(j)
39
 
c          (j = 0,1,...,nq).  on entry for the first step, the first
40
 
c          two columns of yh must be set from the initial values.
41
 
c nyh    = a constant integer .ge. n, the first dimension of yh.
42
 
c yh1    = a one-dimensional array occupying the same space as yh.
43
 
c ewt    = an array of length n containing multiplicative weights
44
 
c          for local error measurements.  local errors in y(i) are
45
 
c          compared to 1.0/ewt(i) in various error tests.
46
 
c savf   = an array of working storage, of length n.
47
 
c          also used for input of yh(*,maxord+2) when jstart = -1
48
 
c          and maxord .lt. the current order nq.
49
 
c acor   = a work array of length n, used for the accumulated
50
 
c          corrections.  on a successful return, acor(i) contains
51
 
c          the estimated one-step local error in y(i).
52
 
c wm,iwm = real and integer work arrays associated with matrix
53
 
c          operations in chord iteration (miter .ne. 0).
54
 
c pjac   = name of routine to evaluate and preprocess jacobian matrix
55
 
c          and p = i - h*el0*jac, if a chord method is being used.
56
 
c slvs   = name of routine to solve linear system in chord iteration.
57
 
c ccmax  = maximum relative change in h*el0 before pjac is called.
58
 
c h      = the step size to be attempted on the next step.
59
 
c          h is altered by the error control algorithm during the
60
 
c          problem.  h can be either positive or negative, but its
61
 
c          sign must remain constant throughout the problem.
62
 
c hmin   = the minimum absolute value of the step size h to be used.
63
 
c hmxi   = inverse of the maximum absolute value of h to be used.
64
 
c          hmxi = 0.0 is allowed and corresponds to an infinite hmax.
65
 
c          hmin and hmxi may be changed at any time, but will not
66
 
c          take effect until the next change of h is considered.
67
 
c tn     = the independent variable. tn is updated on each step taken.
68
 
c jstart = an integer used for input only, with the following
69
 
c          values and meanings..
70
 
c               0  perform the first step.
71
 
c           .gt.0  take a new step continuing from the last.
72
 
c              -1  take the next step with a new value of h, maxord,
73
 
c                    n, meth, miter, and/or matrix parameters.
74
 
c              -2  take the next step with a new value of h,
75
 
c                    but with other inputs unchanged.
76
 
c          on return, jstart is set to 1 to facilitate continuation.
77
 
c kflag  = a completion code with the following meanings..
78
 
c               0  the step was succesful.
79
 
c              -1  the requested error could not be achieved.
80
 
c              -2  corrector convergence could not be achieved.
81
 
c              -3  fatal error in pjac or slvs.
82
 
c          a return with kflag = -1 or -2 means either
83
 
c          abs(h) = hmin or 10 consecutive failures occurred.
84
 
c          on a return with kflag negative, the values of tn and
85
 
c          the yh array are as of the beginning of the last
86
 
c          step, and h is the last step size attempted.
87
 
c maxord = the maximum order of integration method to be allowed.
88
 
c maxcor = the maximum number of corrector iterations allowed.
89
 
c msbp   = maximum number of steps between pjac calls (miter .gt. 0).
90
 
c mxncf  = maximum number of convergence failures allowed.
91
 
c meth/miter = the method flags.  see description in driver.
92
 
c n      = the number of first-order differential equations.
93
 
c-----------------------------------------------------------------------
94
 
      kflag = 0
95
 
      told = tn
96
 
      ncf = 0
97
 
      ierpj = 0
98
 
      iersl = 0
99
 
      jcur = 0
100
 
      icf = 0
101
 
      delp = 0.0d0
102
 
      if (jstart .gt. 0) go to 200
103
 
      if (jstart .eq. -1) go to 100
104
 
      if (jstart .eq. -2) go to 160
105
 
c-----------------------------------------------------------------------
106
 
c on the first call, the order is set to 1, and other variables are
107
 
c initialized.  rmax is the maximum ratio by which h can be increased
108
 
c in a single step.  it is initially 1.e4 to compensate for the small
109
 
c initial h, but then is normally equal to 10.  if a failure
110
 
c occurs (in corrector convergence or error test), rmax is set at 2
111
 
c for the next increase.
112
 
c-----------------------------------------------------------------------
113
 
      lmax = maxord + 1
114
 
      nq = 1
115
 
      l = 2
116
 
      ialth = 2
117
 
      rmax = 10000.0d0
118
 
      rc = 0.0d0
119
 
      el0 = 1.0d0
120
 
      crate = 0.7d0
121
 
      hold = h
122
 
      meo = meth
123
 
      nslp = 0
124
 
      ipup = miter
125
 
      iret = 3
126
 
      go to 140
127
 
c-----------------------------------------------------------------------
128
 
c the following block handles preliminaries needed when jstart = -1.
129
 
c ipup is set to miter to force a matrix update.
130
 
c if an order increase is about to be considered (ialth = 1),
131
 
c ialth is reset to 2 to postpone consideration one more step.
132
 
c if the caller has changed meth, cfode is called to reset
133
 
c the coefficients of the method.
134
 
c if the caller has changed maxord to a value less than the current
135
 
c order nq, nq is reduced to maxord, and a new h chosen accordingly.
136
 
c if h is to be changed, yh must be rescaled.
137
 
c if h or meth is being changed, ialth is reset to l = nq + 1
138
 
c to prevent further changes in h for that many steps.
139
 
c-----------------------------------------------------------------------
140
 
 100  ipup = miter
141
 
      lmax = maxord + 1
142
 
      if (ialth .eq. 1) ialth = 2
143
 
      if (meth .eq. meo) go to 110
144
 
      call cfode (meth, elco, tesco)
145
 
      meo = meth
146
 
      if (nq .gt. maxord) go to 120
147
 
      ialth = l
148
 
      iret = 1
149
 
      go to 150
150
 
 110  if (nq .le. maxord) go to 160
151
 
 120  nq = maxord
152
 
      l = lmax
153
 
      do 125 i = 1,l
154
 
 125    el(i) = elco(i,nq)
155
 
      nqnyh = nq*nyh
156
 
      rc = rc*el(1)/el0
157
 
      el0 = el(1)
158
 
      conit = 0.5d0/dfloat(nq+2)
159
 
      ddn = vnorm (n, savf, ewt)/tesco(1,l)
160
 
      exdn = 1.0d0/dfloat(l)
161
 
      rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
162
 
      rh = dmin1(rhdn,1.0d0)
163
 
      iredo = 3
164
 
      if (h .eq. hold) go to 170
165
 
      rh = dmin1(rh,dabs(h/hold))
166
 
      h = hold
167
 
      go to 175
168
 
c-----------------------------------------------------------------------
169
 
c cfode is called to get all the integration coefficients for the
170
 
c current meth.  then the el vector and related constants are reset
171
 
c whenever the order nq is changed, or at the start of the problem.
172
 
c-----------------------------------------------------------------------
173
 
 140  call cfode (meth, elco, tesco)
174
 
 150  do 155 i = 1,l
175
 
 155    el(i) = elco(i,nq)
176
 
      nqnyh = nq*nyh
177
 
      rc = rc*el(1)/el0
178
 
      el0 = el(1)
179
 
      conit = 0.5d0/dfloat(nq+2)
180
 
      go to (160, 170, 200), iret
181
 
c-----------------------------------------------------------------------
182
 
c if h is being changed, the h ratio rh is checked against
183
 
c rmax, hmin, and hmxi, and the yh array rescaled.  ialth is set to
184
 
c l = nq + 1 to prevent a change of h for that many steps, unless
185
 
c forced by a convergence or error test failure.
186
 
c-----------------------------------------------------------------------
187
 
 160  if (h .eq. hold) go to 200
188
 
      rh = h/hold
189
 
      h = hold
190
 
      iredo = 3
191
 
      go to 175
192
 
 170  rh = dmax1(rh,hmin/dabs(h))
193
 
 175  rh = dmin1(rh,rmax)
194
 
      rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh)
195
 
      r = 1.0d0
196
 
      do 180 j = 2,l
197
 
        r = r*rh
198
 
        do 180 i = 1,n
199
 
 180      yh(i,j) = yh(i,j)*r
200
 
      h = h*rh
201
 
      rc = rc*rh
202
 
      ialth = l
203
 
      if (iredo .eq. 0) go to 690
204
 
c-----------------------------------------------------------------------
205
 
c this section computes the predicted values by effectively
206
 
c multiplying the yh array by the pascal triangle matrix.
207
 
c rc is the ratio of new to old values of the coefficient  h*el(1).
208
 
c when rc differs from 1 by more than ccmax, ipup is set to miter
209
 
c to force pjac to be called, if a jacobian is involved.
210
 
c in any case, pjac is called at least every msbp steps.
211
 
c-----------------------------------------------------------------------
212
 
 200  if (dabs(rc-1.0d0) .gt. ccmax) ipup = miter
213
 
      if (nst .ge. nslp+msbp) ipup = miter
214
 
      tn = tn + h
215
 
      i1 = nqnyh + 1
216
 
      do 215 jb = 1,nq
217
 
        i1 = i1 - nyh
218
 
cdir$ ivdep
219
 
        do 210 i = i1,nqnyh
220
 
 210      yh1(i) = yh1(i) + yh1(i+nyh)
221
 
 215    continue
222
 
c-----------------------------------------------------------------------
223
 
c up to maxcor corrector iterations are taken.  a convergence test is
224
 
c made on the r.m.s. norm of each correction, weighted by the error
225
 
c weight vector ewt.  the sum of the corrections is accumulated in the
226
 
c vector acor(i).  the yh array is not altered in the corrector loop.
227
 
c-----------------------------------------------------------------------
228
 
 220  m = 0
229
 
      do 230 i = 1,n
230
 
 230    y(i) = yh(i,1)
231
 
      call f (neq, tn, y, savf)
232
 
      nfe = nfe + 1
233
 
      if (ipup .le. 0) go to 250
234
 
c-----------------------------------------------------------------------
235
 
c if indicated, the matrix p = i - h*el(1)*j is reevaluated and
236
 
c preprocessed before starting the corrector iteration.  ipup is set
237
 
c to 0 as an indicator that this has been done.
238
 
c-----------------------------------------------------------------------
239
 
      call pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
240
 
      ipup = 0
241
 
      rc = 1.0d0
242
 
      nslp = nst
243
 
      crate = 0.7d0
244
 
      if (ierpj .ne. 0) go to 430
245
 
 250  do 260 i = 1,n
246
 
 260    acor(i) = 0.0d0
247
 
 270  if (miter .ne. 0) go to 350
248
 
c-----------------------------------------------------------------------
249
 
c in the case of functional iteration, update y directly from
250
 
c the result of the last function evaluation.
251
 
c-----------------------------------------------------------------------
252
 
      do 290 i = 1,n
253
 
        savf(i) = h*savf(i) - yh(i,2)
254
 
 290    y(i) = savf(i) - acor(i)
255
 
      del = vnorm (n, y, ewt)
256
 
      do 300 i = 1,n
257
 
        y(i) = yh(i,1) + el(1)*savf(i)
258
 
 300    acor(i) = savf(i)
259
 
      go to 400
260
 
c-----------------------------------------------------------------------
261
 
c in the case of the chord method, compute the corrector error,
262
 
c and solve the linear system with that as right-hand side and
263
 
c p as coefficient matrix.
264
 
c-----------------------------------------------------------------------
265
 
 350  do 360 i = 1,n
266
 
 360    y(i) = h*savf(i) - (yh(i,2) + acor(i))
267
 
      call slvs (wm, iwm, y, savf)
268
 
      if (iersl .lt. 0) go to 430
269
 
      if (iersl .gt. 0) go to 410
270
 
      del = vnorm (n, y, ewt)
271
 
      do 380 i = 1,n
272
 
        acor(i) = acor(i) + y(i)
273
 
 380    y(i) = yh(i,1) + el(1)*acor(i)
274
 
c-----------------------------------------------------------------------
275
 
c test for convergence.  if m.gt.0, an estimate of the convergence
276
 
c rate constant is stored in crate, and this is used in the test.
277
 
c-----------------------------------------------------------------------
278
 
 400  if (m .ne. 0) crate = dmax1(0.2d0*crate,del/delp)
279
 
      dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
280
 
      if (dcon .le. 1.0d0) go to 450
281
 
      m = m + 1
282
 
      if (m .eq. maxcor) go to 410
283
 
      if (m .ge. 2 .and. del .gt. 2.0d0*delp) go to 410
284
 
      delp = del
285
 
      call f (neq, tn, y, savf)
286
 
      nfe = nfe + 1
287
 
      go to 270
288
 
c-----------------------------------------------------------------------
289
 
c the corrector iteration failed to converge.
290
 
c if miter .ne. 0 and the jacobian is out of date, pjac is called for
291
 
c the next try.  otherwise the yh array is retracted to its values
292
 
c before prediction, and h is reduced, if possible.  if h cannot be
293
 
c reduced or mxncf failures have occurred, exit with kflag = -2.
294
 
c-----------------------------------------------------------------------
295
 
 410  if (miter .eq. 0 .or. jcur .eq. 1) go to 430
296
 
      icf = 1
297
 
      ipup = miter
298
 
      go to 220
299
 
 430  icf = 2
300
 
      ncf = ncf + 1
301
 
      rmax = 2.0d0
302
 
      tn = told
303
 
      i1 = nqnyh + 1
304
 
      do 445 jb = 1,nq
305
 
        i1 = i1 - nyh
306
 
cdir$ ivdep
307
 
        do 440 i = i1,nqnyh
308
 
 440      yh1(i) = yh1(i) - yh1(i+nyh)
309
 
 445    continue
310
 
      if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680
311
 
      if (dabs(h) .le. hmin*1.00001d0) go to 670
312
 
      if (ncf .eq. mxncf) go to 670
313
 
      rh = 0.25d0
314
 
      ipup = miter
315
 
      iredo = 1
316
 
      go to 170
317
 
c-----------------------------------------------------------------------
318
 
c the corrector has converged.  jcur is set to 0
319
 
c to signal that the jacobian involved may need updating later.
320
 
c the local error test is made and control passes to statement 500
321
 
c if it fails.
322
 
c-----------------------------------------------------------------------
323
 
 450  jcur = 0
324
 
      if (m .eq. 0) dsm = del/tesco(2,nq)
325
 
      if (m .gt. 0) dsm = vnorm (n, acor, ewt)/tesco(2,nq)
326
 
      if (dsm .gt. 1.0d0) go to 500
327
 
c-----------------------------------------------------------------------
328
 
c after a successful step, update the yh array.
329
 
c consider changing h if ialth = 1.  otherwise decrease ialth by 1.
330
 
c if ialth is then 1 and nq .lt. maxord, then acor is saved for
331
 
c use in a possible order increase on the next step.
332
 
c if a change in h is considered, an increase or decrease in order
333
 
c by one is considered also.  a change in h is made only if it is by a
334
 
c factor of at least 1.1.  if not, ialth is set to 3 to prevent
335
 
c testing for that many steps.
336
 
c-----------------------------------------------------------------------
337
 
      kflag = 0
338
 
      iredo = 0
339
 
      nst = nst + 1
340
 
      hu = h
341
 
      nqu = nq
342
 
      do 470 j = 1,l
343
 
        do 470 i = 1,n
344
 
 470      yh(i,j) = yh(i,j) + el(j)*acor(i)
345
 
      ialth = ialth - 1
346
 
      if (ialth .eq. 0) go to 520
347
 
      if (ialth .gt. 1) go to 700
348
 
      if (l .eq. lmax) go to 700
349
 
      do 490 i = 1,n
350
 
 490    yh(i,lmax) = acor(i)
351
 
      go to 700
352
 
c-----------------------------------------------------------------------
353
 
c the error test failed.  kflag keeps track of multiple failures.
354
 
c restore tn and the yh array to their previous values, and prepare
355
 
c to try the step again.  compute the optimum step size for this or
356
 
c one lower order.  after 2 or more failures, h is forced to decrease
357
 
c by a factor of 0.2 or less.
358
 
c-----------------------------------------------------------------------
359
 
 500  kflag = kflag - 1
360
 
      tn = told
361
 
      i1 = nqnyh + 1
362
 
      do 515 jb = 1,nq
363
 
        i1 = i1 - nyh
364
 
cdir$ ivdep
365
 
        do 510 i = i1,nqnyh
366
 
 510      yh1(i) = yh1(i) - yh1(i+nyh)
367
 
 515    continue
368
 
      rmax = 2.0d0
369
 
      if (dabs(h) .le. hmin*1.00001d0) go to 660
370
 
      if (kflag .le. -3) go to 640
371
 
      iredo = 2
372
 
      rhup = 0.0d0
373
 
      go to 540
374
 
c-----------------------------------------------------------------------
375
 
c regardless of the success or failure of the step, factors
376
 
c rhdn, rhsm, and rhup are computed, by which h could be multiplied
377
 
c at order nq - 1, order nq, or order nq + 1, respectively.
378
 
c in the case of failure, rhup = 0.0 to avoid an order increase.
379
 
c the largest of these is determined and the new order chosen
380
 
c accordingly.  if the order is to be increased, we compute one
381
 
c additional scaled derivative.
382
 
c-----------------------------------------------------------------------
383
 
 520  rhup = 0.0d0
384
 
      if (l .eq. lmax) go to 540
385
 
      do 530 i = 1,n
386
 
 530    savf(i) = acor(i) - yh(i,lmax)
387
 
      dup = vnorm (n, savf, ewt)/tesco(3,nq)
388
 
      exup = 1.0d0/dfloat(l+1)
389
 
      rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
390
 
 540  exsm = 1.0d0/dfloat(l)
391
 
      rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
392
 
      rhdn = 0.0d0
393
 
      if (nq .eq. 1) go to 560
394
 
      ddn = vnorm (n, yh(1,l), ewt)/tesco(1,nq)
395
 
      exdn = 1.0d0/dfloat(nq)
396
 
      rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
397
 
 560  if (rhsm .ge. rhup) go to 570
398
 
      if (rhup .gt. rhdn) go to 590
399
 
      go to 580
400
 
 570  if (rhsm .lt. rhdn) go to 580
401
 
      newq = nq
402
 
      rh = rhsm
403
 
      go to 620
404
 
 580  newq = nq - 1
405
 
      rh = rhdn
406
 
      if (kflag .lt. 0 .and. rh .gt. 1.0d0) rh = 1.0d0
407
 
      go to 620
408
 
 590  newq = l
409
 
      rh = rhup
410
 
      if (rh .lt. 1.1d0) go to 610
411
 
      r = el(l)/dfloat(l)
412
 
      do 600 i = 1,n
413
 
 600    yh(i,newq+1) = acor(i)*r
414
 
      go to 630
415
 
 610  ialth = 3
416
 
      go to 700
417
 
 620  if ((kflag .eq. 0) .and. (rh .lt. 1.1d0)) go to 610
418
 
      if (kflag .le. -2) rh = dmin1(rh,0.2d0)
419
 
c-----------------------------------------------------------------------
420
 
c if there is a change of order, reset nq, l, and the coefficients.
421
 
c in any case h is reset according to rh and the yh array is rescaled.
422
 
c then exit from 690 if the step was ok, or redo the step otherwise.
423
 
c-----------------------------------------------------------------------
424
 
      if (newq .eq. nq) go to 170
425
 
 630  nq = newq
426
 
      l = nq + 1
427
 
      iret = 2
428
 
      go to 150
429
 
c-----------------------------------------------------------------------
430
 
c control reaches this section if 3 or more failures have occured.
431
 
c if 10 failures have occurred, exit with kflag = -1.
432
 
c it is assumed that the derivatives that have accumulated in the
433
 
c yh array have errors of the wrong order.  hence the first
434
 
c derivative is recomputed, and the order is set to 1.  then
435
 
c h is reduced by a factor of 10, and the step is retried,
436
 
c until it succeeds or h reaches hmin.
437
 
c-----------------------------------------------------------------------
438
 
 640  if (kflag .eq. -10) go to 660
439
 
      rh = 0.1d0
440
 
      rh = dmax1(hmin/dabs(h),rh)
441
 
      h = h*rh
442
 
      do 645 i = 1,n
443
 
 645    y(i) = yh(i,1)
444
 
      call f (neq, tn, y, savf)
445
 
      nfe = nfe + 1
446
 
      do 650 i = 1,n
447
 
 650    yh(i,2) = h*savf(i)
448
 
      ipup = miter
449
 
      ialth = 5
450
 
      if (nq .eq. 1) go to 200
451
 
      nq = 1
452
 
      l = 2
453
 
      iret = 3
454
 
      go to 150
455
 
c-----------------------------------------------------------------------
456
 
c all returns are made through this section.  h is saved in hold
457
 
c to allow the caller to change h on the next step.
458
 
c-----------------------------------------------------------------------
459
 
 660  kflag = -1
460
 
      go to 720
461
 
 670  kflag = -2
462
 
      go to 720
463
 
 680  kflag = -3
464
 
      go to 720
465
 
 690  rmax = 10.0d0
466
 
 700  r = 1.0d0/tesco(2,nqu)
467
 
      do 710 i = 1,n
468
 
 710    acor(i) = acor(i)*r
469
 
 720  hold = h
470
 
      jstart = 1
471
 
      return
472
 
c----------------------- end of subroutine stode -----------------------
473
 
      end