| [ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag]=dreort(ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag); |
function [ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag]=dreort(ncomp,y,yp,yhp,niv,w,s,p,ip,stowa,iflag);
%***BEGIN PROLOGUE DREORT
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBVSUP
%***LIBRARY SLATEC
%***TYPE doubleprecision (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 DBVSUP
%***ROUTINES CALLED DDOT, DMGSBV, DSTOR1, DSTWAY
%***COMMON BLOCKS DML15T, DML18J, DML8SZ
%***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 DREORT
%
persistent dnd dndt dx gt ijk j k kk l mflag nfcp srp vnorm wcnd ypnm ;
global dml18j_18; if isempty(dml18j_18), dml18j_18=0; end;
global dml8sz_3; if isempty(dml8sz_3), dml8sz_3=0; end;
if isempty(ijk), ijk=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;
ip_shape=size(ip);ip=reshape(ip,1,[]);
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(k), k=0; end;
if isempty(kk), kk=0; end;
global dml15t_11; if isempty(dml15t_11), dml15t_11=0; end;
global dml15t_12; if isempty(dml15t_12), dml15t_12=0; end;
if isempty(l), l=0; end;
global dml15t_13; if isempty(dml15t_13), dml15t_13=0; end;
if isempty(mflag), mflag=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 dml18j_10; if isempty(dml18j_10), dml18j_10=0; end;
global dml18j_15; if isempty(dml18j_15), dml18j_15=0; end;
global dml8sz_7; if isempty(dml8sz_7), dml8sz_7=0; end;
global dml18j_17; if isempty(dml18j_17), dml18j_17=0; end;
if isempty(nfcp), nfcp=0; end;
global dml18j_5; if isempty(dml18j_5), dml18j_5=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;
if isempty(gt), gt=0; end;
global dml18j_1; if isempty(dml18j_1), dml18j_1=0; end;
global dml8sz_1; if isempty(dml8sz_1), dml8sz_1=0; end;
if isempty(dnd), dnd=0; end;
if isempty(dndt), dndt=0; end;
if isempty(dx), dx=0; end;
p_shape=size(p);p=reshape(p,1,[]);
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,[]);
if isempty(srp), srp=0; end;
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;
if isempty(vnorm), vnorm=0; end;
w_shape=size(w);w=reshape(w,1,[]);
if isempty(wcnd), wcnd=0; end;
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;
global dml8sz_2; if isempty(dml8sz_2), dml8sz_2=0; end;
y_shape=size(y);y=reshape([y(:).',zeros(1,ceil(numel(y)./prod([ncomp])).*prod([ncomp])-numel(y))],ncomp,[]);
yhp_shape=size(yhp);yhp=reshape([yhp(:).',zeros(1,ceil(numel(yhp)./prod([ncomp])).*prod([ncomp])-numel(yhp))],ncomp,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
if isempty(ypnm), ypnm=0; end;
%
% ******************************************************************
%
% common :: ;
%% common /dml8sz/ c , xsav , igofx , inhomo , ivp , ncompd , nfc;
%% 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 , nxpts , nic , nopg , mxnon ,ndisk , ntape , neq , indpvt , integ , nps , ntp ,neqivp , numort , nfcc , 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;
%
% **********************************************************************
% BEGIN BLOCK PERMITTING ...EXITS TO 210
% BEGIN BLOCK PERMITTING ...EXITS TO 10
%***FIRST EXECUTABLE STATEMENT DREORT
nfcp = fix(dml8sz_7 + 1);
%
% CHECK TO SEE IF ORTHONORMALIZATION TEST IS TO BE PERFORMED
%
% ...EXIT
if( iflag==1 )
dml15t_11 = fix(dml15t_11 + 1);
% ...EXIT
if( dml15t_11<dml15t_15 )
% ......EXIT
if((dml15t_6-dml15t_4).*(dml15t_4-dml15t_7)<0.0d0 )
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
return;
end;
end;
end;
[y,yhp,yp,dumvar4]=dstor1(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(dml8sz_7);
ncomp_orig=ncomp; [ncomp,dml8sz_7,y,dumvar4,niv,mflag,s,p,ip,dml8sz_4,yp,w,wcnd]=dmgsbv(ncomp,dml8sz_7,y,ncomp,niv,mflag,s,p,ip,dml8sz_4,yp,w,wcnd); ncomp(dumvar4~=ncomp_orig)=dumvar4(dumvar4~=ncomp_orig);
%
% ************************************************************
%
% CHECK FOR LINEAR DEPENDENCE OF THE SOLUTIONS.
%
if( mflag==0 )
% BEGIN BLOCK PERMITTING ...EXITS TO 190
% BEGIN BLOCK PERMITTING ...EXITS TO 110
%
% ******************************************************
%
% ...EXIT
if( iflag==1 )
%
% TEST FOR ORTHONORMALIZATION
%
% ...EXIT
if( wcnd>=50.0d0.*dml18j_3 )
gt=0;
for ijk = 1 : nfcp;
% ......EXIT
if( s(ijk)>1.0d20 )
gt=1;
end;
end; ijk = fix(nfcp+1);
%
% 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.
%
if(gt==0)
dml15t_15 = fix(dml15t_11);
dml15t_11 = 0;
dml15t_13 = 0;
wcnd = log10(wcnd);
if( wcnd>dml15t_3+3.0d0 )
dml15t_15 = fix(2.*dml15t_15);
end;
if( wcnd<dml15t_2 )
dx = dml15t_4 - dml15t_1;
dnd = dml15t_2 - wcnd;
if( dnd>=4 )
dml15t_15 = fix(fix(dml15t_15./2));
end;
dndt = wcnd - dml15t_3;
if( abs(dx.*dndt)<=dnd.*abs(dml15t_6-dml15t_4) )
dml15t_7 = dml15t_4 + dx.*dndt./dnd;
dml15t_15 = fix(min(dml15t_14,dml15t_15));
dml15t_2 = wcnd;
% ......EXIT
dml15t_1 = dml15t_4;
else;
dml15t_7 = dml15t_6;
dml15t_15 = fix(min(dml15t_14,dml15t_15));
dml15t_2 = wcnd;
dml15t_1 = dml15t_4;
end;
else;
dml15t_7 = dml15t_6;
dml15t_15 = fix(min(dml15t_14,dml15t_15));
dml15t_2 = wcnd;
dml15t_1 = dml15t_4;
end;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
return;
end;
end;
end;
%
% *********************************************************
%
% ORTHONORMALIZATION NECESSARY SO WE NORMALIZE THE
% HOMOGENEOUS SOLUTION VECTORS AND CHANGE W ACCORDINGLY.
%
dml15t_15 = 1;
dml15t_11 = 0;
dml15t_13 = 1;
kk = 1;
l = 1;
for k = 1 : dml18j_17;
% BEGIN BLOCK PERMITTING ...EXITS TO 140
srp = sqrt(p(kk));
if( dml8sz_4==1 )
w(k) = srp.*w(k);
end;
vnorm = 1.0d0./srp;
p(kk) = vnorm;
kk = fix(kk + dml18j_17 + 1 - k);
if( dml8sz_7~=dml18j_17 )
% ......EXIT
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(dml18j_17+1);
%
if( dml8sz_4==1 && dml18j_13~=1 )
%
% NORMALIZE THE PARTICULAR SOLUTION
%
yp_orig=yp; [ypnm ,ncomp,yp,dumvar4,dumvar5]=ddot(ncomp,yp,1,yp,1); yp(dumvar5~=yp_orig)=dumvar5(dumvar5~=yp_orig);
if( ypnm==0.0d0 )
ypnm = 1.0d0;
end;
ypnm = sqrt(ypnm);
s(nfcp) = ypnm;
for j = 1 : ncomp;
yp(j) = yp(j)./ypnm;
end; j = fix(ncomp+1);
for j = 1 : dml18j_17;
w(j) = dml8sz_1.*w(j);
end; j = fix(dml18j_17+1);
end;
%
if( iflag==1 )
[y,yp,yhp,dumvar4,stowa]=dstway(y,yp,yhp,0,stowa);
end;
iflag = 0;
% BEGIN BLOCK PERMITTING ...EXITS TO 40
elseif( iflag==2 ) ;
iflag = 30;
elseif( dml15t_15<=1 && dml15t_13~=0 ) ;
iflag = 30;
else;
%
% RETRIEVE DATA FOR A RESTART AT LAST
% ORTHONORMALIZATION POINT
%
[y,yp,yhp,dumvar4,stowa]=dstway(y,yp,yhp,1,stowa);
dml15t_13 = 1;
dml15t_15 = 1;
dml15t_11 = 0;
dml15t_14 = fix(fix(dml15t_14./2));
dml15t_3 = dml15t_3 + 1.0d0;
% .........EXIT
iflag = 10;
end;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_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;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yhp_shape=zeros(yhp_shape);yhp_shape(:)=yhp(1:numel(yhp_shape));yhp=yhp_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
end %subroutine dreort
%DECK DRF
|
|
