| [ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc]=drkfab(ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc); |
function [ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc]=drkfab(ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc);
%***BEGIN PROLOGUE DRKFAB
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBVSUP
%***LIBRARY SLATEC
%***TYPE doubleprecision (RKFAB-S, DRKFAB-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
%
% subroutine DRKFAB integrates the initial value equations using
% the variable-step Runge-Kutta-Fehlberg integration scheme or
% the variable-order Adams method and orthonormalization
% determined by a linear dependence test.
%
% **********************************************************************
%
%***SEE ALSO DBVSUP
%***ROUTINES CALLED DBVDER, DDEABM, DDERKF, DREORT, DSTOR1
%***COMMON BLOCKS DML15T, DML17B, DML18J, DML8SZ
%***REVISION HISTORY (YYMMDD)
% 750601 DATE WRITTEN
% 890921 Realigned order of variables in certain COMMON blocks.
% (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DRKFAB
%
persistent gt idid ipar j jflag jon kod kopp nfcp1 non xxop ;
global dml18j_18; if isempty(dml18j_18), dml18j_18=0; end;
if isempty(idid), idid=0; end;
global dml8sz_3; if isempty(dml8sz_3), dml8sz_3=0; end;
global dml18j_11; if isempty(dml18j_11), dml18j_11=0; end;
global dml15t_9; if isempty(dml15t_9), dml15t_9=zeros(1,15); end;
global dml8sz_4; if isempty(dml8sz_4), dml8sz_4=0; end;
global dml18j_12; if isempty(dml18j_12), dml18j_12=0; end;
if isempty(ipar), ipar=0; end;
global dml15t_10; if isempty(dml15t_10), dml15t_10=0; end;
global dml8sz_5; if isempty(dml8sz_5), dml8sz_5=0; end;
if isempty(j), j=0; end;
if isempty(jflag), jflag=0; end;
if isempty(jon), jon=0; end;
global dml17b_4; if isempty(dml17b_4), dml17b_4=0; end;
global dml17b_13; if isempty(dml17b_13), dml17b_13=0; end;
global dml17b_14; if isempty(dml17b_14), dml17b_14=0; end;
global dml17b_5; if isempty(dml17b_5), dml17b_5=0; end;
global dml17b_6; if isempty(dml17b_6), dml17b_6=0; end;
global dml17b_7; if isempty(dml17b_7), dml17b_7=0; end;
global dml17b_8; if isempty(dml17b_8), dml17b_8=0; end;
global dml17b_9; if isempty(dml17b_9), dml17b_9=0; end;
global dml17b_10; if isempty(dml17b_10), dml17b_10=0; end;
global dml17b_11; if isempty(dml17b_11), dml17b_11=0; end;
global dml17b_12; if isempty(dml17b_12), dml17b_12=0; end;
global dml17b_17; if isempty(dml17b_17), dml17b_17=0; end;
global dml17b_1; if isempty(dml17b_1), dml17b_1=0; end;
global dml15t_11; if isempty(dml15t_11), dml15t_11=0; end;
if isempty(kod), kod=0; end;
global dml15t_12; if isempty(dml15t_12), dml15t_12=0; end;
if isempty(kopp), kopp=0; end;
global dml17b_15; if isempty(dml17b_15), dml17b_15=0; end;
global dml17b_16; if isempty(dml17b_16), dml17b_16=0; end;
global dml17b_18; if isempty(dml17b_18), dml17b_18=0; end;
global dml15t_13; if isempty(dml15t_13), dml15t_13=0; end;
global dml15t_14; if isempty(dml15t_14), dml15t_14=0; end;
global dml18j_7; if isempty(dml18j_7), dml18j_7=0; end;
global dml8sz_6; if isempty(dml8sz_6), dml8sz_6=0; end;
global dml18j_8; if isempty(dml18j_8), dml18j_8=0; end;
global dml17b_3; if isempty(dml17b_3), dml17b_3=0; end;
global dml17b_2; if isempty(dml17b_2), dml17b_2=0; end;
global dml18j_10; if isempty(dml18j_10), dml18j_10=0; end;
global dml18j_15; if isempty(dml18j_15), dml18j_15=0; end;
global dml18j_17; if isempty(dml18j_17), dml18j_17=0; end;
global dml8sz_7; if isempty(dml8sz_7), dml8sz_7=0; end;
if isempty(nfcp1), nfcp1=0; end;
global dml18j_5; if isempty(dml18j_5), dml18j_5=0; end;
if isempty(non), non=0; end;
global dml18j_6; if isempty(dml18j_6), dml18j_6=0; end;
global dml18j_13; if isempty(dml18j_13), dml18j_13=0; end;
global dml15t_15; if isempty(dml15t_15), dml15t_15=0; end;
global dml18j_9; if isempty(dml18j_9), dml18j_9=0; end;
global dml18j_14; if isempty(dml18j_14), dml18j_14=0; end;
global dml18j_16; if isempty(dml18j_16), dml18j_16=0; end;
global dml18j_4; if isempty(dml18j_4), dml18j_4=0; end;
ip_shape=size(ip);ip=reshape([ip(:).',zeros(1,ceil(numel(ip)./prod([nfcc])).*prod([nfcc])-numel(ip))],nfcc,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
if isempty(gt), gt=zeros(1,3); end;
global dml18j_1; if isempty(dml18j_1), dml18j_1=0; end;
global dml8sz_1; if isempty(dml8sz_1), dml8sz_1=0; end;
g_shape=size(g);g=reshape(g,1,[]);
p_shape=size(p);p=reshape([p(:).',zeros(1,ceil(numel(p)./prod([ntp])).*prod([ntp])-numel(p))],ntp,[]);
global dml15t_2; if isempty(dml15t_2), dml15t_2=0; end;
global dml15t_1; if isempty(dml15t_1), dml15t_1=0; end;
global dml18j_2; if isempty(dml18j_2), dml18j_2=0; end;
s_shape=size(s);s=reshape(s,1,[]);
stowa_shape=size(stowa);stowa=reshape(stowa,1,[]);
global dml15t_3; if isempty(dml15t_3), dml15t_3=0; end;
global dml18j_3; if isempty(dml18j_3), dml18j_3=0; end;
u_shape=size(u);u=reshape([u(:).',zeros(1,ceil(numel(u)./prod([ncomp,nfc])).*prod([ncomp,nfc])-numel(u))],ncomp,nfc,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([ncomp])).*prod([ncomp])-numel(v))],ncomp,[]);
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([nfcc])).*prod([nfcc])-numel(w))],nfcc,[]);
work_shape=size(work);work=reshape(work,1,[]);
global dml15t_4; if isempty(dml15t_4), dml15t_4=0; end;
global dml15t_5; if isempty(dml15t_5), dml15t_5=0; end;
global dml15t_6; if isempty(dml15t_6), dml15t_6=0; end;
global dml15t_8; if isempty(dml15t_8), dml15t_8=0; end;
global dml15t_7; if isempty(dml15t_7), dml15t_7=0; end;
xpts_shape=size(xpts);xpts=reshape(xpts,1,[]);
global dml8sz_2; if isempty(dml8sz_2), dml8sz_2=0; end;
if isempty(xxop), xxop=0; end;
yhp_shape=size(yhp);yhp=reshape([yhp(:).',zeros(1,ceil(numel(yhp)./prod([ncomp])).*prod([ncomp])-numel(yhp))],ncomp,[]);
z_shape=size(z);z=reshape(z,1,[]);
%
% ******************************************************************
%
% common :: ;
%% common /dml8sz/ c , xsav , igofx , inhomo , ivp , ncompd , nfcd;
%% common /dml8sz/ dml8sz_1 , dml8sz_2 , dml8sz_3 , dml8sz_4 , dml8sz_5 , dml8sz_6 , dml8sz_7;
% common :: ;
%% common /dml15t/ px , pwcnd , tnd , x , xbeg , xend , xot , xop ,info(15) , istkop , knswot , kop , lotjp ,mnswot , nswot;
%% common /dml15t/ dml15t_1 , dml15t_2 , dml15t_3 , dml15t_4 , dml15t_5 , dml15t_6 , dml15t_7 , dml15t_8 ,dml15t_9(15) , dml15t_10 , dml15t_11 , dml15t_12 , dml15t_13 ,dml15t_14 , dml15t_15;
% common :: ;
%% common /dml18j/ ae , re , tol , nxptsd , nic , nopg , mxnond ,ndisk , ntape , neq , indpvt , integ , nps ,ntpd , neqivp , numort , nfccd , icoco;
%% common /dml18j/ dml18j_1 , dml18j_2 , dml18j_3 , dml18j_4 , dml18j_5 , dml18j_6 , dml18j_7 ,dml18j_8 , dml18j_9 , dml18j_10 , dml18j_11 , dml18j_12 , dml18j_13 ,dml18j_14 , dml18j_15 , dml18j_16 , dml18j_17 , dml18j_18;
% common :: ;
%% common /dml17b/ kkkzpw , needw , neediw , k1 , k2 , k3 , k4 , k5 ,k6 , k7 , k8 , k9 , k10 , k11 , l1 , l2 , kkkint ,lllint;
%% common /dml17b/ dml17b_1 , dml17b_2 , dml17b_3 , dml17b_4 , dml17b_5 , dml17b_6 , dml17b_7 , dml17b_8 ,dml17b_9 , dml17b_10 , dml17b_11 , dml17b_12 , dml17b_13 , dml17b_14 , dml17b_15 , dml17b_16 , dml17b_17 ,dml17b_18;
%
%
% *****************************************************************
% INITIALIZATION OF COUNTERS AND VARIABLES.
%
% BEGIN BLOCK PERMITTING ...EXITS TO 220
% BEGIN BLOCK PERMITTING ...EXITS TO 10
%***FIRST EXECUTABLE STATEMENT DRKFAB
kod = 1;
non = 1;
dml15t_4 = dml15t_5;
jon = 1;
dml15t_9(1) = 0;
dml15t_9(2) = 0;
dml15t_9(3) = 1;
dml15t_9(4) = 1;
work(1) = dml15t_6;
% ...EXIT
if( dml18j_6~=0 )
dml15t_9(3) = 0;
if( dml15t_4==z(1) )
jon = 2;
end;
end;
nfcp1 = fix(nfc + 1);
%
% ***************************************************************
% *****BEGINNING OF INTEGRATION LOOP AT OUTPUT
% POINTS.******************
% ***************************************************************
%
gt(:)=0;
for kopp = 2 : nxpts;
while (1);
if(gt(2)==0)
if(gt(1)==0)
dml15t_12 = fix(kopp);
dml15t_8 = xpts(dml15t_12);
if( dml18j_8==0 )
kod = fix(dml15t_12);
end;
%
%
% STEP BY STEP INTEGRATION LOOP BETWEEN OUTPUT POINTS.
%
% BEGIN BLOCK PERMITTING ...EXITS TO 190
% BEGIN BLOCK PERMITTING ...EXITS TO 30
end;
gt(1)=0;
xxop = dml15t_8;
% ...EXIT
if( dml18j_6~=0 )
if( dml15t_6>dml15t_5 && dml15t_8>z(jon) )
xxop = z(jon);
end;
if( dml15t_6<dml15t_5 && dml15t_8<z(jon) )
xxop = z(jon);
end;
end;
%
% ******************************************************
% BEGIN BLOCK PERMITTING ...EXITS TO 170
end;
gt(2)=0;
if( dml18j_12==2 )
% DDEABM INTEGRATOR
%
[dumvar1,dml18j_10,dml15t_4,yhp,xxop,dml15t_9,dml18j_2,dml18j_1,idid,work,dml17b_17,iwork,dml17b_18,g,ipar]=ddeabm(@dbvder,dml18j_10,dml15t_4,yhp,xxop,dml15t_9,dml18j_2,dml18j_1,idid,work,dml17b_17,iwork,dml17b_18,g,ipar);
else;
% DDERKF INTEGRATOR
%
[dumvar1,dml18j_10,dml15t_4,yhp,xxop,dml15t_9,dml18j_2,dml18j_1,idid,work,dml17b_17,iwork,dml17b_18,g,ipar]=dderkf(@dbvder,dml18j_10,dml15t_4,yhp,xxop,dml15t_9,dml18j_2,dml18j_1,idid,work,dml17b_17,iwork,dml17b_18,g,ipar);
end;
if( idid>=1 )
%
% ************************************************
% GRAM-SCHMIDT ORTHOGONALIZATION TEST FOR
% ORTHONORMALIZATION (TEMPORARILY USING U AND
% V IN THE TEST)
%
if( dml18j_6==0 )
jflag = 1;
if( dml8sz_4==3 && dml15t_4==dml15t_6 )
jflag = 3;
end;
elseif( xxop==z(jon) ) ;
jflag = 2;
else;
%
% ******************************************
% CONTINUE INTEGRATION IF WE ARE NOT AT
% AN OUTPUT POINT.
%
% ..................EXIT
% .........EXIT
if( idid==1 )
gt(2)=1;
continue;
end;
break;
end;
%
if( dml18j_8==0 )
non = fix(dml18j_16 + 1);
end;
[ncomp,u(sub2ind(size(u),1,1,kod):end),v(sub2ind(size(v),1,kod):end),yhp,niv,w(sub2ind(size(w),1,non):end),s,p(sub2ind(size(p),1,non):end),ip(sub2ind(size(p),1,non):end),stowa,jflag]=dreort(ncomp,u(sub2ind(size(u),1,1,kod):end),v(sub2ind(size(v),1,kod):end),yhp,niv,w(sub2ind(size(w),1,non):end),s,p(sub2ind(size(p),1,non):end),ip(sub2ind(size(p),1,non):end),stowa,jflag);
%
if( jflag==30 )
iflag = 30;
% .....................EXIT
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
%
elseif( jflag==10 ) ;
dml15t_8 = xpts(dml15t_12);
if( dml18j_8==0 )
kod = fix(dml15t_12);
end;
% ............EXIT
gt(1)=1;
continue;
%
elseif( jflag==0 ) ;
%
% ************************************************
% STORE ORTHONORMALIZED VECTORS INTO SOLUTION
% VECTORS.
%
if( dml18j_16>=mxnon )
if( dml15t_4~=dml15t_6 )
iflag = 13;
% .....................EXIT
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end;
end;
%
dml18j_16 = fix(dml18j_16 + 1);
[yhp,u(sub2ind(size(u),1,1,kod):end),dumvar3,v(sub2ind(size(v),1,kod):end),dumvar5,dml18j_8,dml18j_9]=dstor1(yhp,u(sub2ind(size(u),1,1,kod):end),yhp(sub2ind(size(yhp),1,nfcp1):end),v(sub2ind(size(v),1,kod):end),1,dml18j_8,dml18j_9); yhp(sub2ind(size(yhp),1,nfcp1):end)=dumvar3;
%
% ************************************************
% STORE ORTHONORMALIZATION INFORMATION,
% INITIALIZE INTEGRATION FLAG, AND CONTINUE
% INTEGRATION TO THE NEXT ORTHONORMALIZATION
% POINT OR OUTPUT POINT.
%
z(dml18j_16) = dml15t_4;
if( dml8sz_4==1 && dml18j_13==0 )
dml8sz_1 = s(nfcp1).*dml8sz_1;
end;
if( dml18j_8~=0 )
if( dml8sz_4==1 )
for j=(1):(nfcc), disp({w(j,1)}); end;
end;
for j=(1):(nfcc), for j=(1):(ntp), disp({ip(j,1) ,p(j,1)}); end; end;
end;
dml15t_9(1) = 0;
jon = fix(jon + 1);
% ......EXIT
if( dml18j_6==1 && dml15t_4~=dml15t_8 )
gt(1)=1;
continue;
end;
%
% ************************************************
% CONTINUE INTEGRATION IF WE ARE NOT AT AN
% OUTPUT POINT.
%
% ............EXIT
if( idid==1 )
gt(2)=1;
continue;
end;
%
% *********************************************
% CONTINUE INTEGRATION IF WE ARE NOT AT AN
% OUTPUT POINT.
%
% ...............EXIT
elseif( idid==1 ) ;
gt(2)=1;
continue;
% ......EXIT
end;
else;
dml15t_9(1) = 1;
% ......EXIT
if( idid==-1 )
gt(2)=1;
continue;
end;
iflag = fix(20 - idid);
% .....................EXIT
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
end;
%
% STORAGE OF HOMOGENEOUS SOLUTIONS IN U AND THE PARTICULAR
% SOLUTION IN V AT THE OUTPUT POINTS.
%
break;
end;
[u(sub2ind(size(u),1,1,kod):end),yhp,v(sub2ind(size(v),1,kod):end),dumvar4,dumvar5,dml18j_8,dml18j_9]=dstor1(u(sub2ind(size(u),1,1,kod):end),yhp,v(sub2ind(size(v),1,kod):end),yhp(sub2ind(size(yhp),1,nfcp1):end),0,dml18j_8,dml18j_9); yhp(sub2ind(size(yhp),1,nfcp1):end)=dumvar4;
end;
% ***************************************************************
% ***************************************************************
%
iflag = 0;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine drkfab
%DECK DRKFS
|
|
