[dfdy,el,fa,h,ierror,impl,ipvt,matdim,miter,ml,mu,n,nde,nq,t,users,y,yh,ywt,evalfa,save1,save2,a,d,jstate]=ddcor(dfdy,el,fa,h,ierror,impl,ipvt,matdim,miter,ml,mu,n,nde,nq,t,users,y,yh,ywt,evalfa,save1,save2,a,d,jstate);

function [dfdy,el,fa,h,ierror,impl,ipvt,matdim,miter,ml,mu,n,nde,nq,t,users,y,yh,ywt,evalfa,save1,save2,a,d,jstate]=ddcor(dfdy,el,fa,h,ierror,impl,ipvt,matdim,miter,ml,mu,n,nde,nq,t,users,y,yh,ywt,evalfa,save1,save2,a,d,jstate);
%***BEGIN PROLOGUE  DDCOR
%***SUBSIDIARY
%***PURPOSE  Subroutine DDCOR computes corrections to the Y array.
%***LIBRARY   SLATEC (SDRIVE)
%***TYPE      doubleprecision (SDCOR-S, DDCOR-D, CDCOR-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
%
%  In the case of functional iteration, update Y directly from the
%  result of the last call to F.
%  In the case of the chord method, compute the corrector error and
%  solve the linear system with that as right hand side and DFDY as
%  coefficient matrix, using the LU decomposition if MITER is 1, 2, 4,
%  or 5.
%
%***ROUTINES CALLED  DGBSL, DGESL, DNRM2
%***REVISION HISTORY  (YYMMDD)
%   790601  DATE WRITTEN
%   900329  Initial submission to SLATEC.
%***end PROLOGUE  DDCOR
persistent i iflag j mw ; 

if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(j), j=0; end;
if isempty(mw), mw=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([matdim])).*prod([matdim])-numel(a))],matdim,[]);
dfdy_shape=size(dfdy);dfdy=reshape([dfdy(:).',zeros(1,ceil(numel(dfdy)./prod([matdim])).*prod([matdim])-numel(dfdy))],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);
save1_shape=size(save1);save1=reshape(save1,1,[]);
save2_shape=size(save2);save2=reshape(save2,1,[]);
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,[]);
%***FIRST EXECUTABLE STATEMENT  DDCOR
if( miter==0 )
if( ierror==1 || ierror==5 )
for i = 1 : n;
save1(i) =(h.*save2(i)-yh(i,2)-save1(i))./ywt(i);
end; i = fix(n+1);
else;
for i = 1 : n;
save1(i) =(h.*save2(i)-yh(i,2)-save1(i))./max(abs(y(i)),ywt(i));
end; i = fix(n+1);
end;
d = dnrm2(n,save1,1)./sqrt((n));
for i = 1 : n;
save1(i) = h.*save2(i) - yh(i,2);
end; i = fix(n+1);
elseif( miter==1 || miter==2 ) ;
if( impl==0 )
for i = 1 : n;
save2(i) = h.*save2(i) - yh(i,2) - save1(i);
end; i = fix(n+1);
elseif( impl==1 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i);
end; i = fix(n+1);
for j = 1 : n;
for i = 1 : n;
save2(i) = save2(i) - a(i,j).*(yh(j,2)+save1(j));
end; i = fix(n+1);
end; j = fix(n+1);
elseif( impl==2 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i) - a(i,1).*(yh(i,2)+save1(i));
end; i = fix(n+1);
elseif( impl==3 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i);
end; i = fix(n+1);
for j = 1 : nde;
for i = 1 : nde;
save2(i) = save2(i) - a(i,j).*(yh(j,2)+save1(j));
end; i = fix(nde+1);
end; j = fix(nde+1);
end;
[dfdy,matdim,n,ipvt,save2]=dgesl(dfdy,matdim,n,ipvt,save2,0);
if( ierror==1 || ierror==5 )
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./ywt(i);
end; i = fix(n+1);
else;
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./max(abs(y(i)),ywt(i));
end; i = fix(n+1);
end;
d = dnrm2(n,save2,1)./sqrt((n));
elseif( miter==4 || miter==5 ) ;
if( impl==0 )
for i = 1 : n;
save2(i) = h.*save2(i) - yh(i,2) - save1(i);
end; i = fix(n+1);
elseif( impl==1 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i);
end; i = fix(n+1);
mw = fix(ml + 1 + mu);
for j = 1 : n;
for i = max(ml+1,mw+1-j) : min(mw+n-j,mw+ml);
save2(i+j-mw) = save2(i+j-mw) - a(i,j).*(yh(j,2)+save1(j));
end; i = fix(min(mw+n-j,mw+ml)+1);
end; j = fix(n+1);
elseif( impl==2 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i) - a(i,1).*(yh(i,2)+save1(i));
end; i = fix(n+1);
elseif( impl==3 ) ;
if( evalfa )
[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;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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;
else;
evalfa = true;
end;
for i = 1 : n;
save2(i) = h.*save2(i);
end; i = fix(n+1);
mw = fix(ml + 1 + mu);
for j = 1 : nde;
for i = max(ml+1,mw+1-j) : min(mw+nde-j,mw+ml);
save2(i+j-mw) = save2(i+j-mw) - a(i,j).*(yh(j,2)+save1(j));
end; i = fix(min(mw+nde-j,mw+ml)+1);
end; j = fix(nde+1);
end;
[dfdy,matdim,n,ml,mu,ipvt,save2]=dgbsl(dfdy,matdim,n,ml,mu,ipvt,save2,0);
if( ierror==1 || ierror==5 )
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./ywt(i);
end; i = fix(n+1);
else;
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./max(abs(y(i)),ywt(i));
end; i = fix(n+1);
end;
d = dnrm2(n,save2,1)./sqrt((n));
elseif( miter==3 ) ;
iflag = 2;
[y,yh(1,2),ywt,save1,save2,t,h,el(1,nq),impl,n,nde,iflag]=users(y,yh(1,2),ywt,save1,save2,t,h,el(1,nq),impl,n,nde,iflag);
if( n==0 )
jstate = 10;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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( ierror==1 || ierror==5 )
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./ywt(i);
end; i = fix(n+1);
else;
for i = 1 : n;
save1(i) = save1(i) + save2(i);
save2(i) = save2(i)./max(abs(y(i)),ywt(i));
end; i = fix(n+1);
end;
d = dnrm2(n,save2,1)./sqrt((n));
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
dfdy_shape=zeros(dfdy_shape);dfdy_shape(:)=dfdy(1:numel(dfdy_shape));dfdy=dfdy_shape;
el_orig(1:prod(el_shape))=el;el=el_orig;
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;
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 DDCST