| [neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,df,djac,rpar,ipar]=dpjac(neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,df,djac,rpar,ipar); |
function [neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,df,djac,rpar,ipar]=dpjac(neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,df,djac,rpar,ipar);
%***BEGIN PROLOGUE DPJAC
%***SUBSIDIARY
%***PURPOSE Subsidiary to DDEBDF
%***LIBRARY SLATEC
%***TYPE doubleprecision (PJAC-S, DPJAC-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% DPJAC sets up the iteration matrix (involving the Jacobian) for the
% integration package DDEBDF.
%
%***SEE ALSO DDEBDF
%***ROUTINES CALLED DGBFA, DGEFA, 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 DPJAC
%
persistent con di fac hl0 i i1 i2 ii j j1 jj lenp mba mband meb1 meband ml ml3 mu r r0 srur yi yj yjj ;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(i2), i2=0; end;
global ddebd1_12; if isempty(ddebd1_12), ddebd1_12=0; end;
if isempty(ii), ii=0; end;
global ddebd1_10; if isempty(ddebd1_10), ddebd1_10=zeros(1,14); end;
global ddebd1_11; if isempty(ddebd1_11), ddebd1_11=zeros(1,6); end;
if isempty(j), j=0; end;
if isempty(j1), j1=0; end;
if isempty(jj), jj=0; end;
global ddebd1_13; if isempty(ddebd1_13), ddebd1_13=0; end;
global ddebd1_14; if isempty(ddebd1_14), ddebd1_14=0; end;
global ddebd1_15; if isempty(ddebd1_15), ddebd1_15=0; end;
if isempty(lenp), lenp=0; end;
global ddebd1_18; if isempty(ddebd1_18), ddebd1_18=0; end;
if isempty(mba), mba=0; end;
if isempty(mband), mband=0; end;
if isempty(meb1), meb1=0; end;
if isempty(meband), meband=0; end;
global ddebd1_16; if isempty(ddebd1_16), ddebd1_16=0; end;
global ddebd1_17; if isempty(ddebd1_17), ddebd1_17=0; end;
if isempty(ml), ml=0; end;
if isempty(ml3), ml3=0; end;
if isempty(mu), mu=0; end;
global ddebd1_19; if isempty(ddebd1_19), ddebd1_19=0; end;
global ddebd1_22; if isempty(ddebd1_22), ddebd1_22=0; end;
global ddebd1_23; if isempty(ddebd1_23), ddebd1_23=0; end;
global ddebd1_20; if isempty(ddebd1_20), ddebd1_20=0; end;
global ddebd1_24; if isempty(ddebd1_24), ddebd1_24=0; end;
global ddebd1_21; if isempty(ddebd1_21), ddebd1_21=0; end;
if isempty(con), con=0; end;
if isempty(di), di=0; end;
global ddebd1_3; if isempty(ddebd1_3), ddebd1_3=0; end;
if isempty(fac), fac=0; end;
global ddebd1_4; if isempty(ddebd1_4), ddebd1_4=0; end;
if isempty(hl0), hl0=0; end;
global ddebd1_5; if isempty(ddebd1_5), ddebd1_5=0; end;
global ddebd1_6; if isempty(ddebd1_6), ddebd1_6=0; end;
global ddebd1_7; if isempty(ddebd1_7), ddebd1_7=0; end;
if isempty(r), r=0; end;
if isempty(r0), r0=0; end;
global ddebd1_1; if isempty(ddebd1_1), ddebd1_1=0; end;
global ddebd1_2; if isempty(ddebd1_2), ddebd1_2=zeros(1,210); end;
if isempty(srur), srur=0; end;
global ddebd1_8; if isempty(ddebd1_8), ddebd1_8=0; end;
global ddebd1_9; if isempty(ddebd1_9), ddebd1_9=0; end;
if isempty(yi), yi=0; end;
if isempty(yj), yj=0; end;
if isempty(yjj), yjj=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,[]);
ewt_shape=size(ewt);ewt=reshape(ewt,1,[]);
ftem_shape=size(ftem);ftem=reshape(ftem,1,[]);
savf_shape=size(savf);savf=reshape(savf,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 , rowns(210) , el0 , h , hmin , hmxi , hu ,tn , uround , iownd(14) , iowns(6) , ier ,jstart , kflag , l , meth , miter , maxord , n ,nq , nst , nfe , nje , nqu;
%% common /ddebd1/ ddebd1_1 , ddebd1_2(210) , ddebd1_3 , ddebd1_4 , ddebd1_5 , ddebd1_6 , ddebd1_7 ,ddebd1_8 , ddebd1_9 , ddebd1_10(14) , ddebd1_11(6) , ddebd1_12 ,ddebd1_13 , ddebd1_14 , ddebd1_15 , ddebd1_16 , ddebd1_17 , ddebd1_18 , ddebd1_19 ,ddebd1_20 , ddebd1_21 , ddebd1_22 , ddebd1_23 , ddebd1_24;
% ------------------------------------------------------------------
% DPJAC IS CALLED BY DSTOD TO COMPUTE AND PROCESS THE MATRIX
% P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
% HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE DJAC IF
% MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
% IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
% J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
% SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
% OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
% BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
%
% IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
% WITH DPJAC USES THE FOLLOWING..
% Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
% FTEM = WORK ARRAY OF LENGTH N (ACOR IN DSTOD ).
% SAVF = ARRAY CONTAINING DF EVALUATED AT PREDICTED Y.
% WM = doubleprecision WORK SPACE FOR MATRICES. ON OUTPUT IT
% contains THE
% INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU
% DECOMPOSITION OF P IF MITER IS 1, 2 , 4, OR 5.
% STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
% WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
% WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN
% INCREMENTS. WM(2) = H*EL0, SAVED FOR LATER USE IF MITER =
% 3.
% IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
% AT IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS
% THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER
% IS 4 OR 5.
% EL0 = EL(1) (INPUT).
% IER = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .NE. 0 IF
% P MATRIX FOUND TO BE SINGULAR.
% THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
% MITER, N, NFE, AND NJE.
%-----------------------------------------------------------------------
% BEGIN BLOCK PERMITTING ...EXITS TO 240
% BEGIN BLOCK PERMITTING ...EXITS TO 220
% BEGIN BLOCK PERMITTING ...EXITS TO 130
% BEGIN BLOCK PERMITTING ...EXITS TO 70
%***FIRST EXECUTABLE STATEMENT DPJAC
ddebd1_23 = fix(ddebd1_23 + 1);
hl0 = ddebd1_4.*ddebd1_3;
if( ddebd1_17==2 )
% IF MITER = 2, MAKE N CALLS TO DF TO APPROXIMATE J.
% --------------------
[fac ,ddebd1_19,savf,ewt]=dvnrms(ddebd1_19,savf,ewt);
r0 = 1000.0d0.*abs(ddebd1_4).*ddebd1_9.*ddebd1_19.*fac;
if( r0==0.0d0 )
r0 = 1.0d0;
end;
srur = wm(1);
j1 = 2;
for j = 1 : ddebd1_19;
yj = y(j);
r = max(srur.*abs(yj),r0.*ewt(j));
y(j) = y(j) + r;
fac = -hl0./r;
[ddebd1_8,y,ftem,rpar,ipar]=df(ddebd1_8,y,ftem,rpar,ipar);
for i = 1 : ddebd1_19;
wm(i+j1) =(ftem(i)-savf(i)).*fac;
end; i = fix(ddebd1_19+1);
y(j) = yj;
j1 = fix(j1 + ddebd1_19);
end; j = fix(ddebd1_19+1);
ddebd1_22 = fix(ddebd1_22 + ddebd1_19);
elseif( ddebd1_17==3 ) ;
% IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND
% P. ---------
wm(2) = hl0;
ddebd1_12 = 0;
r = ddebd1_3.*0.1d0;
for i = 1 : ddebd1_19;
y(i) = y(i) + r.*(ddebd1_4.*savf(i)-yh(i,2));
end; i = fix(ddebd1_19+1);
[ddebd1_8,y,wm(3:end),rpar,ipar]=df(ddebd1_8,y,wm(3:end),rpar,ipar);
ddebd1_22 = fix(ddebd1_22 + 1);
for i = 1 : ddebd1_19;
r0 = ddebd1_4.*savf(i) - yh(i,2);
di = 0.1d0.*r0 - ddebd1_4.*(wm(i+2)-savf(i));
wm(i+2) = 1.0d0;
if( abs(r0)>=ddebd1_9.*ewt(i) )
% .........EXIT
if( abs(di)==0.0d0 )
ddebd1_12 = -1;
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;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
ftem_shape=zeros(ftem_shape);ftem_shape(:)=ftem(1:numel(ftem_shape));ftem=ftem_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_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;
return;
else;
wm(i+2) = 0.1d0.*r0./di;
end;
% .........EXIT
end;
end; i = fix(ddebd1_19+1);
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;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
ftem_shape=zeros(ftem_shape);ftem_shape(:)=ftem(1:numel(ftem_shape));ftem=ftem_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_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;
return;
else;
while (1);
if( ddebd1_17==4 )
% IF MITER = 4, CALL DJAC AND MULTIPLY BY SCALAR.
% -----------------------
ml = fix(iwm(1));
mu = fix(iwm(2));
ml3 = 3;
mband = fix(ml + mu + 1);
meband = fix(mband + ml);
lenp = fix(meband.*ddebd1_19);
for i = 1 : lenp;
wm(i+2) = 0.0d0;
end; i = fix(lenp+1);
[ddebd1_8,y,wm(ml3:end),meband,rpar,ipar]=djac(ddebd1_8,y,wm(ml3:end),meband,rpar,ipar);
con = -hl0;
for i = 1 : lenp;
wm(i+2) = wm(i+2).*con;
% ...EXIT
end; i = fix(lenp+1);
elseif( ddebd1_17==5 ) ;
% IF MITER = 5, MAKE MBAND CALLS TO DF TO APPROXIMATE J.
% ----------------
ml = fix(iwm(1));
mu = fix(iwm(2));
mband = fix(ml + mu + 1);
mba = fix(min(mband,ddebd1_19));
meband = fix(mband + ml);
meb1 = fix(meband - 1);
srur = wm(1);
[fac ,ddebd1_19,savf,ewt]=dvnrms(ddebd1_19,savf,ewt);
r0 = 1000.0d0.*abs(ddebd1_4).*ddebd1_9.*ddebd1_19.*fac;
if( r0==0.0d0 )
r0 = 1.0d0;
end;
for j = 1 : mba;
for i = j : mband: ddebd1_19 ;
yi = y(i);
r = max(srur.*abs(yi),r0.*ewt(i));
y(i) = y(i) + r;
end; i = fix(ddebd1_19 +1);
[ddebd1_8,y,ftem,rpar,ipar]=df(ddebd1_8,y,ftem,rpar,ipar);
for jj = j : mband: ddebd1_19 ;
y(jj) = yh(jj,1);
yjj = y(jj);
r = max(srur.*abs(yjj),r0.*ewt(jj));
fac = -hl0./r;
i1 = fix(max(jj-mu,1));
i2 = fix(min(jj+ml,ddebd1_19));
ii = fix(jj.*meb1 - ml + 2);
for i = i1 : i2;
wm(ii+i) =(ftem(i)-savf(i)).*fac;
end; i = fix(i2+1);
end; jj = fix(ddebd1_19 +1);
end; j = fix(mba+1);
ddebd1_22 = fix(ddebd1_22 + mba);
else;
% IF MITER = 1, CALL DJAC AND MULTIPLY BY SCALAR.
% -----------------------
lenp = fix(ddebd1_19.*ddebd1_19);
for i = 1 : lenp;
wm(i+2) = 0.0d0;
end; i = fix(lenp+1);
[ddebd1_8,y,wm(3:end),ddebd1_19,rpar,ipar]=djac(ddebd1_8,y,wm(3:end),ddebd1_19,rpar,ipar);
con = -hl0;
for i = 1 : lenp;
wm(i+2) = wm(i+2).*con;
% ...EXIT
end; i = fix(lenp+1);
break;
end;
% ADD IDENTITY MATRIX.
% -------------------------------------------------
ii = fix(mband + 2);
for i = 1 : ddebd1_19;
wm(ii) = wm(ii) + 1.0d0;
ii = fix(ii + meband);
end; i = fix(ddebd1_19+1);
% DO LU DECOMPOSITION OF P.
% --------------------------------------------
[wm(sub2ind(size(wm),max(3,1)):end),meband,ddebd1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),ddebd1_12]=dgbfa(wm(sub2ind(size(wm),max(3,1)):end),meband,ddebd1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),ddebd1_12);
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;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
ftem_shape=zeros(ftem_shape);ftem_shape(:)=ftem(1:numel(ftem_shape));ftem=ftem_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_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;
return;
break;
end;
end;
% ADD IDENTITY MATRIX.
% -------------------------------------------------
j = 3;
for i = 1 : ddebd1_19;
wm(j) = wm(j) + 1.0d0;
j = fix(j +(ddebd1_19+1));
end; i = fix(ddebd1_19+1);
% DO LU DECOMPOSITION ON P.
% --------------------------------------------
ddebd1_19_orig=ddebd1_19; [wm(sub2ind(size(wm),max(3,1)):end),ddebd1_19,dumvar3,iwm(sub2ind(size(iwm),max(21,1)):end),ddebd1_12]=dgefa(wm(sub2ind(size(wm),max(3,1)):end),ddebd1_19,ddebd1_19,iwm(sub2ind(size(iwm),max(21,1)):end),ddebd1_12); ddebd1_19(dumvar3~=ddebd1_19_orig)=dumvar3(dumvar3~=ddebd1_19_orig);
% .........EXIT
% ----------------------- end OF SUBROUTINE DPJAC
% -----------------------
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;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
ftem_shape=zeros(ftem_shape);ftem_shape(:)=ftem(1:numel(ftem_shape));ftem=ftem_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_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 dpjac
%DECK DPLINT
|
|
