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,
77 IMPLICIT DOUBLE PRECISION(a-h,o-z)
79 dimension phi(neq,*),savr(*),delta(*),e(*)
81 dimension psi(*),alpha(*),beta(*),gamma(*),sigma(*)
83 EXTERNAL res, jac, psol, nls
86 parameter(lnst=11, letf=14, lcfn=15)
101 IF(jstart .NE. 0)
GO TO 120
128 IF(h.NE.hold.OR.k .NE. kold) ns = 0
131 IF(kp1 .LT. ns)
GO TO 230
141 beta(i)=beta(i-1)*psi(i-1)/temp2
144 sigma(i)=(i-1)*sigma(i-1)*alpha(i)
145 gamma(i)=gamma(i-1)+alpha(i-1)/h
155 alphas = alphas - 1.0d0/i
156 alpha0 = alpha0 - alpha(i)
166 ck =
abs(alpha(kp1) + alphas - alpha0)
167 ck =
max(ck,alpha(kp1))
171 IF(kp1 .LT. nsp1)
GO TO 280
174 260 phi(i,j)=beta(j)*phi(i,j)
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)
202 IF(iernls .NE. 0)
GO TO 600
218 enorm = ddwnrm(neq,e,vt,rpar,ipar)
219 erk = sigma(k+1)*enorm
223 IF(k .EQ. 1)
GO TO 430
225 405 delta(i) = phi(i,kp1) + e(i)
226 erkm1=sigma(k)*ddwnrm(neq,delta,vt,rpar,ipar)
228 IF(k .GT. 2)
GO TO 410
229 IF(terkm1 .LE. 0.5*terk)
GO TO 420
233 415 delta(i) = phi(i,k) + delta(i)
234 erkm2=sigma(k-1)*ddwnrm(neq,delta,vt,rpar,ipar)
236 IF(
max(terkm1,terkm2).GT.terk)
GO TO 430
250 IF(err .GT. 1.0d0)
GO TO 600
264 iwm(lnst)=iwm(lnst)+1
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
282 510 delta(i)=e(i)-phi(i,kp2)
283 erkp1 = (1.0d0/(k+2))*ddwnrm(neq,delta,vt,rpar,ipar)
286 IF(terkp1.GE.0.5d0*terk)
GO TO 550
288 520
IF(terkm1.LE.
min(terk,terkp1))
GO TO 540
289 IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))
GO TO 550
317 r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
318 IF(r .LT. 2.0d0)
GO TO 555
321 555
IF(r .GT. 1.0d0)
GO TO 560
322 r =
max(0.5d0,
min(0.9d0,r))
330 IF(kold.EQ.iwm(lmxord))
GO TO 585
335 590 phi(i,kp1)=phi(i,kp1)+e(i)
339 595 phi(i,j)=phi(i,j)+phi(i,j+1)
360 IF(kp1.LT.nsp1)
GO TO 630
364 610 phi(i,j)=temp1*phi(i,j)
368 640 psi(i-1)=psi(i)-h
374 IF(iernls .EQ. 0)
GO TO 660
375 iwm(lcfn)=iwm(lcfn)+1
383 IF (iernls .LT. 0)
GO TO 675
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
398 iwm(letf)=iwm(letf)+1
399 IF (nef .GT. 1)
GO TO 665
407 r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
408 r =
max(0.25d0,
min(0.9d0,r))
410 IF (
abs(h) .GE. hmin)
GO TO 690
418 665
IF (nef .GT. 2)
GO TO 670
422 IF (
abs(h) .GE. hmin)
GO TO 690
432 IF (
abs(h) .GE. hmin)
GO TO 690
448 CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
450 IF (idid .GE. 0) idid = -7
457 690
IF (kold .EQ. 0)
THEN 460 695 phi(i,2) = r*phi(i,2)
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
charNDArray min(char d, const charNDArray &m)