| [neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,f,jac,rpar,ipar]=pjac(neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,f,jac,rpar,ipar); |
function [neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,f,jac,rpar,ipar]=pjac(neq,y,yh,nyh,ewt,ftem,savf,wm,iwm,f,jac,rpar,ipar);
%***BEGIN PROLOGUE PJAC
%***SUBSIDIARY
%***PURPOSE Subsidiary to DEBDF
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (PJAC-S, DPJAC-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% PJAC sets up the iteration matrix (involving the Jacobian) for the
% integration package DEBDF.
%
%***SEE ALSO DEBDF
%***ROUTINES CALLED SGBFA, SGEFA, VNWRMS
%***COMMON BLOCKS DEBDF1
%***REVISION HISTORY (YYMMDD)
% 800901 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (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 PJAC
%
%LLL. OPTIMIZE
persistent con di fac hl0 i i1 i2 igo 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 debdf1_12; if isempty(debdf1_12), debdf1_12=0; end;
if isempty(ii), ii=0; end;
global debdf1_10; if isempty(debdf1_10), debdf1_10=zeros(1,14); end;
global debdf1_11; if isempty(debdf1_11), debdf1_11=zeros(1,6); end;
if isempty(j), j=0; end;
if isempty(j1), j1=0; end;
if isempty(jj), jj=0; end;
global debdf1_13; if isempty(debdf1_13), debdf1_13=0; end;
global debdf1_14; if isempty(debdf1_14), debdf1_14=0; end;
global debdf1_15; if isempty(debdf1_15), debdf1_15=0; end;
if isempty(lenp), lenp=0; end;
global debdf1_18; if isempty(debdf1_18), debdf1_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 debdf1_16; if isempty(debdf1_16), debdf1_16=0; end;
global debdf1_17; if isempty(debdf1_17), debdf1_17=0; end;
if isempty(ml), ml=0; end;
if isempty(ml3), ml3=0; end;
if isempty(mu), mu=0; end;
global debdf1_19; if isempty(debdf1_19), debdf1_19=0; end;
global debdf1_22; if isempty(debdf1_22), debdf1_22=0; end;
global debdf1_23; if isempty(debdf1_23), debdf1_23=0; end;
global debdf1_20; if isempty(debdf1_20), debdf1_20=0; end;
global debdf1_24; if isempty(debdf1_24), debdf1_24=0; end;
global debdf1_21; if isempty(debdf1_21), debdf1_21=0; end;
if isempty(igo), igo=0; end;
global debdf1_1; if isempty(debdf1_1), debdf1_1=0; end;
global debdf1_2; if isempty(debdf1_2), debdf1_2=zeros(1,210); end;
global debdf1_3; if isempty(debdf1_3), debdf1_3=0; end;
global debdf1_4; if isempty(debdf1_4), debdf1_4=0; end;
global debdf1_5; if isempty(debdf1_5), debdf1_5=0; end;
global debdf1_6; if isempty(debdf1_6), debdf1_6=0; end;
global debdf1_7; if isempty(debdf1_7), debdf1_7=0; end;
global debdf1_8; if isempty(debdf1_8), debdf1_8=0; end;
global debdf1_9; if isempty(debdf1_9), debdf1_9=0; end;
if isempty(con), con=0; end;
if isempty(di), di=0; end;
if isempty(fac), fac=0; end;
if isempty(hl0), hl0=0; end;
if isempty(r), r=0; end;
if isempty(r0), r0=0; end;
if isempty(srur), srur=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 /debdf1/ 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 /debdf1/ debdf1_1 , debdf1_2(210) , debdf1_3 , debdf1_4 , debdf1_5 , debdf1_6 , debdf1_7 ,debdf1_8 , debdf1_9 , debdf1_10(14) , debdf1_11(6) , debdf1_12 ,debdf1_13 , debdf1_14 , debdf1_15 , debdf1_16 , debdf1_17 , debdf1_18 , debdf1_19 ,debdf1_20 , debdf1_21 , debdf1_22 , debdf1_23 , debdf1_24;
%-----------------------------------------------------------------------
% PJAC IS CALLED BY STOD 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 JAC 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 SGEFA IF MITER = 1 OR 2, AND BY SGBFA IF MITER = 4 OR 5.
%
% IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
% WITH PJAC USES THE FOLLOWING..
% Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
% FTEM = WORK ARRAY OF LENGTH N (ACOR IN STOD ).
% SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
% WM = REAL 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.
%-----------------------------------------------------------------------
%***FIRST EXECUTABLE STATEMENT PJAC
debdf1_23 = fix(debdf1_23 + 1);
hl0 = debdf1_4.*debdf1_3;
igo=0;
if( debdf1_17==2 )
% IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
[fac ,debdf1_19,savf,ewt]=vnwrms(debdf1_19,savf,ewt);
r0 = 1000.0e0.*abs(debdf1_4).*debdf1_9.*debdf1_19.*fac;
if( r0==0.0e0 )
r0 = 1.0e0;
end;
srur = wm(1);
j1 = 2;
for j = 1 : debdf1_19;
yj = y(j);
r = max(srur.*abs(yj),r0.*ewt(j));
y(j) = y(j) + r;
fac = -hl0./r;
[debdf1_8,y,ftem,rpar,ipar]=f(debdf1_8,y,ftem,rpar,ipar);
for i = 1 : debdf1_19;
wm(i+j1) =(ftem(i)-savf(i)).*fac;
end; i = fix(debdf1_19+1);
y(j) = yj;
j1 = fix(j1 + debdf1_19);
end; j = fix(debdf1_19+1);
debdf1_22 = fix(debdf1_22 + debdf1_19);
elseif( debdf1_17==3 ) ;
% IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
wm(2) = hl0;
debdf1_12 = 0;
r = debdf1_3.*0.1e0;
for i = 1 : debdf1_19;
y(i) = y(i) + r.*(debdf1_4.*savf(i)-yh(i,2));
end; i = fix(debdf1_19+1);
[debdf1_8,y,wm(3:end),rpar,ipar]=f(debdf1_8,y,wm(3:end),rpar,ipar);
debdf1_22 = fix(debdf1_22 + 1);
for i = 1 : debdf1_19;
r0 = debdf1_4.*savf(i) - yh(i,2);
di = 0.1e0.*r0 - debdf1_4.*(wm(i+2)-savf(i));
wm(i+2) = 1.0e0;
if( abs(r0)>=debdf1_9.*ewt(i) )
if( abs(di)==0.0e0 )
debdf1_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.1e0.*r0./di;
end;
end;
end; i = fix(debdf1_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;
if( debdf1_17==4 )
% IF MITER = 4, CALL JAC 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.*debdf1_19);
for i = 1 : lenp;
wm(i+2) = 0.0e0;
end; i = fix(lenp+1);
[debdf1_8,y,wm(ml3:end),meband,rpar,ipar]=jac(debdf1_8,y,wm(ml3:end),meband,rpar,ipar);
con = -hl0;
for i = 1 : lenp;
wm(i+2) = wm(i+2).*con;
end; i = fix(lenp+1);
elseif( debdf1_17==5 ) ;
% IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
ml = fix(iwm(1));
mu = fix(iwm(2));
mband = fix(ml + mu + 1);
mba = fix(min(mband,debdf1_19));
meband = fix(mband + ml);
meb1 = fix(meband - 1);
srur = wm(1);
[fac ,debdf1_19,savf,ewt]=vnwrms(debdf1_19,savf,ewt);
r0 = 1000.0e0.*abs(debdf1_4).*debdf1_9.*debdf1_19.*fac;
if( r0==0.0e0 )
r0 = 1.0e0;
end;
for j = 1 : mba;
for i = j : mband: debdf1_19 ;
yi = y(i);
r = max(srur.*abs(yi),r0.*ewt(i));
y(i) = y(i) + r;
end; i = fix(debdf1_19 +1);
[debdf1_8,y,ftem,rpar,ipar]=f(debdf1_8,y,ftem,rpar,ipar);
for jj = j : mband: debdf1_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,debdf1_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(debdf1_19 +1);
end; j = fix(mba+1);
debdf1_22 = fix(debdf1_22 + mba);
else;
% IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
lenp = fix(debdf1_19.*debdf1_19);
for i = 1 : lenp;
wm(i+2) = 0.0e0;
end; i = fix(lenp+1);
[debdf1_8,y,wm(3:end),debdf1_19,rpar,ipar]=jac(debdf1_8,y,wm(3:end),debdf1_19,rpar,ipar);
con = -hl0;
for i = 1 : lenp;
wm(i+2) = wm(i+2).*con;
end; i = fix(lenp+1);
igo=1;
end;
if(igo==0)
% ADD IDENTITY MATRIX. -------------------------------------------------
ii = fix(mband + 2);
for i = 1 : debdf1_19;
wm(ii) = wm(ii) + 1.0e0;
ii = fix(ii + meband);
end; i = fix(debdf1_19+1);
% DO LU DECOMPOSITION OF P. --------------------------------------------
[wm(sub2ind(size(wm),max(3,1)):end),meband,debdf1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),debdf1_12]=sgbfa(wm(sub2ind(size(wm),max(3,1)):end),meband,debdf1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),debdf1_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;
end;
end;
% ADD IDENTITY MATRIX. -------------------------------------------------
j = 3;
for i = 1 : debdf1_19;
wm(j) = wm(j) + 1.0e0;
j = fix(j +(debdf1_19+1));
end; i = fix(debdf1_19+1);
% DO LU DECOMPOSITION ON P. --------------------------------------------
debdf1_19_orig=debdf1_19; [wm(sub2ind(size(wm),max(3,1)):end),debdf1_19,dumvar3,iwm(sub2ind(size(iwm),max(21,1)):end),debdf1_12]=sgefa(wm(sub2ind(size(wm),max(3,1)):end),debdf1_19,debdf1_19,iwm(sub2ind(size(iwm),max(21,1)):end),debdf1_12); debdf1_19(dumvar3~=debdf1_19_orig)=dumvar3(dumvar3~=debdf1_19_orig);
%----------------------- end OF SUBROUTINE PJAC -----------------------
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 pjac
%DECK PNNZRS
|
|
