| [ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc]=rkfab(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]=rkfab(ncomp,xpts,nxpts,nfc,iflag,z,mxnon,p,ntp,ip,yhp,niv,u,v,w,s,stowa,g,work,iwork,nfcc);
persistent gt idid ipar j jflag jon kod kopp nfcp1 non xxop ;
global ml18jr_1; if isempty(ml18jr_1), ml18jr_1=0; end;
global ml8sz_1; if isempty(ml8sz_1), ml8sz_1=0; end;
global ml15to_2; if isempty(ml15to_2), ml15to_2=0; end;
global ml15to_1; if isempty(ml15to_1), ml15to_1=0; end;
global ml18jr_2; if isempty(ml18jr_2), ml18jr_2=0; end;
global ml15to_3; if isempty(ml15to_3), ml15to_3=0; end;
global ml18jr_3; if isempty(ml18jr_3), ml18jr_3=0; end;
global ml15to_4; if isempty(ml15to_4), ml15to_4=0; end;
global ml15to_5; if isempty(ml15to_5), ml15to_5=0; end;
global ml15to_6; if isempty(ml15to_6), ml15to_6=0; end;
global ml15to_8; if isempty(ml15to_8), ml15to_8=0; end;
global ml15to_7; if isempty(ml15to_7), ml15to_7=0; end;
global ml8sz_2; if isempty(ml8sz_2), ml8sz_2=0; end;
if isempty(xxop), xxop=0; end;
global ml18jr_18; if isempty(ml18jr_18), ml18jr_18=0; end;
if isempty(idid), idid=0; end;
global ml8sz_3; if isempty(ml8sz_3), ml8sz_3=0; end;
global ml18jr_11; if isempty(ml18jr_11), ml18jr_11=0; end;
global ml15to_9; if isempty(ml15to_9), ml15to_9=zeros(1,15); end;
global ml8sz_4; if isempty(ml8sz_4), ml8sz_4=0; end;
global ml18jr_12; if isempty(ml18jr_12), ml18jr_12=0; end;
if isempty(ipar), ipar=0; end;
global ml15to_10; if isempty(ml15to_10), ml15to_10=0; end;
global ml8sz_5; if isempty(ml8sz_5), ml8sz_5=0; end;
if isempty(j), j=0; end;
if isempty(jflag), jflag=0; end;
if isempty(jon), jon=0; end;
global ml17bw_4; if isempty(ml17bw_4), ml17bw_4=0; end;
global ml17bw_13; if isempty(ml17bw_13), ml17bw_13=0; end;
global ml17bw_14; if isempty(ml17bw_14), ml17bw_14=0; end;
global ml17bw_5; if isempty(ml17bw_5), ml17bw_5=0; end;
global ml17bw_6; if isempty(ml17bw_6), ml17bw_6=0; end;
global ml17bw_7; if isempty(ml17bw_7), ml17bw_7=0; end;
global ml17bw_8; if isempty(ml17bw_8), ml17bw_8=0; end;
global ml17bw_9; if isempty(ml17bw_9), ml17bw_9=0; end;
global ml17bw_10; if isempty(ml17bw_10), ml17bw_10=0; end;
global ml17bw_11; if isempty(ml17bw_11), ml17bw_11=0; end;
global ml17bw_12; if isempty(ml17bw_12), ml17bw_12=0; end;
global ml17bw_17; if isempty(ml17bw_17), ml17bw_17=0; end;
global ml17bw_1; if isempty(ml17bw_1), ml17bw_1=0; end;
global ml15to_11; if isempty(ml15to_11), ml15to_11=0; end;
if isempty(kod), kod=0; end;
global ml15to_12; if isempty(ml15to_12), ml15to_12=0; end;
if isempty(kopp), kopp=0; end;
global ml17bw_15; if isempty(ml17bw_15), ml17bw_15=0; end;
global ml17bw_16; if isempty(ml17bw_16), ml17bw_16=0; end;
global ml17bw_18; if isempty(ml17bw_18), ml17bw_18=0; end;
global ml15to_13; if isempty(ml15to_13), ml15to_13=0; end;
global ml15to_14; if isempty(ml15to_14), ml15to_14=0; end;
global ml18jr_7; if isempty(ml18jr_7), ml18jr_7=0; end;
global ml8sz_6; if isempty(ml8sz_6), ml8sz_6=0; end;
global ml18jr_8; if isempty(ml18jr_8), ml18jr_8=0; end;
global ml17bw_3; if isempty(ml17bw_3), ml17bw_3=0; end;
global ml17bw_2; if isempty(ml17bw_2), ml17bw_2=0; end;
global ml18jr_10; if isempty(ml18jr_10), ml18jr_10=0; end;
global ml18jr_15; if isempty(ml18jr_15), ml18jr_15=0; end;
global ml18jr_17; if isempty(ml18jr_17), ml18jr_17=0; end;
global ml8sz_7; if isempty(ml8sz_7), ml8sz_7=0; end;
if isempty(nfcp1), nfcp1=0; end;
global ml18jr_5; if isempty(ml18jr_5), ml18jr_5=0; end;
if isempty(non), non=0; end;
global ml18jr_6; if isempty(ml18jr_6), ml18jr_6=0; end;
global ml18jr_13; if isempty(ml18jr_13), ml18jr_13=0; end;
global ml15to_15; if isempty(ml15to_15), ml15to_15=0; end;
global ml18jr_9; if isempty(ml18jr_9), ml18jr_9=0; end;
global ml18jr_14; if isempty(ml18jr_14), ml18jr_14=0; end;
global ml18jr_16; if isempty(ml18jr_16), ml18jr_16=0; end;
global ml18jr_4; if isempty(ml18jr_4), ml18jr_4=0; end;
if isempty(gt), gt=zeros(1,2); end;
%***BEGIN PROLOGUE RKFAB
%***SUBSIDIARY
%***PURPOSE Subsidiary to BVSUP
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (RKFAB-S, DRKFAB-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
%
% subroutine RKFAB 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 BVSUP
%***ROUTINES CALLED BVDER, DEABM, DERKF, REORT, STOR1
%***COMMON BLOCKS ML15TO, ML17BW, ML18JR, ML8SZ
%***REVISION HISTORY (YYMMDD)
% 750601 DATE WRITTEN
% 890831 Modified array declarations. (WRB)
% 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 RKFAB
%
p_shape=size(p);p=reshape([p(:).',zeros(1,ceil(numel(p)./prod([ntp])).*prod([ntp])-numel(p))],ntp,[]);
ip_shape=size(ip);ip=reshape([ip(:).',zeros(1,ceil(numel(ip)./prod([nfcc])).*prod([nfcc])-numel(ip))],nfcc,[]);
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,[]);
z_shape=size(z);z=reshape(z,1,[]);
yhp_shape=size(yhp);yhp=reshape([yhp(:).',zeros(1,ceil(numel(yhp)./prod([ncomp])).*prod([ncomp])-numel(yhp))],ncomp,[]);
xpts_shape=size(xpts);xpts=reshape(xpts,1,[]);
s_shape=size(s);s=reshape(s,1,[]);
stowa_shape=size(stowa);stowa=reshape(stowa,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
g_shape=size(g);g=reshape(g,1,[]);
%
% **********************************************************************
%
% common :: ;
%% common /ml8sz / c , xsav , igofx , inhomo , ivp , ncompd , nfcd;
%% common /ml8sz / ml8sz_1 , ml8sz_2 , ml8sz_3 , ml8sz_4 , ml8sz_5 , ml8sz_6 , ml8sz_7;
% common :: ;
%% common /ml15to/ px , pwcnd , tnd , x , xbeg , xend , xot , xop ,info(15) , istkop , knswot , kop , lotjp ,mnswot , nswot;
%% common /ml15to/ ml15to_1 , ml15to_2 , ml15to_3 , ml15to_4 , ml15to_5 , ml15to_6 , ml15to_7 , ml15to_8 ,ml15to_9(15) , ml15to_10 , ml15to_11 , ml15to_12 , ml15to_13 ,ml15to_14 , ml15to_15;
% common :: ;
%% common /ml18jr/ ae , re , tol , nxptsd , nic , nopg , mxnond ,ndisk , ntape , neq , indpvt , integ , nps ,ntpd , neqivp , numort , nfccd , icoco;
%% common /ml18jr/ ml18jr_1 , ml18jr_2 , ml18jr_3 , ml18jr_4 , ml18jr_5 , ml18jr_6 , ml18jr_7 ,ml18jr_8 , ml18jr_9 , ml18jr_10 , ml18jr_11 , ml18jr_12 , ml18jr_13 ,ml18jr_14 , ml18jr_15 , ml18jr_16 , ml18jr_17 , ml18jr_18;
% common :: ;
%% common /ml17bw/ kkkzpw , needw , neediw , k1 , k2 , k3 , k4 , k5 ,k6 , k7 , k8 , k9 , k10 , k11 , l1 , l2 , kkkint ,lllint;
%% common /ml17bw/ ml17bw_1 , ml17bw_2 , ml17bw_3 , ml17bw_4 , ml17bw_5 , ml17bw_6 , ml17bw_7 , ml17bw_8 ,ml17bw_9 , ml17bw_10 , ml17bw_11 , ml17bw_12 , ml17bw_13 , ml17bw_14 , ml17bw_15 , ml17bw_16 , ml17bw_17 ,ml17bw_18;
%
%
% **********************************************************************
% INITIALIZATION OF COUNTERS AND VARIABLES.
%
%***FIRST EXECUTABLE STATEMENT RKFAB
kod = 1;
non = 1;
ml15to_4 = ml15to_5;
jon = 1;
ml15to_9(1) = 0;
ml15to_9(2) = 0;
ml15to_9(3) = 1;
ml15to_9(4) = 1;
work(1) = ml15to_6;
if( ml18jr_6~=0 )
ml15to_9(3) = 0;
if( ml15to_4==z(1) )
jon = 2;
end;
end;
nfcp1 = fix(nfc + 1);
%
% **********************************************************************
% *****BEGINNING OF INTEGRATION LOOP AT OUTPUT POINTS.******************
% **********************************************************************
%
for kopp = 2 : nxpts;
ml15to_12 = fix(kopp);
%
while( true );
ml15to_8 = xpts(ml15to_12);
if( ml18jr_8==0 )
kod = fix(ml15to_12);
end;
%
% STEP BY STEP INTEGRATION LOOP BETWEEN OUTPUT POINTS.
%
while( true );
xxop = ml15to_8;
if( ml18jr_6~=0 )
if( ml15to_6>ml15to_5 && ml15to_8>z(jon) )
xxop = z(jon);
end;
if( ml15to_6<ml15to_5 && ml15to_8<z(jon) )
xxop = z(jon);
end;
end;
%
% **********************************************************************
while( true );
gt(:)=0;
if( ml18jr_12==2 )
% DEABM INTEGRATOR
%
[dumvar1,ml18jr_10,ml15to_4,yhp,xxop,ml15to_9,ml18jr_2,ml18jr_1,idid,work,ml17bw_17,iwork,ml17bw_18,g,ipar]=deabm(@bvder,ml18jr_10,ml15to_4,yhp,xxop,ml15to_9,ml18jr_2,ml18jr_1,idid,work,ml17bw_17,iwork,ml17bw_18,g,ipar);
else;
% DERKF INTEGRATOR
%
[dumvar1,ml18jr_10,ml15to_4,yhp,xxop,ml15to_9,ml18jr_2,ml18jr_1,idid,work,ml17bw_17,iwork,ml17bw_18,g,ipar]=derkf(@bvder,ml18jr_10,ml15to_4,yhp,xxop,ml15to_9,ml18jr_2,ml18jr_1,idid,work,ml17bw_17,iwork,ml17bw_18,g,ipar);
end;
if( idid>=1 )
%
% **********************************************************************
% GRAM-SCHMIDT ORTHOGONALIZATION TEST FOR ORTHONORMALIZATION
% (TEMPORARILY USING U AND V IN THE TEST)
%
while (1);
if( ml18jr_6==0 )
jflag = 1;
if( ml8sz_4==3 && ml15to_4==ml15to_6 )
jflag = 3;
end;
elseif( xxop==z(jon) ) ;
jflag = 2;
else;
break;
end;
%
if( ml18jr_8==0 )
non = fix(ml18jr_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]=reort(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;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_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;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_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;
return;
else;
%
if( jflag==10 )
gt(1)=1;
break;
end;
%
if( jflag==0 )
%
% **********************************************************************
% STORE ORTHONORMALIZED VECTORS INTO SOLUTION VECTORS.
%
if( ml18jr_16>=mxnon )
if( ml15to_4~=ml15to_6 )
iflag = 13;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_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;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_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;
return;
end;
end;
%
ml18jr_16 = fix(ml18jr_16 + 1);
[yhp,u(sub2ind(size(u),1,1,kod):end),dumvar3,v(sub2ind(size(v),1,kod):end),dumvar5,ml18jr_8,ml18jr_9]=stor1(yhp,u(sub2ind(size(u),1,1,kod):end),yhp(sub2ind(size(yhp),1,nfcp1):end),v(sub2ind(size(v),1,kod):end),1,ml18jr_8,ml18jr_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(ml18jr_16) = ml15to_4;
if( ml8sz_4==1 && ml18jr_13==0 )
ml8sz_1 = s(nfcp1).*ml8sz_1;
end;
if( ml18jr_8~=0 )
if( ml8sz_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;
ml15to_9(1) = 0;
jon = fix(jon + 1);
if( ml18jr_6==1 && ml15to_4~=ml15to_8 )
break;
end;
end;
end;
%
% **********************************************************************
% CONTINUE INTEGRATION IF WE ARE NOT AT AN OUTPUT POINT.
%
break;
end;
if(gt(1)==1)
break;
end;
if( idid~=1 )
gt(2)=1;
break;
end;
else;
ml15to_9(1) = 1;
if( idid~=-1 )
iflag = fix(20 - idid);
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_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;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_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;
return;
end;
end;
end;
if(any(gt==1))
break;
end;
end;
if(gt(2)==1)
break;
end;
end;
%
% STORAGE OF HOMOGENEOUS SOLUTIONS IN U AND THE PARTICULAR
% SOLUTION IN V AT THE OUTPUT POINTS.
%
[u(sub2ind(size(u),1,1,kod):end),yhp,v(sub2ind(size(v),1,kod):end),dumvar4,dumvar5,ml18jr_8,ml18jr_9]=stor1(u(sub2ind(size(u),1,1,kod):end),yhp,v(sub2ind(size(v),1,kod):end),yhp(sub2ind(size(yhp),1,nfcp1):end),0,ml18jr_8,ml18jr_9); yhp(sub2ind(size(yhp),1,nfcp1):end)=dumvar4;
end;
% **********************************************************************
% **********************************************************************
%
iflag = 0;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_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;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_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;
end %subroutine rkfab
%DECK RPQR79
|
|
