| [m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd]=dmgsbv(m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd); |
function [m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd]=dmgsbv(m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd);
%***BEGIN PROLOGUE DMGSBV
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBVSUP
%***LIBRARY SLATEC
%***TYPE doubleprecision (MGSBV-S, DMGSBV-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
% Orthogonalize a set of N doubleprecision vectors and determine their
% rank.
%
% **********************************************************************
% INPUT
% **********************************************************************
% M = dimension of vectors.
% N = no. of vectors.
% A = array whose first N cols contain the vectors.
% IA = first dimension of array A (col length).
% NIV = number of independent vectors needed.
% INHOMO = 1 corresponds to having a non-zero particular solution.
% V = particular solution vector (not included in the pivoting).
% INDPVT = 1 means pivoting will not be used.
%
% **********************************************************************
% OUTPUT
% **********************************************************************
% NIV = no. of linear independent vectors in input set.
% A = matrix whose first NIV cols. contain NIV orthogonal vectors
% which span the vector space determined by the input vectors.
% IFLAG
% = 0 success
% = 1 incorrect input
% = 2 rank of new vectors less than N
% P = decomposition matrix. P is upper triangular and
% (old vectors) = (new vectors) * P.
% The old vectors will be reordered due to pivoting.
% The dimension of P must be .GE. N*(N+1)/2.
% ( N*(2*N+1) when N .NE. NFCC )
% IP = pivoting vector. The dimension of IP must be .GE. N.
% ( 2*N when N .NE. NFCC )
% S = square of norms of incoming vectors.
% V = vector which is orthogonal to the vectors of A.
% W = orthogonalization information for the vector V.
% WCND = worst case (smallest) norm decrement value of the
% vectors being orthogonalized (represents a test
% for linear dependence of the vectors).
% **********************************************************************
%
%***SEE ALSO DBVSUP
%***ROUTINES CALLED DDOT, DPRVEC
%***COMMON BLOCKS DML18J, DML5MC
%***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)
% 890921 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DMGSBV
%
persistent dot gt i ip1 ix iz j jk jp jq jy jz k kd kj kp l lix lr m2 nivn nmnr nn np1 nr nrm1 pjp psave ry sv t vl vnorm y ;
if isempty(i), i=0; end;
global dml18j_18; if isempty(dml18j_18), dml18j_18=0; end;
global dml18j_11; if isempty(dml18j_11), dml18j_11=0; end;
global dml18j_12; if isempty(dml18j_12), dml18j_12=0; end;
ip_shape=size(ip);ip=reshape(ip,1,[]);
if isempty(ip1), ip1=0; end;
if isempty(ix), ix=0; end;
if isempty(iz), iz=0; end;
if isempty(j), j=0; end;
if isempty(jk), jk=0; end;
if isempty(jp), jp=0; end;
if isempty(jq), jq=0; end;
if isempty(jy), jy=0; end;
if isempty(jz), jz=0; end;
if isempty(k), k=0; end;
if isempty(kd), kd=0; end;
if isempty(kj), kj=0; end;
if isempty(kp), kp=0; end;
if isempty(l), l=0; end;
if isempty(lix), lix=0; end;
global dml5mc_7; if isempty(dml5mc_7), dml5mc_7=0; end;
if isempty(lr), lr=0; end;
if isempty(m2), m2=0; end;
global dml18j_7; if isempty(dml18j_7), dml18j_7=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 dml18j_17; if isempty(dml18j_17), dml18j_17=0; end;
global dml18j_5; if isempty(dml18j_5), dml18j_5=0; end;
if isempty(nivn), nivn=0; end;
if isempty(nmnr), nmnr=0; end;
if isempty(nn), nn=0; end;
global dml18j_6; if isempty(dml18j_6), dml18j_6=0; end;
if isempty(np1), np1=0; end;
global dml18j_13; if isempty(dml18j_13), dml18j_13=0; end;
if isempty(nr), nr=0; end;
if isempty(nrm1), nrm1=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;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([ia])).*prod([ia])-numel(a))],ia,[]);
global dml18j_1; if isempty(dml18j_1), dml18j_1=0; end;
if isempty(dot), dot=0; end;
global dml5mc_3; if isempty(dml5mc_3), dml5mc_3=0; end;
global dml5mc_6; if isempty(dml5mc_6), dml5mc_6=0; end;
p_shape=size(p);p=reshape(p,1,[]);
if isempty(pjp), pjp=0; end;
if isempty(psave), psave=0; end;
global dml18j_2; if isempty(dml18j_2), dml18j_2=0; end;
if isempty(ry), ry=0; end;
s_shape=size(s);s=reshape(s,1,[]);
global dml5mc_4; if isempty(dml5mc_4), dml5mc_4=0; end;
global dml5mc_2; if isempty(dml5mc_2), dml5mc_2=0; end;
if isempty(sv), sv=0; end;
if isempty(t), t=0; end;
global dml18j_3; if isempty(dml18j_3), dml18j_3=0; end;
global dml5mc_5; if isempty(dml5mc_5), dml5mc_5=0; end;
global dml5mc_1; if isempty(dml5mc_1), dml5mc_1=0; end;
v_shape=size(v);v=reshape(v,1,[]);
if isempty(vl), vl=0; end;
if isempty(vnorm), vnorm=0; end;
w_shape=size(w);w=reshape(w,1,[]);
if isempty(y), y=0; end;
%
%
% 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;
%
% common :: ;
%% common /dml5mc/ uro , sru , eps , sqovfl , twou , fouru , lpar;
%% common /dml5mc/ dml5mc_1 , dml5mc_2 , dml5mc_3 , dml5mc_4 , dml5mc_5 , dml5mc_6 , dml5mc_7;
%
%***FIRST EXECUTABLE STATEMENT DMGSBV
gt=0;
if( m>0 && n>0 && ia>=m )
% BEGIN BLOCK PERMITTING ...EXITS TO 270
% BEGIN BLOCK PERMITTING ...EXITS TO 260
%
jp = 0;
iflag = 0;
np1 = fix(n + 1);
y = 0.0d0;
m2 = fix(fix(m./2));
%
% CALCULATE SQUARE OF NORMS OF INCOMING VECTORS AND SEARCH
% FOR VECTOR WITH LARGEST MAGNITUDE
%
j = 0;
for i = 1 : n;
[vl ,m,dumvar2,dumvar4,dumvar2]=ddot(m,a(sub2ind(size(a),1,i):end),1,a(sub2ind(size(a),1,i):end),1); a(sub2ind(size(a),1,i):end)=dumvar2;
s(i) = vl;
if( n~=dml18j_17 )
j = fix(2.*i - 1);
p(j) = vl;
ip(j) = fix(j);
end;
j = fix(j + 1);
p(j) = vl;
ip(j) = fix(j);
if( vl>y )
y = vl;
ix = fix(i);
end;
end; i = fix(n+1);
if( dml18j_11==1 )
ix = 1;
y = p(1);
end;
lix = fix(ix);
if( n~=dml18j_17 )
lix = fix(2.*ix - 1);
end;
p(lix) = p(1);
s(np1) = 0.0d0;
if( inhomo==1 )
v_orig=v; [ s(np1) ,m,v,dumvar4,dumvar5]=ddot(m,v,1,v,1); v(dumvar5~=v_orig)=dumvar5(dumvar5~=v_orig);
end;
wcnd = 1.0d0;
nivn = fix(niv);
niv = 0;
%
% ...EXIT
if( y~=0.0d0 )
% *********************************************************
for nr = 1 : n;
% BEGIN BLOCK PERMITTING ...EXITS TO 230
% ......EXIT
if( nivn==niv )
break;
end;
niv = fix(nr);
if( ix~=nr )
%
% PIVOTING OF COLUMNS OF P MATRIX
%
nn = fix(n);
lix = fix(ix);
lr = fix(nr);
if( n~=dml18j_17 )
nn = fix(dml18j_17);
lix = fix(2.*ix - 1);
lr = fix(2.*nr - 1);
end;
if( nr~=1 )
kd = fix(lix - lr);
kj = fix(lr);
nrm1 = fix(lr - 1);
for j = 1 : nrm1;
psave = p(kj);
jk = fix(kj + kd);
p(kj) = p(jk);
p(jk) = psave;
kj = fix(kj + nn - j);
end; j = fix(nrm1+1);
jy = fix(jk + nmnr);
jz = fix(jy - kd);
p(jy) = p(jz);
end;
iz = fix(ip(lix));
ip(lix) = fix(ip(lr));
ip(lr) = fix(iz);
sv = s(ix);
s(ix) = s(nr);
s(nr) = sv;
if( n~=dml18j_17 )
if( nr~=1 )
kj = fix(lr + 1);
for k = 1 : nrm1;
psave = p(kj);
jk = fix(kj + kd);
p(kj) = p(jk);
p(jk) = psave;
kj = fix(kj + dml18j_17 - k);
end; k = fix(nrm1+1);
end;
iz = fix(ip(lix+1));
ip(lix+1) = fix(ip(lr+1));
ip(lr+1) = fix(iz);
end;
%
% PIVOTING OF COLUMNS OF VECTORS
%
for l = 1 : m;
t = a(l,ix);
a(l,ix) = a(l,nr);
a(l,nr) = t;
end; l = fix(m+1);
end;
%
% CALCULATE P(NR,NR) AS NORM SQUARED OF PIVOTAL
% VECTOR
%
jp = fix(jp + 1);
p(jp) = y;
ry = 1.0d0./y;
nmnr = fix(n - nr);
if( n~=dml18j_17 )
nmnr = fix(dml18j_17 -(2.*nr-1));
jp = fix(jp + 1);
p(jp) = 0.0d0;
kp = fix(jp + nmnr);
p(kp) = y;
end;
if( nr~=n && nivn~=niv )
%
% CALCULATE ORTHOGONAL PROJECTION VECTORS AND
% SEARCH FOR LARGEST NORM
%
y = 0.0d0;
ip1 = fix(nr + 1);
ix = fix(ip1);
% ************************************************
for j = ip1 : n;
[dot ,m,dumvar2,dumvar4,dumvar4]=ddot(m,a(sub2ind(size(a),1,nr):end),1,a(sub2ind(size(a),1,j):end),1); a(sub2ind(size(a),1,nr):end)=dumvar2; a(sub2ind(size(a),1,j):end)=dumvar4;
jp = fix(jp + 1);
jq = fix(jp + nmnr);
if( n~=dml18j_17 )
jq = fix(jq + nmnr - 1);
end;
p(jq) = p(jp) - dot.*(dot.*ry);
p(jp) = dot.*ry;
for i = 1 : m;
a(i,j) = a(i,j) - p(jp).*a(i,nr);
end; i = fix(m+1);
if( n~=dml18j_17 )
kp = fix(jp + nmnr);
jp = fix(jp + 1);
pjp = ry.*dprvec(m,a(sub2ind(size(a),1,nr):end),a(sub2ind(size(a),1,j):end));
p(jp) = pjp;
p(kp) = -pjp;
kp = fix(kp + 1);
p(kp) = ry.*dot;
for k = 1 : m2;
l = fix(m2 + k);
a(k,j) = a(k,j) - pjp.*a(l,nr);
a(l,j) = a(l,j) + pjp.*a(k,nr);
end; k = fix(m2+1);
p(jq) = p(jq) - pjp.*(pjp./ry);
end;
%
% TEST FOR CANCELLATION IN RECURRENCE RELATION
%
if( p(jq)<=s(j).*dml5mc_2 )
[ p(jq) ,m,dumvar2,dumvar4,dumvar2]=ddot(m,a(sub2ind(size(a),1,j):end),1,a(sub2ind(size(a),1,j):end),1); a(sub2ind(size(a),1,j):end)=dumvar2;
end;
if( p(jq)>y )
y = p(jq);
ix = fix(j);
end;
end; j = fix(n+1);
if( n~=dml18j_17 )
jp = fix(kp);
end;
% ************************************************
if( dml18j_11==1 )
ix = fix(ip1);
end;
%
% RECOMPUTE NORM SQUARED OF PIVOTAL VECTOR WITH
% SCALAR PRODUCT
%
[y ,m,dumvar2,dumvar4,dumvar2]=ddot(m,a(sub2ind(size(a),1,ix):end),1,a(sub2ind(size(a),1,ix):end),1); a(sub2ind(size(a),1,ix):end)=dumvar2;
% ............EXIT
if( y<=dml5mc_3.*s(ix) )
gt=1;
break;
end;
wcnd = min(wcnd,y./s(ix));
end;
%
% COMPUTE ORTHOGONAL PROJECTION OF PARTICULAR
% SOLUTION
%
% ...EXIT
if( inhomo==1 )
lr = fix(nr);
if( n~=dml18j_17 )
lr = fix(2.*nr - 1);
end;
w(lr) = ddot(m,a(sub2ind(size(a),1,nr):end),1,v,1).*ry;
for i = 1 : m;
v(i) = v(i) - w(lr).*a(i,nr);
end; i = fix(m+1);
% ...EXIT
if( n~=dml18j_17 )
lr = fix(2.*nr);
w(lr) = ry.*dprvec(m,v,a(sub2ind(size(a),1,nr):end));
for k = 1 : m2;
l = fix(m2 + k);
v(k) = v(k) + w(lr).*a(l,nr);
v(l) = v(l) - w(lr).*a(k,nr);
end; k = fix(m2+1);
end;
end;
end;
% *********************************************************
%
% TEST FOR LINEAR DEPENDENCE OF PARTICULAR SOLUTION
%
% ......EXIT
if(gt==0)
if( inhomo~=1 )
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
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;
return;
end;
if((n>1) &&(s(np1)<1.0) )
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
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;
return;
end;
v_orig=v; [vnorm ,m,v,dumvar4,dumvar5]=ddot(m,v,1,v,1); v(dumvar5~=v_orig)=dumvar5(dumvar5~=v_orig);
if( s(np1)~=0.0d0 )
wcnd = min(wcnd,vnorm./s(np1));
end;
% ......EXIT
if( vnorm>=dml5mc_3.*s(np1) )
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
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;
return;
end;
end;
end;
iflag = 2;
wcnd = dml5mc_3;
else;
iflag = 1;
end;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
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;
end %subroutine dmgsbv
%DECK DMOUT
|
|
