| [n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5]=bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5); |
function [n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5]=bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5);
%***BEGIN PROLOGUE BISECT
%***PURPOSE Compute the eigenvalues of a symmetric tridiagonal matrix
% in a given interval using Sturm sequencing.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4A5, D4C2A
%***TYPE SINGLE PRECISION (BISECT-S)
%***KEYWORDS EIGENVALUES, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the bisection technique
% in the ALGOL procedure TRISTURM by Peters and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
%
% This subroutine finds those eigenvalues of a TRIDIAGONAL
% SYMMETRIC matrix which lie in a specified interval,
% using bisection.
%
% On INPUT
%
% N is the order of the matrix. N is an INTEGER variable.
%
% EPS1 is an absolute error tolerance for the computed
% eigenvalues. If the input EPS1 is non-positive,
% it is reset for each submatrix to a default value,
% namely, minus the product of the relative machine
% precision and the 1-norm of the submatrix.
% EPS1 is a REAL variable.
%
% D contains the diagonal elements of the input matrix.
% D is a one-dimensional REAL array, dimensioned D(N).
%
% E contains the subdiagonal elements of the input matrix
% in its last N-1 positions. E(1) is arbitrary.
% E is a one-dimensional REAL array, dimensioned E(N).
%
% E2 contains the squares of the corresponding elements of E.
% E2(1) is arbitrary. E2 is a one-dimensional REAL array,
% dimensioned E2(N).
%
% LB and UB define the interval to be searched for eigenvalues.
% If LB is not less than UB, no eigenvalues will be found.
% LB and UB are REAL variables.
%
% MM should be set to an upper bound for the number of
% eigenvalues in the interval. WARNING - If more than
% MM eigenvalues are determined to lie in the interval,
% an error return is made with no eigenvalues found.
% MM is an INTEGER variable.
%
% On OUTPUT
%
% EPS1 is unaltered unless it has been reset to its
% (last) default value.
%
% D and E are unaltered.
%
% Elements of E2, corresponding to elements of E regarded
% as negligible, have been replaced by zero causing the
% matrix to split into a direct sum of submatrices.
% E2(1) is also set to zero.
%
% M is the number of eigenvalues determined to lie in (LB,UB).
% M is an INTEGER variable.
%
% W contains the M eigenvalues in ascending order.
% W is a one-dimensional REAL array, dimensioned W(MM).
%
% IND contains in its first M positions the submatrix indices
% associated with the corresponding eigenvalues in W --
% 1 for eigenvalues belonging to the first submatrix from
% the top, 2 for those belonging to the second submatrix, etc.
% IND is an one-dimensional INTEGER array, dimensioned IND(MM).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% 3*N+1 if M exceeds MM. In this case, M contains the
% number of eigenvalues determined to lie in
% (LB,UB).
%
% RV4 and RV5 are one-dimensional REAL arrays used for temporary
% storage, dimensioned RV4(N) and RV5(N).
%
% The ALGOL procedure STURMCNT contained in TRISTURM
% appears in BISECT in-line.
%
% Note that subroutine TQL1 or IMTQL1 is generally faster than
% BISECT, if more than N/4 eigenvalues are to be found.
%
% Questions and comments should be directed to B. S. Garbow,
% Applied Mathematics Division, ARGONNE NATIONAL LABORATORY
% ------------------------------------------------------------------
%
%***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
% Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
% system Routines - EISPACK Guide, Springer-Verlag,
% 1976.
%***ROUTINES CALLED R1MACH
%***REVISION HISTORY (YYMMDD)
% 760101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE BISECT
%
persistent first firstCall gt200 gt300 gt400 gt500 i igo ii isturm j k l m1 m2 machep p q r s s1 s2 t1 t2 tag u v x0 x1 xu ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(p), p=0; end;
if isempty(q), q=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(ii), ii=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(tag), tag=0; end;
if isempty(isturm), isturm=0; end;
if isempty(igo), igo=0; end;
if isempty(gt200), gt200=0; end;
if isempty(gt300), gt300=0; end;
if isempty(gt400), gt400=0; end;
if isempty(gt500), gt500=0; end;
d_shape=size(d);d=reshape(d,1,[]);
e_shape=size(e);e=reshape(e,1,[]);
e2_shape=size(e2);e2=reshape(e2,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
rv4_shape=size(rv4);rv4=reshape(rv4,1,[]);
rv5_shape=size(rv5);rv5=reshape(rv5,1,[]);
if isempty(u), u=0; end;
if isempty(v), v=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(xu), xu=0; end;
if isempty(x0), x0=0; end;
if isempty(x1), x1=0; end;
if isempty(machep), machep=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
ind_shape=size(ind);ind=reshape(ind,1,[]);
if isempty(first), first=false; end;
%
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT BISECT
if( first )
[ machep ]=r1mach(4);
end;
first = false;
%
ierr = 0;
tag = 0;
t1 = lb;
t2 = ub;
% .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES ..........
for i = 1 : n;
if( i~=1 )
s1 = abs(d(i)) + abs(d(i-1));
s2 = s1 + abs(e(i));
if( s2>s1 )
continue;
end;
end;
e2(i) = 0.0e0;
end; i = fix(n+1);
% .......... DETERMINE THE NUMBER OF EIGENVALUES
% IN THE INTERVAL ..........
p = 1;
q = fix(n);
x1 = ub;
isturm = 1;
gt200=0;
gt300=0;
gt400=1;
% .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
% INTERVAL BY THE GERSCHGORIN BOUNDS ..........
igo=1;
while(igo==1);
igo=0;
%goes to 500
while (1);
gt500=0;
while(gt400==0);
if(gt200==0)
if( r==m )
lb = t1;
ub = t2;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv4_shape=zeros(rv4_shape);rv4_shape(:)=rv4(1:numel(rv4_shape));rv4=rv4_shape;
rv5_shape=zeros(rv5_shape);rv5_shape(:)=rv5(1:numel(rv5_shape));rv5=rv5_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
tag = fix(tag + 1);
p = fix(q + 1);
xu = d(p);
x0 = d(p);
u = 0.0e0;
%
for q = p : n;
x1 = u;
u = 0.0e0;
v = 0.0e0;
if( q~=n )
u = abs(e(q+1));
v = e2(q+1);
end;
xu = min(d(q)-(x1+u),xu);
x0 = max(d(q)+(x1+u),x0);
if( v==0.0e0 )
break;
end;
end;
%
x1 = max(abs(xu),abs(x0)).*machep;
if( eps1<=0.0e0 )
eps1 = -x1;
end;
if( p~=q )
x1 = x1.*(q-p+1);
lb = max(t1,xu-x1);
ub = min(t2,x0+x1);
x1 = lb;
isturm = 3;
%gt400
break;
else;
% .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL ..........
%old 600
if( t1>d(p) || d(p)>=t2 )
if( q<n )
igo=1;
continue;
end;
lb = t1;
ub = t2;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv4_shape=zeros(rv4_shape);rv4_shape(:)=rv4(1:numel(rv4_shape));rv4=rv4_shape;
rv5_shape=zeros(rv5_shape);rv5_shape(:)=rv5(1:numel(rv5_shape));rv5=rv5_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
m1 = fix(p);
m2 = fix(p);
rv5(p) = d(p);
gt500=1;
%to 500
break;
end;
%gt200
end;
gt200=0;
%old 200
while (1);
xu = lb;
% .......... FOR I=K STEP -1 UNTIL M1 DO -- ..........
for ii = m1 : k;
i = fix(m1 + k - ii);
if( xu<rv4(i) )
xu = rv4(i);
break;
end;
end;
%
if( x0>rv5(k) )
x0 = rv5(k);
end;
% .......... NEXT BISECTION STEP ..........
x1 =(xu+x0).*0.5e0;
s1 = 2.0e0.*(abs(xu)+abs(x0)+abs(eps1));
s2 = s1 + abs(x0-xu);
if( s2==s1 )
% .......... K-TH EIGENVALUE FOUND ..........
rv5(k) = x1;
k = fix(k - 1);
if( k<m1 )
gt500=1;
%to 500
break;
end;
continue;
end;
break;
%old 200
end;
% .......... IN-LINE PROCEDURE FOR STURM SEQUENCE ..........
break;
%gt400
end;
gt400=0;
if(gt500==0)
%old 400
while (1);
s = fix(p - 1);
u = 1.0e0;
%
for i = p : q;
if( u~=0.0e0 )
v = e2(i)./u;
else;
v = abs(e(i))./machep;
if( e2(i)==0.0e0 )
v = 0.0e0;
end;
end;
u = d(i) - x1 - v;
if( u<0.0e0 )
s = fix(s + 1);
end;
end; i = fix(q+1);
%
if( isturm==1 )
m = fix(s);
x1 = lb;
isturm = 2;
continue;
elseif( isturm==2 ) ;
m = fix(m - s);
if( m>mm )
% .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF
% EIGENVALUES IN INTERVAL ..........
ierr = fix(3.*n + 1);
lb = t1;
ub = t2;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv4_shape=zeros(rv4_shape);rv4_shape(:)=rv4(1:numel(rv4_shape));rv4=rv4_shape;
rv5_shape=zeros(rv5_shape);rv5_shape(:)=rv5(1:numel(rv5_shape));rv5=rv5_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
else;
q = 0;
r = 0;
igo=1;
continue;
end;
elseif( isturm==3 ) ;
m1 = fix(s + 1);
x1 = ub;
isturm = 4;
continue;
elseif( isturm==4 ) ;
m2 = fix(s);
%old 600
if( m1>m2 )
if( q<n )
igo=1;
continue;
end;
lb = t1;
ub = t2;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv4_shape=zeros(rv4_shape);rv4_shape(:)=rv4(1:numel(rv4_shape));rv4=rv4_shape;
rv5_shape=zeros(rv5_shape);rv5_shape(:)=rv5(1:numel(rv5_shape));rv5=rv5_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
% .......... FIND ROOTS BY BISECTION ..........
x0 = ub;
isturm = 5;
%
for i = m1 : m2;
rv5(i) = ub;
rv4(i) = lb;
end; i = fix(m2+1);
% .......... LOOP FOR K-TH EIGENVALUE
% FOR K=M2 STEP -1 UNTIL M1 DO --
% (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) ..........
k = fix(m2);
gt200=1;
% gt200
break;
else;
% .......... REFINE INTERVALS ..........
if( s>=k )
x0 = x1;
else;
xu = x1;
if( s>=m1 )
rv4(s+1) = x1;
if( rv5(s)>x1 )
rv5(s) = x1;
end;
else;
rv4(m1) = x1;
end;
end;
%gt300
x1 =(xu+x0).*0.5e0;
s1 = 2.0e0.*(abs(xu)+abs(x0)+abs(eps1));
s2 = s1 + abs(x0-xu);
if( s2==s1 )
% .......... K-TH EIGENVALUE FOUND ..........
rv5(k) = x1;
k = fix(k - 1);
if( k<m1 )
gt500=1;
%to 500
break;
end;
gt200=1;
break;
end;
end;
break;
%old 400
end;
if(gt500==1)
break;
end;
if(gt200==1)
continue;
end;
%gt500
end;
break;
%do after igo
end;
gt500=0;
% .......... ORDER EIGENVALUES TAGGED WITH THEIR
% SUBMATRIX ASSOCIATIONS ..........
s = fix(r);
r = fix(r + m2 - m1 + 1);
j = 1;
k = fix(m1);
%
for l = 1 : r;
if( j<=s )
if( k>m2 )
break;
end;
if( rv5(k)>=w(l) )
j = fix(j + 1);
continue;
else;
%
for ii = j : s;
i = fix(l + s - ii);
w(i+1) = w(i);
ind(i+1) = fix(ind(i));
end; ii = fix(s+1);
end;
end;
%
w(l) = rv5(k);
ind(l) = fix(tag);
k = fix(k + 1);
end;
%
if( q<n )
igo=1;
continue;
end;
lb = t1;
ub = t2;
%igo
end;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv4_shape=zeros(rv4_shape);rv4_shape(:)=rv4(1:numel(rv4_shape));rv4=rv4_shape;
rv5_shape=zeros(rv5_shape);rv5_shape(:)=rv5(1:numel(rv5_shape));rv5=rv5_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
end %subroutine bisect
%DECK BKIAS
|
|
