| [ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag]=reort(ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag); |
function [ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag]=reort(ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag);
persistent dnd dndt dx gt ijk j k kk l mflag nfcp srp vnorm wcnd ypnm ;
global ml18jr_1; if isempty(ml18jr_1), ml18jr_1=0; end;
global ml8sz_1; if isempty(ml8sz_1), ml8sz_1=0; end;
if isempty(dnd), dnd=0; end;
if isempty(dndt), dndt=0; end;
if isempty(dx), dx=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;
if isempty(srp), srp=0; end;
global ml15to_3; if isempty(ml15to_3), ml15to_3=0; end;
global ml18jr_3; if isempty(ml18jr_3), ml18jr_3=0; end;
if isempty(vnorm), vnorm=0; end;
if isempty(wcnd), wcnd=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(ypnm), ypnm=0; end;
global ml18jr_18; if isempty(ml18jr_18), ml18jr_18=0; end;
global ml8sz_3; if isempty(ml8sz_3), ml8sz_3=0; end;
if isempty(ijk), ijk=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;
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(k), k=0; end;
if isempty(kk), kk=0; end;
global ml15to_11; if isempty(ml15to_11), ml15to_11=0; end;
global ml15to_12; if isempty(ml15to_12), ml15to_12=0; end;
if isempty(l), l=0; end;
global ml15to_13; if isempty(ml15to_13), ml15to_13=0; end;
if isempty(mflag), mflag=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 ml18jr_10; if isempty(ml18jr_10), ml18jr_10=0; end;
global ml18jr_15; if isempty(ml18jr_15), ml18jr_15=0; end;
global ml8sz_7; if isempty(ml8sz_7), ml8sz_7=0; end;
global ml18jr_17; if isempty(ml18jr_17), ml18jr_17=0; end;
if isempty(nfcp), nfcp=0; end;
global ml18jr_5; if isempty(ml18jr_5), ml18jr_5=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=0; end;
%***BEGIN PROLOGUE REORT
%***SUBSIDIARY
%***PURPOSE Subsidiary to BVSUP
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (REORT-S, DREORT-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
% INPUT
% *********
% Y, YP and YHP = homogeneous solution matrix and particular
% solution vector to be orthonormalized.
% IFLAG = 1 -- store YHP into Y and YP, test for
% reorthonormalization, orthonormalize if needed,
% savemlv restart data.
% 2 -- store YHP into Y and YP, reorthonormalization,
% no restarts.
% (preset orthonormalization mode)
% 3 -- store YHP into Y and YP, reorthonormalization
% (when INHOMO=3 and X=XEND).
% **********************************************************************
% OUTPUT
% *********
% Y, YP = orthonormalized solutions.
% NIV = number of independent vectors returned from DMGSBV.
% IFLAG = 0 -- reorthonormalization was performed.
% 10 -- solution process must be restarted at the last
% orthonormalization point.
% 30 -- solutions are linearly dependent, problem must
% be restarted from the beginning.
% W, P, IP = orthonormalization information.
% **********************************************************************
%
%***SEE ALSO BVSUP
%***ROUTINES CALLED MGSBV, SDOT, STOR1, STWAY
%***COMMON BLOCKS ML15TO, ML18JR, ML8SZ
%***REVISION HISTORY (YYMMDD)
% 750601 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 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 REORT
%
y_shape=size(y);y=reshape([y(:).',zeros(1,ceil(numel(y)./prod([ncomp])).*prod([ncomp])-numel(y))],ncomp,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
s_shape=size(s);s=reshape(s,1,[]);
p_shape=size(p);p=reshape(p,1,[]);
ip_shape=size(ip);ip=reshape(ip,1,[]);
stowa_shape=size(stowa);stowa=reshape(stowa,1,[]);
yhp_shape=size(yhp);yhp=reshape([yhp(:).',zeros(1,ceil(numel(yhp)./prod([ncomp])).*prod([ncomp])-numel(yhp))],ncomp,[]);
%
% **********************************************************************
%
% common :: ;
%% common /ml8sz / c , xsav , igofx , inhomo , ivp , ncompd , nfc;
%% 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 , nxpts , nic , nopg , mxnon ,ndisk , ntape , neq , indpvt , integ , nps , ntp ,neqivp , numort , nfcc , 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;
%
% **********************************************************************
%***FIRST EXECUTABLE STATEMENT REORT
nfcp = fix(ml8sz_7 + 1);
%
% CHECK TO SEE IF ORTHONORMALIZATION TEST IS TO BE PERFORMED
%
if( iflag==1 )
ml15to_11 = fix(ml15to_11 + 1);
if( ml15to_11<ml15to_15 )
if((ml15to_6-ml15to_4).*(ml15to_4-ml15to_7)<0. )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
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;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
return;
end;
end;
end;
[y,yhp,yp,dumvar4]=stor1(y,yhp,yp,yhp(sub2ind(size(yhp),1,nfcp):end),1,0,0); yhp(sub2ind(size(yhp),1,nfcp):end)=dumvar4;
%
% ****************************************
%
% ORTHOGONALIZE THE HOMOGENEOUS SOLUTIONS Y
% AND PARTICULAR SOLUTION YP.
%
niv = fix(ml8sz_7);
ncomp_orig=ncomp; [ncomp,ml8sz_7,y,dumvar4,niv,mflag,s,p,ip,ml8sz_4,yp,w,wcnd]=mgsbv(ncomp,ml8sz_7,y,ncomp,niv,mflag,s,p,ip,ml8sz_4,yp,w,wcnd); ncomp(dumvar4~=ncomp_orig)=dumvar4(dumvar4~=ncomp_orig);
%
% ****************************************
%
% CHECK FOR LINEAR DEPENDENCE OF THE SOLUTIONS.
%
if( mflag==0 )
%
% ****************************************
%
while (1);
gt=0;
if( iflag==1 )
%
% TEST FOR ORTHONORMALIZATION
%
if( wcnd>=50..*ml18jr_3 )
for ijk = 1 : nfcp;
if( s(ijk)>1.0e+20 )
gt=1;
break;
end;
end;
if(gt==1)
break;
end;
%
% use LINEAR EXTRAPOLATION ON LOGARITHMIC VALUES OF THE NORM
% DECREMENTS TO DETERMINE NEXT ORTHONORMALIZATION CHECKPOINT.
% OTHER CONTROLS ON THE NUMBER OF STEPS TO THE NEXT CHECKPOINT
% ARE ADDED FOR SAFETY PURPOSES.
%
ml15to_15 = fix(ml15to_11);
ml15to_11 = 0;
ml15to_13 = 0;
wcnd = log10(wcnd);
if( wcnd>ml15to_3+3. )
ml15to_15 = fix(2.*ml15to_15);
end;
if( wcnd>=ml15to_2 )
ml15to_7 = ml15to_6;
else;
dx = ml15to_4 - ml15to_1;
dnd = ml15to_2 - wcnd;
if( dnd>=4 )
ml15to_15 = fix(fix(ml15to_15./2));
end;
dndt = wcnd - ml15to_3;
if( abs(dx.*dndt)>dnd.*abs(ml15to_6-ml15to_4) )
ml15to_7 = ml15to_6;
else;
ml15to_7 = ml15to_4 + dx.*dndt./dnd;
end;
end;
ml15to_15 = fix(min(ml15to_14,ml15to_15));
ml15to_2 = wcnd;
ml15to_1 = ml15to_4;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
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;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
return;
end;
end;
break;
end;
%
% ****************************************
%
% ORTHONORMALIZATION NECESSARY SO WE NORMALIZE THE HOMOGENEOUS
% SOLUTION VECTORS AND CHANGE W ACCORDINGLY.
%
ml15to_15 = 1;
ml15to_11 = 0;
ml15to_13 = 1;
kk = 1;
l = 1;
for k = 1 : ml18jr_17;
srp = sqrt(p(kk));
if( ml8sz_4==1 )
w(k) = srp.*w(k);
end;
vnorm = 1../srp;
p(kk) = vnorm;
kk = fix(kk + ml18jr_17 + 1 - k);
if( ml8sz_7~=ml18jr_17 )
if( l~=fix(k./2) )
continue;
end;
end;
for j = 1 : ncomp;
y(j,l) = y(j,l).*vnorm;
end; j = fix(ncomp+1);
l = fix(l + 1);
end; k = fix(ml18jr_17+1);
%
if( ml8sz_4==1 && ml18jr_13~=1 )
%
% NORMALIZE THE PARTICULAR SOLUTION
%
yp_orig=yp; [ypnm ,ncomp,yp,dumvar4,dumvar5]=sdot(ncomp,yp,1,yp,1); yp(dumvar5~=yp_orig)=dumvar5(dumvar5~=yp_orig);
if( ypnm==0.0 )
ypnm = 1.0;
end;
ypnm = sqrt(ypnm);
s(nfcp) = ypnm;
for j = 1 : ncomp;
yp(j) = yp(j)./ypnm;
end; j = fix(ncomp+1);
for j = 1 : ml18jr_17;
w(j) = ml8sz_1.*w(j);
end; j = fix(ml18jr_17+1);
end;
%
if( iflag==1 )
[y,yp,yhp,dumvar4,stowa]=stway(y,yp,yhp,0,stowa);
end;
iflag = 0;
else;
if( iflag~=2 )
if( ml15to_15>1 || ml15to_13==0 )
%
% RETRIEVE DATA FOR A RESTART AT LAST ORTHONORMALIZATION POINT
%
[y,yp,yhp,dumvar4,stowa]=stway(y,yp,yhp,1,stowa);
ml15to_13 = 1;
ml15to_15 = 1;
ml15to_11 = 0;
ml15to_14 = fix(fix(ml15to_14./2));
ml15to_3 = ml15to_3 + 1.;
iflag = 10;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
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;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
return;
end;
end;
iflag = 30;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
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;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
return;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
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;
stowa_shape=zeros(stowa_shape);stowa_shape(:)=stowa(1:numel(stowa_shape));stowa=stowa_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
end %subroutine reort
% Structured by FOR_STRUCT, v2.2, on 08/23/2006 at 08:57:10
% Options SET: alpha=l cl=cmt=t cnst=9 dc=k dup=rs i=60,1
% labf=100,100 labx=100,100 nosp=a1c1e0l1r1 o=.fs off=i
% s=9 sp=dl tof90=r v
%DECK RF
|
|
