GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ddstp.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE ddstp(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT,
6  * jstart,idid,rpar,ipar,phi,savr,delta,e,wm,iwm,
7  * alpha,beta,gamma,psi,sigma,cj,cjold,hold,s,hmin,uround,
8  * epli,sqrtn,rsqrtn,epcon,iphase,jcalc,jflg,k,kold,ns,nonneg,
9  * ntype,nls)
10 C
11 C***BEGIN PROLOGUE DDSTP
12 C***REFER TO DDASPK
13 C***DATE WRITTEN 890101 (YYMMDD)
14 C***REVISION DATE 900926 (YYMMDD)
15 C***REVISION DATE 940909 (YYMMDD) (Reset PSI(1), PHI(*,2) at 690)
16 C
17 C
18 C-----------------------------------------------------------------------
19 C***DESCRIPTION
20 C
21 C DDSTP solves a system of differential/algebraic equations of
22 C the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H).
23 C
24 C The methods used are modified divided difference, fixed leading
25 C coefficient forms of backward differentiation formulas.
26 C The code adjusts the stepsize and order to control the local error
27 C per step.
28 C
29 C
30 C The parameters represent
31 C X -- Independent variable.
32 C Y -- Solution vector at X.
33 C YPRIME -- Derivative of solution vector
34 C after successful step.
35 C NEQ -- Number of equations to be integrated.
36 C RES -- External user-supplied subroutine
37 C to evaluate the residual. See RES description
38 C in DDASPK prologue.
39 C JAC -- External user-supplied routine to update
40 C Jacobian or preconditioner information in the
41 C nonlinear solver. See JAC description in DDASPK
42 C prologue.
43 C PSOL -- External user-supplied routine to solve
44 C a linear system using preconditioning.
45 C (This is optional). See PSOL in DDASPK prologue.
46 C H -- Appropriate step size for next step.
47 C Normally determined by the code.
48 C WT -- Vector of weights for error criterion used in Newton test.
49 C VT -- Masked vector of weights used in error test.
50 C JSTART -- Integer variable set 0 for
51 C first step, 1 otherwise.
52 C IDID -- Completion code returned from the nonlinear solver.
53 C See IDID description in DDASPK prologue.
54 C RPAR,IPAR -- Real and integer parameter arrays that
55 C are used for communication between the
56 C calling program and external user routines.
57 C They are not altered by DNSK
58 C PHI -- Array of divided differences used by
59 C DDSTP. The length is NEQ*(K+1), where
60 C K is the maximum order.
61 C SAVR -- Work vector for DDSTP of length NEQ.
62 C DELTA,E -- Work vectors for DDSTP of length NEQ.
63 C WM,IWM -- Real and integer arrays storing
64 C information required by the linear solver.
65 C
66 C The other parameters are information
67 C which is needed internally by DDSTP to
68 C continue from step to step.
69 C
70 C-----------------------------------------------------------------------
71 C***ROUTINES CALLED
72 C NLS, DDWNRM, DDATRP
73 C
74 C***END PROLOGUE DDSTP
75 C
76 C
77  IMPLICIT DOUBLE PRECISION(a-h,o-z)
78  dimension y(*),yprime(*),wt(*),vt(*)
79  dimension phi(neq,*),savr(*),delta(*),e(*)
80  dimension wm(*),iwm(*)
81  dimension psi(*),alpha(*),beta(*),gamma(*),sigma(*)
82  dimension rpar(*),ipar(*)
83  EXTERNAL res, jac, psol, nls
84 C
85  parameter(lmxord=3)
86  parameter(lnst=11, letf=14, lcfn=15)
87 C
88 C
89 C-----------------------------------------------------------------------
90 C BLOCK 1.
91 C Initialize. On the first call, set
92 C the order to 1 and initialize
93 C other variables.
94 C-----------------------------------------------------------------------
95 C
96 C Initializations for all calls
97 C
98  xold=x
99  ncf=0
100  nef=0
101  IF(jstart .NE. 0) go to 120
102 C
103 C If this is the first step, perform
104 C other initializations
105 C
106  k=1
107  kold=0
108  hold=0.0d0
109  psi(1)=h
110  cj = 1.d0/h
111  iphase = 0
112  ns=0
113 120 CONTINUE
114 C
115 C
116 C
117 C
118 C
119 C-----------------------------------------------------------------------
120 C BLOCK 2
121 C Compute coefficients of formulas for
122 C this step.
123 C-----------------------------------------------------------------------
124 200 CONTINUE
125  kp1=k+1
126  kp2=k+2
127  km1=k-1
128  IF(h.NE.hold.OR.k .NE. kold) ns = 0
129  ns=min0(ns+1,kold+2)
130  nsp1=ns+1
131  IF(kp1 .LT. ns)go to 230
132 C
133  beta(1)=1.0d0
134  alpha(1)=1.0d0
135  temp1=h
136  gamma(1)=0.0d0
137  sigma(1)=1.0d0
138  DO 210 i=2,kp1
139  temp2=psi(i-1)
140  psi(i-1)=temp1
141  beta(i)=beta(i-1)*psi(i-1)/temp2
142  temp1=temp2+h
143  alpha(i)=h/temp1
144  sigma(i)=(i-1)*sigma(i-1)*alpha(i)
145  gamma(i)=gamma(i-1)+alpha(i-1)/h
146 210 CONTINUE
147  psi(kp1)=temp1
148 230 CONTINUE
149 C
150 C Compute ALPHAS, ALPHA0
151 C
152  alphas = 0.0d0
153  alpha0 = 0.0d0
154  DO 240 i = 1,k
155  alphas = alphas - 1.0d0/i
156  alpha0 = alpha0 - alpha(i)
157 240 CONTINUE
158 C
159 C Compute leading coefficient CJ
160 C
161  cjlast = cj
162  cj = -alphas/h
163 C
164 C Compute variable stepsize error coefficient CK
165 C
166  ck = abs(alpha(kp1) + alphas - alpha0)
167  ck = max(ck,alpha(kp1))
168 C
169 C Change PHI to PHI STAR
170 C
171  IF(kp1 .LT. nsp1) go to 280
172  DO 270 j=nsp1,kp1
173  DO 260 i=1,neq
174 260 phi(i,j)=beta(j)*phi(i,j)
175 270 CONTINUE
176 280 CONTINUE
177 C
178 C Update time
179 C
180  x=x+h
181 C
182 C Initialize IDID to 1
183 C
184  idid = 1
185 C
186 C
187 C
188 C
189 C
190 C-----------------------------------------------------------------------
191 C BLOCK 3
192 C Call the nonlinear system solver to obtain the solution and
193 C derivative.
194 C-----------------------------------------------------------------------
195 C
196  CALL nls(x,y,yprime,neq,
197  * res,jac,psol,h,wt,jstart,idid,rpar,ipar,phi,gamma,
198  * savr,delta,e,wm,iwm,cj,cjold,cjlast,s,
199  * uround,epli,sqrtn,rsqrtn,epcon,jcalc,jflg,kp1,
200  * nonneg,ntype,iernls)
201 C
202  IF(iernls .NE. 0)go to 600
203 C
204 C
205 C
206 C
207 C
208 C-----------------------------------------------------------------------
209 C BLOCK 4
210 C Estimate the errors at orders K,K-1,K-2
211 C as if constant stepsize was used. Estimate
212 C the local error at order K and test
213 C whether the current step is successful.
214 C-----------------------------------------------------------------------
215 C
216 C Estimate errors at orders K,K-1,K-2
217 C
218  enorm = ddwnrm(neq,e,vt,rpar,ipar)
219  erk = sigma(k+1)*enorm
220  terk = (k+1)*erk
221  est = erk
222  knew=k
223  IF(k .EQ. 1)go to 430
224  DO 405 i = 1,neq
225 405 delta(i) = phi(i,kp1) + e(i)
226  erkm1=sigma(k)*ddwnrm(neq,delta,vt,rpar,ipar)
227  terkm1 = k*erkm1
228  IF(k .GT. 2)go to 410
229  IF(terkm1 .LE. 0.5*terk)go to 420
230  go to 430
231 410 CONTINUE
232  DO 415 i = 1,neq
233 415 delta(i) = phi(i,k) + delta(i)
234  erkm2=sigma(k-1)*ddwnrm(neq,delta,vt,rpar,ipar)
235  terkm2 = (k-1)*erkm2
236  IF(max(terkm1,terkm2).GT.terk)go to 430
237 C
238 C Lower the order
239 C
240 420 CONTINUE
241  knew=k-1
242  est = erkm1
243 C
244 C
245 C Calculate the local error for the current step
246 C to see if the step was successful
247 C
248 430 CONTINUE
249  err = ck * enorm
250  IF(err .GT. 1.0d0)go to 600
251 C
252 C
253 C
254 C
255 C
256 C-----------------------------------------------------------------------
257 C BLOCK 5
258 C The step is successful. Determine
259 C the best order and stepsize for
260 C the next step. Update the differences
261 C for the next step.
262 C-----------------------------------------------------------------------
263  idid=1
264  iwm(lnst)=iwm(lnst)+1
265  kdiff=k-kold
266  kold=k
267  hold=h
268 C
269 C
270 C Estimate the error at order K+1 unless
271 C already decided to lower order, or
272 C already using maximum order, or
273 C stepsize not constant, or
274 C order raised in previous step
275 C
276  IF(knew.EQ.km1.OR.k.EQ.iwm(lmxord))iphase=1
277  IF(iphase .EQ. 0)go to 545
278  IF(knew.EQ.km1)go to 540
279  IF(k.EQ.iwm(lmxord)) go to 550
280  IF(kp1.GE.ns.OR.kdiff.EQ.1)go to 550
281  DO 510 i=1,neq
282 510 delta(i)=e(i)-phi(i,kp2)
283  erkp1 = (1.0d0/(k+2))*ddwnrm(neq,delta,vt,rpar,ipar)
284  terkp1 = (k+2)*erkp1
285  IF(k.GT.1)go to 520
286  IF(terkp1.GE.0.5d0*terk)go to 550
287  go to 530
288 520 IF(terkm1.LE.min(terk,terkp1))go to 540
289  IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))go to 550
290 C
291 C Raise order
292 C
293 530 k=kp1
294  est = erkp1
295  go to 550
296 C
297 C Lower order
298 C
299 540 k=km1
300  est = erkm1
301  go to 550
302 C
303 C If IPHASE = 0, increase order by one and multiply stepsize by
304 C factor two
305 C
306 545 k = kp1
307  hnew = h*2.0d0
308  h = hnew
309  go to 575
310 C
311 C
312 C Determine the appropriate stepsize for
313 C the next step.
314 C
315 550 hnew=h
316  temp2=k+1
317  r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
318  IF(r .LT. 2.0d0) go to 555
319  hnew = 2.0d0*h
320  go to 560
321 555 IF(r .GT. 1.0d0) go to 560
322  r = max(0.5d0,min(0.9d0,r))
323  hnew = h*r
324 560 h=hnew
325 C
326 C
327 C Update differences for next step
328 C
329 575 CONTINUE
330  IF(kold.EQ.iwm(lmxord))go to 585
331  DO 580 i=1,neq
332 580 phi(i,kp2)=e(i)
333 585 CONTINUE
334  DO 590 i=1,neq
335 590 phi(i,kp1)=phi(i,kp1)+e(i)
336  DO 595 j1=2,kp1
337  j=kp1-j1+1
338  DO 595 i=1,neq
339 595 phi(i,j)=phi(i,j)+phi(i,j+1)
340  jstart = 1
341  RETURN
342 C
343 C
344 C
345 C
346 C
347 C-----------------------------------------------------------------------
348 C BLOCK 6
349 C The step is unsuccessful. Restore X,PSI,PHI
350 C Determine appropriate stepsize for
351 C continuing the integration, or exit with
352 C an error flag if there have been many
353 C failures.
354 C-----------------------------------------------------------------------
355 600 iphase = 1
356 C
357 C Restore X,PHI,PSI
358 C
359  x=xold
360  IF(kp1.LT.nsp1)go to 630
361  DO 620 j=nsp1,kp1
362  temp1=1.0d0/beta(j)
363  DO 610 i=1,neq
364 610 phi(i,j)=temp1*phi(i,j)
365 620 CONTINUE
366 630 CONTINUE
367  DO 640 i=2,kp1
368 640 psi(i-1)=psi(i)-h
369 C
370 C
371 C Test whether failure is due to nonlinear solver
372 C or error test
373 C
374  IF(iernls .EQ. 0)go to 660
375  iwm(lcfn)=iwm(lcfn)+1
376 C
377 C
378 C The nonlinear solver failed to converge.
379 C Determine the cause of the failure and take appropriate action.
380 C If IERNLS .LT. 0, then return. Otherwise, reduce the stepsize
381 C and try again, unless too many failures have occurred.
382 C
383  IF (iernls .LT. 0) go to 675
384  ncf = ncf + 1
385  r = 0.25d0
386  h = h*r
387  IF (ncf .LT. 10 .AND. abs(h) .GE. hmin) go to 690
388  IF (idid .EQ. 1) idid = -7
389  IF (nef .GE. 3) idid = -9
390  go to 675
391 C
392 C
393 C The nonlinear solver converged, and the cause
394 C of the failure was the error estimate
395 C exceeding the tolerance.
396 C
397 660 nef=nef+1
398  iwm(letf)=iwm(letf)+1
399  IF (nef .GT. 1) go to 665
400 C
401 C On first error test failure, keep current order or lower
402 C order by one. Compute new stepsize based on differences
403 C of the solution.
404 C
405  k = knew
406  temp2 = k + 1
407  r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
408  r = max(0.25d0,min(0.9d0,r))
409  h = h*r
410  IF (abs(h) .GE. hmin) go to 690
411  idid = -6
412  go to 675
413 C
414 C On second error test failure, use the current order or
415 C decrease order by one. Reduce the stepsize by a factor of
416 C one quarter.
417 C
418 665 IF (nef .GT. 2) go to 670
419  k = knew
420  r = 0.25d0
421  h = r*h
422  IF (abs(h) .GE. hmin) go to 690
423  idid = -6
424  go to 675
425 C
426 C On third and subsequent error test failures, set the order to
427 C one, and reduce the stepsize by a factor of one quarter.
428 C
429 670 k = 1
430  r = 0.25d0
431  h = r*h
432  IF (abs(h) .GE. hmin) go to 690
433  idid = -6
434  go to 675
435 C
436 C
437 C
438 C
439 C For all crashes, restore Y to its last value,
440 C interpolate to find YPRIME at last X, and return.
441 C
442 C Before returning, verify that the user has not set
443 C IDID to a nonnegative value. If the user has set IDID
444 C to a nonnegative value, then reset IDID to be -7, indicating
445 C a failure in the nonlinear system solver.
446 C
447 675 CONTINUE
448  CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
449  jstart = 1
450  IF (idid .GE. 0) idid = -7
451  RETURN
452 C
453 C
454 C Go back and try this step again.
455 C If this is the first step, reset PSI(1) and rescale PHI(*,2).
456 C
457 690 IF (kold .EQ. 0) THEN
458  psi(1) = h
459  DO 695 i = 1,neq
460 695 phi(i,2) = r*phi(i,2)
461  ENDIF
462  go to 200
463 C
464 C------END OF SUBROUTINE DDSTP------------------------------------------
465  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2