| [eps,f,fa,hmax,hold,impl,jtask,matdim,maxord,mint,miter,ml,mu,n,nde,save1,t,uround,users,y,ywt,h,mntold,mtrold,nfe,rc,yh,a,convrg,el,fac,ier,ipvt,nq,nwait,rh,rmax,save2,tq,trend,iswflg,jstate]=sdntl(eps,f,fa,hmax,hold,impl,jtask,matdim,maxord,mint,miter,m |
function [eps,f,fa,hmax,hold,impl,jtask,matdim,maxord,mint,miter,ml,mu,n,nde,save1,t,uround,users,y,ywt,h,mntold,mtrold,nfe,rc,yh,a,convrg,el,fac,ier,ipvt,nq,nwait,rh,rmax,save2,tq,trend,iswflg,jstate]=sdntl(eps,f,fa,hmax,hold,impl,jtask,matdim,maxord,mint,miter,ml,mu,n,nde,save1,t,uround,users,y,ywt,h,mntold,mtrold,nfe,rc,yh,a,convrg,el,fac,ier,ipvt,nq,nwait,rh,rmax,save2,tq,trend,iswflg,jstate);
%***BEGIN PROLOGUE SDNTL
%***SUBSIDIARY
%***PURPOSE Subroutine SDNTL is called to set parameters on the first
% call to SDSTP, on an internal restart, or when the user has
% altered MINT, MITER, and/or H.
%***LIBRARY SLATEC (SDRIVE)
%***TYPE SINGLE PRECISION (SDNTL-S, DDNTL-D, CDNTL-C)
%***AUTHOR Kahaner, D. K., (NIST)
% National Institute of Standards and Technology
% Gaithersburg, MD 20899
% Sutherland, C. D., (LANL)
% Mail Stop D466
% Los Alamos National Laboratory
% Los Alamos, NM 87545
%***DESCRIPTION
%
% On the first call, the order is set to 1 and the initial derivatives
% are calculated. RMAX is the maximum ratio by which H can be
% increased in one step. It is initially RMINIT to compensate
% for the small initial H, but then is normally equal to RMNORM.
% If a failure occurs (in corrector convergence or error test), RMAX
% is set at RMFAIL for the next increase.
% If the caller has changed MINT, or if JTASK = 0, SDCST is called
% to set the coefficients of the method. If the caller has changed H,
% YH must be rescaled. If H or MINT has been changed, NWAIT is
% reset to NQ + 2 to prevent further increases in H for that many
% steps. Also, RC is reset. RC is the ratio of new to old values of
% the coefficient L(0)*H. If the caller has changed MITER, RC is
% set to 0 to force the partials to be updated, if partials are used.
%***ROUTINES CALLED SDCST, SDSCL, SGBFA, SGBSL, SGEFA, SGESL, SNRM2
%***REVISION HISTORY (YYMMDD)
% 790601 DATE WRITTEN
% 900329 Initial submission to SLATEC.
%***end PROLOGUE SDNTL
persistent i iflag info oldl0 rminit summlv ;
if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(info), info=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([matdim])).*prod([matdim])-numel(a))],matdim,[]);
el_orig=el;el_shape=[13,12];el=reshape([el_orig(1:min(prod(el_shape),numel(el_orig))),zeros(1,max(0,prod(el_shape)-numel(el_orig)))],el_shape);
fac_shape=size(fac);fac=reshape(fac,1,[]);
if isempty(oldl0), oldl0=0; end;
save1_shape=size(save1);save1=reshape(save1,1,[]);
save2_shape=size(save2);save2=reshape(save2,1,[]);
if isempty(summlv), summlv=0; end;
tq_orig=tq;tq_shape=[3,12];tq=reshape([tq_orig(1:min(prod(tq_shape),numel(tq_orig))),zeros(1,max(0,prod(tq_shape)-numel(tq_orig)))],tq_shape);
y_shape=size(y);y=reshape(y,1,[]);
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([n])).*prod([n])-numel(yh))],n,[]);
ywt_shape=size(ywt);ywt=reshape(ywt,1,[]);
ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
if isempty(rminit), rminit=10000.0e0 ; end;
%***FIRST EXECUTABLE STATEMENT SDNTL
ier = false;
if( jtask>=0 )
if( jtask==0 )
[maxord,mint,iswflg,el,tq]=sdcst(maxord,mint,iswflg,el,tq);
rmax = rminit;
end;
rc = 0.0e0;
convrg = false;
trend = 1.0e0;
nq = 1;
nwait = 3;
[n,t,y,save2]=f(n,t,y,save2);
if( n==0 )
jstate = 6;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
nfe = fix(nfe + 1);
if( impl~=0 )
if( miter==3 )
iflag = 0;
[y,yh,ywt,save1,save2,t,h,el,impl,n,nde,iflag]=users(y,yh,ywt,save1,save2,t,h,el,impl,n,nde,iflag);
if( iflag==-1 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
if( n==0 )
jstate = 10;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
elseif( impl==1 ) ;
if( miter==1 || miter==2 )
[n,t,y,a,matdim,ml,mu,nde]=fa(n,t,y,a,matdim,ml,mu,nde);
if( n==0 )
jstate = 9;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,n,ipvt,info]=sgefa(a,matdim,n,ipvt,info);
if( info~=0 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,n,ipvt,save2]=sgesl(a,matdim,n,ipvt,save2,0);
elseif( miter==4 || miter==5 ) ;
[n,t,y,a(ml+1,1),matdim,ml,mu,nde]=fa(n,t,y,a(ml+1,1),matdim,ml,mu,nde);
if( n==0 )
jstate = 9;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,n,ml,mu,ipvt,info]=sgbfa(a,matdim,n,ml,mu,ipvt,info);
if( info~=0 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,n,ml,mu,ipvt,save2]=sgbsl(a,matdim,n,ml,mu,ipvt,save2,0);
end;
elseif( impl==2 ) ;
[n,t,y,a,matdim,ml,mu,nde]=fa(n,t,y,a,matdim,ml,mu,nde);
if( n==0 )
jstate = 9;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
for i = 1 : nde;
if( a(i,1)==0.0e0 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
else;
save2(i) = save2(i)./a(i,1);
end;
end; i = fix(nde+1);
for i = nde + 1 : n;
a(i,1) = 0.0e0;
end; i = fix(n+1);
elseif( impl==3 ) ;
if( miter==1 || miter==2 )
[n,t,y,a,matdim,ml,mu,nde]=fa(n,t,y,a,matdim,ml,mu,nde);
if( n==0 )
jstate = 9;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,nde,ipvt,info]=sgefa(a,matdim,nde,ipvt,info);
if( info~=0 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,nde,ipvt,save2]=sgesl(a,matdim,nde,ipvt,save2,0);
elseif( miter==4 || miter==5 ) ;
[n,t,y,a(ml+1,1),matdim,ml,mu,nde]=fa(n,t,y,a(ml+1,1),matdim,ml,mu,nde);
if( n==0 )
jstate = 9;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,nde,ml,mu,ipvt,info]=sgbfa(a,matdim,nde,ml,mu,ipvt,info);
if( info~=0 )
ier = true;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
return;
end;
[a,matdim,nde,ml,mu,ipvt,save2]=sgbsl(a,matdim,nde,ml,mu,ipvt,save2,0);
end;
end;
end;
for i = 1 : nde;
save1(i) = save2(i)./max(1.0e0,ywt(i));
end; i = fix(nde+1);
summlv = snrm2(nde,save1,1)./sqrt(real(nde));
if( summlv>eps./abs(h) )
h = (abs(eps./summlv).*sign(h));
end;
for i = 1 : n;
yh(i,2) = h.*save2(i);
end; i = fix(n+1);
if( miter==2 || miter==5 || iswflg==3 )
for i = 1 : n;
fac(i) = sqrt(uround);
end; i = fix(n+1);
end;
else;
if( miter~=mtrold )
mtrold = fix(miter);
rc = 0.0e0;
convrg = false;
end;
if( mint~=mntold )
mntold = fix(mint);
oldl0 = el(1,nq);
[maxord,mint,iswflg,el,tq]=sdcst(maxord,mint,iswflg,el,tq);
rc = rc.*el(1,nq)./oldl0;
nwait = fix(nq + 2);
end;
if( h~=hold )
nwait = fix(nq + 2);
rh = h./hold;
[hmax,n,nq,rmax,hold,rc,rh,yh]=sdscl(hmax,n,nq,rmax,hold,rc,rh,yh);
end;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
fac_shape=zeros(fac_shape);fac_shape(:)=fac(1:numel(fac_shape));fac=fac_shape;
save1_shape=zeros(save1_shape);save1_shape(:)=save1(1:numel(save1_shape));save1=save1_shape;
save2_shape=zeros(save2_shape);save2_shape(:)=save2(1:numel(save2_shape));save2=save2_shape;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
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;
ywt_shape=zeros(ywt_shape);ywt_shape(:)=ywt(1:numel(ywt_shape));ywt=ywt_shape;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
end
%DECK SDNTP
|
|
