| [n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5]=tridib(n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5); |
function [n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5]=tridib(n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5);
%***BEGIN PROLOGUE TRIDIB
%***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 (TRIDIB-S)
%***KEYWORDS EIGENVALUES OF A REAL SYMMETRIC MATRIX, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure BISECT,
% NUM. MATH. 9, 386-393(1967) by Barth, Martin, and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971).
%
% This subroutine finds those eigenvalues of a TRIDIAGONAL
% SYMMETRIC matrix between specified boundary indices,
% 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 eigen-
% values. 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 symmetric tridiagonal
% matrix. D is a one-dimensional REAL array, dimensioned D(N).
%
% E contains the subdiagonal elements of the symmetric
% tridiagonal 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).
%
% M11 specifies the lower boundary index for the set of desired
% eigenvalues. M11 is an INTEGER variable.
%
% M specifies the number of eigenvalues desired. The upper
% boundary index M22 is then obtained as M22=M11+M-1.
% M 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.
%
% LB and UB define an interval containing exactly the desired
% eigenvalues. LB and UB are REAL variables.
%
% W contains, in its first M positions, the eigenvalues
% between indices M11 and M22 in ascending order.
% W is a one-dimensional REAL array, dimensioned W(M).
%
% 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(M).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% 3*N+1 if multiple eigenvalues at index M11 make
% unique selection of LB impossible,
% 3*N+2 if multiple eigenvalues at index M22 make
% unique selection of UB impossible.
%
% RV4 and RV5 are one-dimensional REAL arrays used for temporary
% storage of the lower and upper bounds for the eigenvalues in
% the bisection process. RV4 and RV5 are dimensioned RV4(N)
% and RV5(N).
%
% Note that subroutine TQL1, IMTQL1, or TQLRAT is generally faster
% than TRIDIB, 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)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE TRIDIB
%
persistent first firstCall gt200 gt300 gt400 gt500 gt600 gt700 gt800 gt900 i ii isturm j k l m1 m2 m22 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(m22), m22=0; end;
if isempty(tag), tag=0; end;
if isempty(isturm), isturm=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;
if isempty(gt600), gt600=0; end;
if isempty(gt700), gt700=0; end;
if isempty(gt800), gt800=0; end;
if isempty(gt900), gt900=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 TRIDIB
if( first )
[ machep ]=r1mach(4);
end;
first = false;
%
ierr = 0;
tag = 0;
xu = d(1);
x0 = d(1);
u = 0.0e0;
% .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DETERMINE AN
% INTERVAL CONTAINING ALL THE EIGENVALUES ..........
for i = 1 : n;
x1 = u;
u = 0.0e0;
if( i~=n )
u = abs(e(i+1));
end;
xu = min(d(i)-(x1+u),xu);
x0 = max(d(i)+(x1+u),x0);
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);
%
x1 = max(abs(xu),abs(x0)).*machep.*n;
xu = xu - x1;
t1 = xu;
x0 = x0 + x1;
t2 = x0;
% .......... DETERMINE AN INTERVAL CONTAINING EXACTLY
% THE DESIRED EIGENVALUES ..........
p = 1;
q = fix(n);
m1 = fix(m11 - 1);
gt200=0;
gt300=0;
gt400=0;
gt500=0;
gt600=0;
gt700=0;
gt800=0;
gt900=0;
if( m1==0 )
gt200=1;
else;
isturm = 1;
end;
while (1);
if(gt900==0)
if(gt800==0)
if(gt700==0)
if(gt600==0)
if(gt500==0)
if(gt400==0)
if(gt300==0)
if(gt200==0)
v = x1;
x1 = xu +(x0-xu).*0.5e0;
if( x1~=v )
gt700=1;
continue;
end;
% .......... SET ERROR -- INTERVAL CANNOT BE FOUND CONTAINING
% EXACTLY THE DESIRED EIGENVALUES ..........
ierr = fix(3.*n + isturm);
%to 1000
break;
%200
end;
gt200=0;
m22 = fix(m1 + m);
if( m22~=n )
x0 = t2;
isturm = 2;
%to 100
continue;
end;
%300
end;
gt300=0;
q = 0;
r = 0;
% .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
% INTERVAL BY THE GERSCHGORIN BOUNDS ..........
%400
end;
gt400=0;
%to 1000
if( r==m )
break;
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;
gt700=1;
continue;
else;
% .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL ..........
if( t1>d(p) || d(p)>=t2 )
gt900=1;
continue;
end;
m1 = fix(p);
m2 = fix(p);
rv5(p) = d(p);
gt800=1;
continue;
end;
%500
end;
gt500=0;
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 ..........
%600
end;
gt600=0;
x1 =(xu+x0).*0.5e0;
s1 = abs(xu) + abs(x0) + abs(eps1);
s2 = s1 + abs(x0-xu)./2.0e0;
if( s2==s1 )
% .......... K-TH EIGENVALUE FOUND ..........
rv5(k) = x1;
k = fix(k - 1);
if( k<m1 )
gt800=1;
continue;
end;
gt500=1;
continue;
end;
% .......... IN-LINE PROCEDURE FOR STURM SEQUENCE ..........
%700
end;
gt700=0;
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 )
if( s<m1 )
xu = x1;
elseif( s==m1 ) ;
xu = x1;
t1 = x1;
gt200=1;
continue;
else;
x0 = x1;
end;
continue;
elseif( isturm==2 ) ;
if( s<m22 )
xu = x1;
elseif( s==m22 ) ;
t2 = x1;
gt300=1;
continue;
else;
x0 = x1;
end;
continue;
elseif( isturm==3 ) ;
m1 = fix(s + 1);
x1 = ub;
isturm = 4;
gt700=1;
continue;
elseif( isturm==4 ) ;
m2 = fix(s);
if( m1>m2 )
gt900=1;
continue;
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);
gt500=1;
continue;
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;
gt600=1;
continue;
end;
% .......... ORDER EIGENVALUES TAGGED WITH THEIR
% SUBMATRIX ASSOCIATIONS ..........
%800
end;
gt800=0;
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;
%
%900
end;
gt900=0;
if( q<n )
gt400=1;
continue;
end;
break;
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 %subroutine tridib
%DECK TRIDQ
|
|
