| [neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,df,djac,rpar,ipar]=dstod(neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,df,djac,rpar,ipar); |
function [neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,df,djac,rpar,ipar]=dstod(neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,df,djac,rpar,ipar);
%***BEGIN PROLOGUE DSTOD
%***SUBSIDIARY
%***PURPOSE Subsidiary to DDEBDF
%***LIBRARY SLATEC
%***TYPE doubleprecision (STOD-S, DSTOD-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% DSTOD integrates a system of first order odes over one step in the
% integrator package DDEBDF.
% ----------------------------------------------------------------------
% DSTOD performs one step of the integration of an initial value
% problem for a system of ordinary differential equations.
% Note.. DSTOD is independent of the value of the iteration method
% indicator MITER, when this is .NE. 0, and hence is independent
% of the type of chord method used, or the Jacobian structure.
% Communication with DSTOD is done with the following variables..
%
% Y = An array of length .GE. N used as the Y argument in
% all calls to DF and DJAC.
% NEQ = Integer array containing problem size in NEQ(1), and
% passed as the NEQ argument in all calls to DF and DJAC.
% YH = An NYH by LMAX array containing the dependent variables
% and their approximate scaled derivatives, where
% LMAX = MAXORD + 1. YH(I,J+1) contains the approximate
% J-th derivative of Y(I), scaled by H**J/FACTORIAL(J)
% (J = 0,1,...,NQ). On entry for the first step, the first
% two columns of YH must be set from the initial values.
% NYH = A constant integer .GE. N, the first dimension of YH.
% YH1 = A one-dimensional array occupying the same space as YH.
% EWT = An array of N elements with which the estimated local
% errors in YH are compared.
% SAVF = An array of working storage, of length N.
% ACOR = A work array of length N, used for the accumulated
% corrections. On a successful return, ACOR(I) contains
% the estimated one-step local error in Y(I).
% WM,IWM = doubleprecision and INTEGER work arrays associated with
% matrix operations in chord iteration (MITER .NE. 0).
% DPJAC = Name of routine to evaluate and preprocess Jacobian matrix
% if a chord method is being used.
% DSLVS = Name of routine to solve linear system in chord iteration.
% H = The step size to be attempted on the next step.
% H is altered by the error control algorithm during the
% problem. H can be either positive or negative, but its
% sign must remain constant throughout the problem.
% HMIN = The minimum absolute value of the step size H to be used.
% HMXI = Inverse of the maximum absolute value of H to be used.
% HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
% HMIN and HMXI may be changed at any time, but will not
% take effect until the next change of H is considered.
% TN = The independent variable. TN is updated on each step taken.
% JSTART = An integer used for input only, with the following
% values and meanings..
% 0 Perform the first step.
% .GT.0 Take a new step continuing from the last.
% -1 Take the next step with a new value of H, MAXORD,
% N, METH, MITER, and/or matrix parameters.
% -2 Take the next step with a new value of H,
% but with other inputs unchanged.
% On return, JSTART is set to 1 to facilitate continuation.
% KFLAG = a completion code with the following meanings..
% 0 The step was successful.
% -1 The requested error could not be achieved.
% -2 Corrector convergence could not be achieved.
% A return with KFLAG = -1 or -2 means either
% ABS(H) = HMIN or 10 consecutive failures occurred.
% On a return with KFLAG negative, the values of TN and
% the YH array are as of the beginning of the last
% step, and H is the last step size attempted.
% MAXORD = The maximum order of integration method to be allowed.
% METH/MITER = The method flags. See description in driver.
% N = The number of first-order differential equations.
% ----------------------------------------------------------------------
%
%***SEE ALSO DDEBDF
%***ROUTINES CALLED DCFOD, DPJAC, DSLVS, DVNRMS
%***COMMON BLOCKS DDEBD1
%***REVISION HISTORY (YYMMDD)
% 820301 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 910722 Updated AUTHOR section. (ALS)
% 920422 Changed DIMENSION statement. (WRB)
%***end PROLOGUE DSTOD
%
persistent dcon ddn del delp dsm dup exdn exsm exup gt i i1 iredo iret j jb m ncf newq r rh rhdn rhsm rhup told ;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
global ddebd1_20; if isempty(ddebd1_20), ddebd1_20=0; end;
global ddebd1_26; if isempty(ddebd1_26), ddebd1_26=0; end;
global ddebd1_19; if isempty(ddebd1_19), ddebd1_19=zeros(1,6); end;
global ddebd1_17; if isempty(ddebd1_17), ddebd1_17=zeros(1,7); end;
global ddebd1_21; if isempty(ddebd1_21), ddebd1_21=0; end;
if isempty(iredo), iredo=0; end;
if isempty(iret), iret=0; end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
global ddebd1_27; if isempty(ddebd1_27), ddebd1_27=0; end;
global ddebd1_28; if isempty(ddebd1_28), ddebd1_28=0; end;
global ddebd1_18; if isempty(ddebd1_18), ddebd1_18=0; end;
global ddebd1_29; if isempty(ddebd1_29), ddebd1_29=0; end;
global ddebd1_22; if isempty(ddebd1_22), ddebd1_22=0; end;
if isempty(m), m=0; end;
global ddebd1_32; if isempty(ddebd1_32), ddebd1_32=0; end;
global ddebd1_23; if isempty(ddebd1_23), ddebd1_23=0; end;
global ddebd1_30; if isempty(ddebd1_30), ddebd1_30=0; end;
global ddebd1_31; if isempty(ddebd1_31), ddebd1_31=0; end;
global ddebd1_33; if isempty(ddebd1_33), ddebd1_33=0; end;
if isempty(ncf), ncf=0; end;
if isempty(newq), newq=0; end;
global ddebd1_36; if isempty(ddebd1_36), ddebd1_36=0; end;
global ddebd1_37; if isempty(ddebd1_37), ddebd1_37=0; end;
global ddebd1_34; if isempty(ddebd1_34), ddebd1_34=0; end;
global ddebd1_24; if isempty(ddebd1_24), ddebd1_24=0; end;
global ddebd1_38; if isempty(ddebd1_38), ddebd1_38=0; end;
global ddebd1_35; if isempty(ddebd1_35), ddebd1_35=0; end;
global ddebd1_25; if isempty(ddebd1_25), ddebd1_25=0; end;
if isempty(gt), gt=zeros(1,9); end;
global ddebd1_2; if isempty(ddebd1_2), ddebd1_2=0; end;
global ddebd1_3; if isempty(ddebd1_3), ddebd1_3=0; end;
if isempty(dcon), dcon=0; end;
if isempty(ddn), ddn=0; end;
if isempty(del), del=0; end;
if isempty(delp), delp=0; end;
if isempty(dsm), dsm=0; end;
if isempty(dup), dup=0; end;
global ddebd1_4; if isempty(ddebd1_4), ddebd1_4=zeros(1,13); end;
global ddebd1_10; if isempty(ddebd1_10), ddebd1_10=0; end;
global ddebd1_5; if isempty(ddebd1_5), ddebd1_5=zeros(13,12); end;
if isempty(exdn), exdn=0; end;
if isempty(exsm), exsm=0; end;
if isempty(exup), exup=0; end;
global ddebd1_11; if isempty(ddebd1_11), ddebd1_11=0; end;
global ddebd1_12; if isempty(ddebd1_12), ddebd1_12=0; end;
global ddebd1_13; if isempty(ddebd1_13), ddebd1_13=0; end;
global ddebd1_6; if isempty(ddebd1_6), ddebd1_6=0; end;
global ddebd1_14; if isempty(ddebd1_14), ddebd1_14=0; end;
if isempty(r), r=0; end;
global ddebd1_7; if isempty(ddebd1_7), ddebd1_7=0; end;
if isempty(rh), rh=0; end;
if isempty(rhdn), rhdn=0; end;
if isempty(rhsm), rhsm=0; end;
if isempty(rhup), rhup=0; end;
global ddebd1_8; if isempty(ddebd1_8), ddebd1_8=0; end;
global ddebd1_1; if isempty(ddebd1_1), ddebd1_1=0; end;
global ddebd1_9; if isempty(ddebd1_9), ddebd1_9=zeros(3,12); end;
global ddebd1_15; if isempty(ddebd1_15), ddebd1_15=0; end;
if isempty(told), told=0; end;
global ddebd1_16; if isempty(ddebd1_16), ddebd1_16=0; end;
%
y_shape=size(y);y=reshape(y,1,[]);
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([nyh])).*prod([nyh])-numel(yh))],nyh,[]);
yh1_shape=size(yh1);yh1=reshape(yh1,1,[]);
ewt_shape=size(ewt);ewt=reshape(ewt,1,[]);
savf_shape=size(savf);savf=reshape(savf,1,[]);
acor_shape=size(acor);acor=reshape(acor,1,[]);
wm_shape=size(wm);wm=reshape(wm,1,[]);
iwm_shape=size(iwm);iwm=reshape(iwm,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
% common :: ;
%% common /ddebd1/ rownd , conit , crate , el(13) , elco(13,12) ,hold , rc , rmax , tesco(3,12) , el0 , h , hmin ,hmxi , hu , tn , uround , iownd(7) , ksteps ,iod(6) , ialth , ipup , lmax , meo , nqnyh ,nstepj , ier , jstart , kflag , l , meth , miter ,maxord , n , nq , nst , nfe , nje , nqu;
%% common /ddebd1/ ddebd1_1 , ddebd1_2 , ddebd1_3 , ddebd1_4(13) , ddebd1_5(13,12) ,ddebd1_6 , ddebd1_7 , ddebd1_8 , ddebd1_9(3,12) , ddebd1_10 , ddebd1_11 , ddebd1_12 ,ddebd1_13 , ddebd1_14 , ddebd1_15 , ddebd1_16 , ddebd1_17(7) , ddebd1_18 ,ddebd1_19(6) , ddebd1_20 , ddebd1_21 , ddebd1_22 , ddebd1_23 , ddebd1_24 ,ddebd1_25 , ddebd1_26 , ddebd1_27 , ddebd1_28 , ddebd1_29 , ddebd1_30 , ddebd1_31 ,ddebd1_32 , ddebd1_33 , ddebd1_34 , ddebd1_35 , ddebd1_36 , ddebd1_37 , ddebd1_38;
%
%
% BEGIN BLOCK PERMITTING ...EXITS TO 690
% BEGIN BLOCK PERMITTING ...EXITS TO 60
%***FIRST EXECUTABLE STATEMENT DSTOD
gt(:)=0;
ddebd1_28 = 0;
told = ddebd1_15;
ncf = 0;
while (1);
if( ddebd1_27>0 )
gt(4)=1;
break;
end;
if( ddebd1_27==-1 )
% BEGIN BLOCK PERMITTING ...EXITS TO 30
% ------------------------------------------------------
% THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN
% JSTART = -1. IPUP IS SET TO MITER TO FORCE A MATRIX
% UPDATE. IF AN ORDER INCREASE IS ABOUT TO BE
% CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO
% POSTPONE CONSIDERATION ONE MORE STEP. IF THE CALLER
% HAS CHANGED METH, DCFOD IS CALLED TO RESET THE
% COEFFICIENTS OF THE METHOD. IF THE CALLER HAS
% CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
% ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN
% ACCORDINGLY. IF H IS TO BE CHANGED, YH MUST BE
% RESCALED. IF H OR METH IS BEING CHANGED, IALTH IS
% RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H
% FOR THAT MANY STEPS.
% ------------------------------------------------------
ddebd1_21 = fix(ddebd1_31);
ddebd1_22 = fix(ddebd1_32 + 1);
if( ddebd1_20==1 )
ddebd1_20 = 2;
end;
if( ddebd1_30~=ddebd1_23 )
[ddebd1_30,ddebd1_5,ddebd1_9]=dcfod(ddebd1_30,ddebd1_5,ddebd1_9);
ddebd1_23 = fix(ddebd1_30);
% ......EXIT
if( ddebd1_34<=ddebd1_32 )
ddebd1_20 = fix(ddebd1_29);
iret = 1;
% ............EXIT
break;
end;
elseif( ddebd1_34<=ddebd1_32 ) ;
gt(2)=1;
break;
end;
ddebd1_34 = fix(ddebd1_32);
ddebd1_29 = fix(ddebd1_22);
for i = 1 : ddebd1_29;
ddebd1_4(i) = ddebd1_5(i,ddebd1_34);
end; i = fix(ddebd1_29+1);
ddebd1_24 = fix(ddebd1_34.*nyh);
ddebd1_7 = ddebd1_7.*ddebd1_4(1)./ddebd1_10;
ddebd1_10 = ddebd1_4(1);
ddebd1_2 = 0.5d0./(ddebd1_34+2);
ddn = dvnrms(ddebd1_33,savf,ewt)./ddebd1_9(1,ddebd1_29);
exdn = 1.0d0./ddebd1_29;
rhdn = 1.0d0./(1.3d0.*ddn.^exdn+0.0000013d0);
rh = min(rhdn,1.0d0);
iredo = 3;
if( ddebd1_11==ddebd1_6 )
rh = max(rh,ddebd1_12./abs(ddebd1_11));
else;
rh = min(rh,abs(ddebd1_11./ddebd1_6));
ddebd1_11 = ddebd1_6;
end;
gt(3)=1;
break;
else;
if( ddebd1_27==-2 )
gt(2)=1;
break;
end;
% ---------------------------------------------------------
% ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER
% VARIABLES ARE INITIALIZED. RMAX IS THE MAXIMUM RATIO BY
% WHICH H CAN BE INCREASED IN A SINGLE STEP. IT IS
% INITIALLY 1.0E4 TO COMPENSATE FOR THE SMALL INITIAL H,
% BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE OCCURS
% (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT
% 2 FOR THE NEXT INCREASE.
% ---------------------------------------------------------
ddebd1_22 = fix(ddebd1_32 + 1);
ddebd1_34 = 1;
ddebd1_29 = 2;
ddebd1_20 = 2;
ddebd1_8 = 10000.0d0;
ddebd1_7 = 0.0d0;
ddebd1_10 = 1.0d0;
ddebd1_3 = 0.7d0;
delp = 0.0d0;
ddebd1_6 = ddebd1_11;
ddebd1_23 = fix(ddebd1_30);
ddebd1_25 = 0;
iret = 3;
% ------------------------------------------------------------
% DCFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS
% FOR THE CURRENT METH. THEN THE EL VECTOR AND RELATED
% CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT
% THE START OF THE PROBLEM.
% ------------------------------------------------------------
[ddebd1_30,ddebd1_5,ddebd1_9]=dcfod(ddebd1_30,ddebd1_5,ddebd1_9);
end;
break;
end;
% BEGIN BLOCK PERMITTING ...EXITS TO 680
while (1);
if(gt(9)==0)
if(gt(8)==0)
if(gt(7)==0)
if(gt(6)==0)
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
for i = 1 : ddebd1_29;
ddebd1_4(i) = ddebd1_5(i,ddebd1_34);
end; i = fix(ddebd1_29+1);
ddebd1_24 = fix(ddebd1_34.*nyh);
ddebd1_7 = ddebd1_7.*ddebd1_4(1)./ddebd1_10;
ddebd1_10 = ddebd1_4(1);
ddebd1_2 = 0.5d0./(ddebd1_34+2);
if( iret==2 )
rh = max(rh,ddebd1_12./abs(ddebd1_11));
gt(3)=1;
continue;
elseif( iret==3 ) ;
gt(4)=1;
continue;
end;
% ---------------------------------------------------------
% IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
% RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH
% IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT
% MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST
% FAILURE.
% ---------------------------------------------------------
end;
gt(2)=0;
if( ddebd1_11==ddebd1_6 )
gt(4)=1;
continue;
end;
rh = ddebd1_11./ddebd1_6;
ddebd1_11 = ddebd1_6;
iredo = 3;
end;
gt(3)=0;
rh = min(rh,ddebd1_8);
rh = rh./max(1.0d0,abs(ddebd1_11).*ddebd1_13.*rh);
r = 1.0d0;
for j = 2 : ddebd1_29;
r = r.*rh;
for i = 1 : ddebd1_33;
yh(i,j) = yh(i,j).*r;
end; i = fix(ddebd1_33+1);
end; j = fix(ddebd1_29+1);
ddebd1_11 = ddebd1_11.*rh;
ddebd1_7 = ddebd1_7.*rh;
ddebd1_20 = fix(ddebd1_29);
if( iredo==0 )
ddebd1_8 = 10.0d0;
r = 1.0d0./ddebd1_9(2,ddebd1_38);
for i = 1 : ddebd1_33;
acor(i) = acor(i).*r;
end; i = fix(ddebd1_33+1);
% ...............EXIT
break;
end;
% ------------------------------------------------------
% THIS SECTION COMPUTES THE PREDICTED VALUES BY
% EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL
% TRIANGLE MATRIX. RC IS THE RATIO OF NEW TO OLD
% VALUES OF THE COEFFICIENT H*EL(1). WHEN RC DIFFERS
% FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
% TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS
% INVOLVED. IN ANY CASE, DPJAC IS CALLED AT LEAST
% EVERY 20-TH STEP.
% ------------------------------------------------------
% BEGIN BLOCK PERMITTING ...EXITS TO 610
% BEGIN BLOCK PERMITTING ...EXITS TO 490
end;
gt(4)=0;
if( abs(ddebd1_7-1.0d0)>0.3d0 )
ddebd1_21 = fix(ddebd1_31);
end;
if( ddebd1_35>=ddebd1_25+20 )
ddebd1_21 = fix(ddebd1_31);
end;
ddebd1_15 = ddebd1_15 + ddebd1_11;
i1 = fix(ddebd1_24 + 1);
for jb = 1 : ddebd1_34;
i1 = fix(i1 - nyh);
for i = i1 : ddebd1_24;
yh1(i) = yh1(i) + yh1(i+nyh);
end; i = fix(ddebd1_24+1);
end; jb = fix(ddebd1_34+1);
ddebd1_18 = fix(ddebd1_18 + 1);
% ---------------------------------------------
% UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A
% CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
% OF EACH CORRECTION, WEIGHTED BY THE ERROR
% WEIGHT VECTOR EWT. THE SUM OF THE
% CORRECTIONS IS ACCUMULATED IN THE VECTOR
% ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE
% CORRECTOR LOOP.
% ---------------------------------------------
end;
gt(5)=0;
m = 0;
for i = 1 : ddebd1_33;
y(i) = yh(i,1);
end; i = fix(ddebd1_33+1);
[ddebd1_15,y,savf,rpar,ipar]=df(ddebd1_15,y,savf,rpar,ipar);
ddebd1_36 = fix(ddebd1_36 + 1);
if( ddebd1_21>0 )
% ---------------------------------------
% IF INDICATED, THE MATRIX P = I -
% H*EL(1)*J IS REEVALUATED AND
% PREPROCESSED BEFORE STARTING THE
% CORRECTOR ITERATION. IPUP IS SET TO 0
% AS AN INDICATOR THAT THIS HAS BEEN
% DONE.
% ---------------------------------------
ddebd1_21 = 0;
ddebd1_7 = 1.0d0;
ddebd1_25 = fix(ddebd1_35);
ddebd1_3 = 0.7d0;
[neq,y,yh,nyh,ewt,acor,savf,wm,iwm,df,djac,rpar,ipar]=dpjac(neq,y,yh,nyh,ewt,acor,savf,wm,iwm,df,djac,rpar,ipar);
% ......EXIT
if( ddebd1_26~=0 )
gt(8)=1;
continue;
end;
end;
for i = 1 : ddebd1_33;
acor(i) = 0.0d0;
end; i = fix(ddebd1_33+1);
end;
gt(6)=0;
if( ddebd1_31~=0 )
% ------------------------------------
% IN THE CASE OF THE CHORD METHOD,
% COMPUTE THE CORRECTOR ERROR, AND
% SOLVE THE LINEAR SYSTEM WITH THAT
% AS RIGHT-HAND SIDE AND P AS
% COEFFICIENT MATRIX.
% ------------------------------------
for i = 1 : ddebd1_33;
y(i) = ddebd1_11.*savf(i) -(yh(i,2)+acor(i));
end; i = fix(ddebd1_33+1);
[wm,iwm,y,savf]=dslvs(wm,iwm,y,savf);
% ......EXIT
if( ddebd1_26~=0 )
gt(7)=1;
continue;
end;
[del ,ddebd1_33,y,ewt]=dvnrms(ddebd1_33,y,ewt);
for i = 1 : ddebd1_33;
acor(i) = acor(i) + y(i);
y(i) = yh(i,1) + ddebd1_4(1).*acor(i);
end; i = fix(ddebd1_33+1);
else;
% ------------------------------------
% IN THE CASE OF FUNCTIONAL
% ITERATION, UPDATE Y DIRECTLY FROM
% THE RESULT OF THE LAST FUNCTION
% EVALUATION.
% ------------------------------------
for i = 1 : ddebd1_33;
savf(i) = ddebd1_11.*savf(i) - yh(i,2);
y(i) = savf(i) - acor(i);
end; i = fix(ddebd1_33+1);
[del ,ddebd1_33,y,ewt]=dvnrms(ddebd1_33,y,ewt);
for i = 1 : ddebd1_33;
y(i) = yh(i,1) + ddebd1_4(1).*savf(i);
acor(i) = savf(i);
end; i = fix(ddebd1_33+1);
end;
% ---------------------------------------
% TEST FOR CONVERGENCE. IF M.GT.0, AN
% ESTIMATE OF THE CONVERGENCE RATE
% CONSTANT IS STORED IN CRATE, AND THIS
% IS USED IN THE TEST.
% ---------------------------------------
if( m~=0 )
ddebd1_3 = max(0.2d0.*ddebd1_3,del./delp);
end;
dcon = del.*min(1.0d0,1.5d0.*ddebd1_3)./(ddebd1_9(2,ddebd1_34).*ddebd1_2);
if( dcon>1.0d0 )
m = fix(m + 1);
% ...EXIT
if( m~=3 )
% ...EXIT
if( m<2 || del<=2.0d0.*delp )
delp = del;
[ddebd1_15,y,savf,rpar,ipar]=df(ddebd1_15,y,savf,rpar,ipar);
ddebd1_36 = fix(ddebd1_36 + 1);
gt(6)=1;
continue;
end;
end;
else;
% ------------------------------------
% THE CORRECTOR HAS CONVERGED. IPUP
% IS SET TO -1 IF MITER .NE. 0, TO
% SIGNAL THAT THE JACOBIAN INVOLVED
% MAY NEED UPDATING LATER. THE LOCAL
% ERROR TEST IS MADE AND CONTROL
% PASSES TO STATEMENT 500 IF IT
% FAILS.
% ------------------------------------
if( ddebd1_31~=0 )
ddebd1_21 = -1;
end;
if( m==0 )
dsm = del./ddebd1_9(2,ddebd1_34);
end;
if( m>0 )
dsm = dvnrms(ddebd1_33,acor,ewt)./ddebd1_9(2,ddebd1_34);
end;
if( dsm>1.0d0 )
% ------------------------------------
% THE ERROR TEST FAILED. KFLAG KEEPS
% TRACK OF MULTIPLE FAILURES.
% RESTORE TN AND THE YH ARRAY TO
% THEIR PREVIOUS VALUES, AND PREPARE
% TO TRY THE STEP AGAIN. COMPUTE THE
% OPTIMUM STEP SIZE FOR THIS OR ONE
% LOWER ORDER. AFTER 2 OR MORE
% FAILURES, H IS FORCED TO DECREASE
% BY A FACTOR OF 0.2 OR LESS.
% ------------------------------------
ddebd1_28 = fix(ddebd1_28 - 1);
ddebd1_15 = told;
i1 = fix(ddebd1_24 + 1);
for jb = 1 : ddebd1_34;
i1 = fix(i1 - nyh);
for i = i1 : ddebd1_24;
yh1(i) = yh1(i) - yh1(i+nyh);
end; i = fix(ddebd1_24+1);
end; jb = fix(ddebd1_34+1);
ddebd1_8 = 2.0d0;
if( abs(ddebd1_11)<=ddebd1_12.*1.00001d0 )
% ---------------------------------
% ALL RETURNS ARE MADE THROUGH
% THIS SECTION. H IS SAVED IN
% HOLD TO ALLOW THE CALLER TO
% CHANGE H ON THE NEXT STEP.
% ---------------------------------
ddebd1_28 = -1;
% .................................EXIT
break;
% ...............EXIT
elseif( ddebd1_28>-3 ) ;
iredo = 2;
rhup = 0.0d0;
% ............EXIT
gt(9)=1;
continue;
% ---------------------------------------------------
% CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES
% HAVE OCCURRED. IF 10 FAILURES HAVE OCCURRED, EXIT
% WITH KFLAG = -1. IT IS ASSUMED THAT THE
% DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY
% HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
% DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO
% 1. THEN H IS REDUCED BY A FACTOR OF 10, AND THE
% STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES
% HMIN.
% ---------------------------------------------------
elseif( ddebd1_28~=-10 ) ;
rh = 0.1d0;
rh = max(ddebd1_12./abs(ddebd1_11),rh);
ddebd1_11 = ddebd1_11.*rh;
for i = 1 : ddebd1_33;
y(i) = yh(i,1);
end; i = fix(ddebd1_33+1);
[ddebd1_15,y,savf,rpar,ipar]=df(ddebd1_15,y,savf,rpar,ipar);
ddebd1_36 = fix(ddebd1_36 + 1);
for i = 1 : ddebd1_33;
yh(i,2) = ddebd1_11.*savf(i);
end; i = fix(ddebd1_33+1);
ddebd1_21 = fix(ddebd1_31);
ddebd1_20 = 5;
% ......EXIT
if( ddebd1_34==1 )
gt(4)=1;
continue;
end;
ddebd1_34 = 1;
ddebd1_29 = 2;
iret = 3;
continue;
else;
% ------------------------------------------------
% ALL RETURNS ARE MADE THROUGH THIS SECTION. H
% IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE
% H ON THE NEXT STEP.
% ------------------------------------------------
ddebd1_28 = -1;
% ..................EXIT
break;
end;
else;
% BEGIN BLOCK
% PERMITTING ...EXITS TO 360
% ------------------------------
% AFTER A SUCCESSFUL STEP,
% UPDATE THE YH ARRAY.
% CONSIDER CHANGING H IF IALTH
% = 1. OTHERWISE DECREASE
% IALTH BY 1. IF IALTH IS THEN
% 1 AND NQ .LT. MAXORD, THEN
% ACOR IS SAVED FOR USE IN A
% POSSIBLE ORDER INCREASE ON
% THE NEXT STEP. IF A CHANGE
% IN H IS CONSIDERED, AN
% INCREASE OR DECREASE IN ORDER
% BY ONE IS CONSIDERED ALSO. A
% CHANGE IN H IS MADE ONLY IF
% IT IS BY A FACTOR OF AT LEAST
% 1.1. IF NOT, IALTH IS SET TO
% 3 TO PREVENT TESTING FOR THAT
% MANY STEPS.
% ------------------------------
ddebd1_28 = 0;
iredo = 0;
ddebd1_35 = fix(ddebd1_35 + 1);
ddebd1_14 = ddebd1_11;
ddebd1_38 = fix(ddebd1_34);
for j = 1 : ddebd1_29;
for i = 1 : ddebd1_33;
yh(i,j) = yh(i,j) + ddebd1_4(j).*acor(i);
end; i = fix(ddebd1_33+1);
end; j = fix(ddebd1_29+1);
ddebd1_20 = fix(ddebd1_20 - 1);
if( ddebd1_20~=0 )
% ...EXIT
if( ddebd1_20<=1 )
% ...EXIT
if( ddebd1_29~=ddebd1_22 )
for i = 1 : ddebd1_33;
yh(i,ddebd1_22) = acor(i);
end; i = fix(ddebd1_33+1);
end;
end;
r = 1.0d0./ddebd1_9(2,ddebd1_38);
for i = 1 : ddebd1_33;
acor(i) = acor(i).*r;
end; i = fix(ddebd1_33+1);
% .................................EXIT
break;
else;
% ---------------------------
% REGARDLESS OF THE SUCCESS
% OR FAILURE OF THE STEP,
% FACTORS RHDN, RHSM, AND
% RHUP ARE COMPUTED, BY
% WHICH H COULD BE
% MULTIPLIED AT ORDER NQ -
% 1, ORDER NQ, OR ORDER NQ +
% 1, RESPECTIVELY. IN THE
% CASE OF FAILURE, RHUP =
% 0.0 TO AVOID AN ORDER
% INCREASE. THE LARGEST OF
% THESE IS DETERMINED AND
% THE NEW ORDER CHOSEN
% ACCORDINGLY. IF THE ORDER
% IS TO BE INCREASED, WE
% COMPUTE ONE ADDITIONAL
% SCALED DERIVATIVE.
% ---------------------------
rhup = 0.0d0;
% .....................EXIT
if( ddebd1_29~=ddebd1_22 )
for i = 1 : ddebd1_33;
savf(i) = acor(i) - yh(i,ddebd1_22);
end; i = fix(ddebd1_33+1);
dup = dvnrms(ddebd1_33,savf,ewt)./ddebd1_9(3,ddebd1_34);
exup = 1.0d0./(ddebd1_29+1);
% .....................EXIT
rhup = 1.0d0./(1.4d0.*dup.^exup+0.0000014d0);
end;
gt(9)=1;
continue;
end;
end;
end;
% ------------------------------------------
% THE CORRECTOR ITERATION FAILED TO
% CONVERGE IN 3 TRIES. IF MITER .NE. 0 AND
% THE JACOBIAN IS OUT OF DATE, DPJAC IS
% CALLED FOR THE NEXT TRY. OTHERWISE THE
% YH ARRAY IS RETRACTED TO ITS VALUES
% BEFORE PREDICTION, AND H IS REDUCED, IF
% POSSIBLE. IF H CANNOT BE REDUCED OR 10
% FAILURES HAVE OCCURRED, EXIT WITH KFLAG =
% -2.
% ------------------------------------------
% ...EXIT
end;
gt(7)=0;
if( ddebd1_21~=0 )
ddebd1_21 = fix(ddebd1_31);
gt(5)=1;
continue;
end;
end;
gt(8)=0;
ddebd1_15 = told;
ncf = fix(ncf + 1);
ddebd1_8 = 2.0d0;
i1 = fix(ddebd1_24 + 1);
for jb = 1 : ddebd1_34;
i1 = fix(i1 - nyh);
for i = i1 : ddebd1_24;
yh1(i) = yh1(i) - yh1(i+nyh);
end; i = fix(ddebd1_24+1);
end; jb = fix(ddebd1_34+1);
if( abs(ddebd1_11)<=ddebd1_12.*1.00001d0 )
ddebd1_28 = -2;
% ........................EXIT
break;
elseif( ncf~=10 ) ;
rh = 0.25d0;
ddebd1_21 = fix(ddebd1_31);
iredo = 1;
% .........EXIT
rh = max(rh,ddebd1_12./abs(ddebd1_11));
gt(3)=1;
continue;
else;
ddebd1_28 = -2;
% ........................EXIT
break;
end;
end;
gt(9)=0;
exsm = 1.0d0./ddebd1_29;
rhsm = 1.0d0./(1.2d0.*dsm.^exsm+0.0000012d0);
rhdn = 0.0d0;
if( ddebd1_34~=1 )
ddn = dvnrms(ddebd1_33,yh(sub2ind(size(yh),1,ddebd1_29):end),ewt)./ddebd1_9(1,ddebd1_34);
exdn = 1.0d0./ddebd1_34;
rhdn = 1.0d0./(1.3d0.*ddn.^exdn+0.0000013d0);
end;
if( rhsm>=rhup )
if( rhsm>=rhdn )
newq = fix(ddebd1_34);
rh = rhsm;
if( ddebd1_28==0 && rh<1.1d0 )
ddebd1_20 = 3;
r = 1.0d0./ddebd1_9(2,ddebd1_38);
for i = 1 : ddebd1_33;
acor(i) = acor(i).*r;
end; i = fix(ddebd1_33+1);
% .....................EXIT
break;
else;
if( ddebd1_28<=-2 )
rh = min(rh,0.2d0);
end;
% ------------------------------------------
% IF THERE IS A CHANGE OF ORDER, RESET NQ,
% L, AND THE COEFFICIENTS. IN ANY CASE H
% IS RESET ACCORDING TO RH AND THE YH ARRAY
% IS RESCALED. THEN EXIT FROM 680 IF THE
% STEP WAS OK, OR REDO THE STEP OTHERWISE.
% ------------------------------------------
% ............EXIT
if( newq==ddebd1_34 )
rh = max(rh,ddebd1_12./abs(ddebd1_11));
gt(3)=1;
continue;
else;
ddebd1_34 = fix(newq);
ddebd1_29 = fix(ddebd1_34 + 1);
iret = 2;
% ..................EXIT
continue;
end;
end;
end;
elseif( rhup>rhdn ) ;
newq = fix(ddebd1_29);
rh = rhup;
if( rh>=1.1d0 )
r = ddebd1_4(ddebd1_29)./ddebd1_29;
for i = 1 : ddebd1_33;
yh(i,newq+1) = acor(i).*r;
end; i = fix(ddebd1_33+1);
ddebd1_34 = fix(newq);
ddebd1_29 = fix(ddebd1_34 + 1);
iret = 2;
% ..................EXIT
continue;
else;
ddebd1_20 = 3;
r = 1.0d0./ddebd1_9(2,ddebd1_38);
for i = 1 : ddebd1_33;
acor(i) = acor(i).*r;
end; i = fix(ddebd1_33+1);
% ...........................EXIT
break;
end;
end;
newq = fix(ddebd1_34 - 1);
rh = rhdn;
if( ddebd1_28<0 && rh>1.0d0 )
rh = 1.0d0;
end;
if( ddebd1_28==0 && rh<1.1d0 )
ddebd1_20 = 3;
r = 1.0d0./ddebd1_9(2,ddebd1_38);
for i = 1 : ddebd1_33;
acor(i) = acor(i).*r;
% ..................EXIT
end; i = fix(ddebd1_33+1);
else;
if( ddebd1_28<=-2 )
rh = min(rh,0.2d0);
end;
% ---------------------------------------------
% IF THERE IS A CHANGE OF ORDER, RESET NQ, L,
% AND THE COEFFICIENTS. IN ANY CASE H IS
% RESET ACCORDING TO RH AND THE YH ARRAY IS
% RESCALED. THEN EXIT FROM 680 IF THE STEP
% WAS OK, OR REDO THE STEP OTHERWISE.
% ---------------------------------------------
% .........EXIT
if( newq==ddebd1_34 )
rh = max(rh,ddebd1_12./abs(ddebd1_11));
gt(3)=1;
continue;
else;
ddebd1_34 = fix(newq);
ddebd1_29 = fix(ddebd1_34 + 1);
iret = 2;
% ...............EXIT
continue;
end;
end;
break;
end;
ddebd1_6 = ddebd1_11;
ddebd1_27 = 1;
% ----------------------- end OF SUBROUTINE DSTOD
% -----------------------
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
yh1_shape=zeros(yh1_shape);yh1_shape(:)=yh1(1:numel(yh1_shape));yh1=yh1_shape;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_shape;
acor_shape=zeros(acor_shape);acor_shape(:)=acor(1:numel(acor_shape));acor=acor_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
end %subroutine dstod
%DECK DSTOR1
|
|
