| [nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,ierr,rv1,rv2,rv3,rv4,rv5,rv6]=tsturm(nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,ierr,rv1,rv2,rv3,rv4,rv5,rv6); |
function [nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,ierr,rv1,rv2,rv3,rv4,rv5,rv6]=tsturm(nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,ierr,rv1,rv2,rv3,rv4,rv5,rv6);
%***BEGIN PROLOGUE TSTURM
%***PURPOSE Find those eigenvalues of a symmetric tridiagonal matrix
% in a given interval and their associated eigenvectors by
% Sturm sequencing.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4A5, D4C2A
%***TYPE SINGLE PRECISION (TSTURM-S)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine finds those eigenvalues of a TRIDIAGONAL
% SYMMETRIC matrix which lie in a specified interval and their
% associated eigenvectors, using bisection and inverse iteration.
%
% On Input
%
% NM must be set to the row dimension of the two-dimensional
% array parameter, Z, as declared in the calling program
% dimension statement. NM is an INTEGER variable.
%
% N is the order of the matrix. N is an INTEGER variable.
% N must be less than or equal to NM.
%
% EPS1 is an absolute error tolerance for the computed eigen-
% values. It should be chosen so that the accuracy of these
% eigenvalues is commensurate with relative perturbations of
% the order of the relative machine precision in the matrix
% elements. 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).
%
% 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. MM is an INTEGER variable.
% WARNING - If more than MM eigenvalues are determined to lie
% in the interval, an error return is made with no values or
% vectors found.
%
% 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 if the matrix
% does not split. If the matrix splits, the eigenvalues are
% in ascending order for each submatrix. If a vector error
% exit is made, W contains those values already found. W is a
% one-dimensional REAL array, dimensioned W(MM).
%
% Z contains the associated set of orthonormal eigenvectors.
% If an error exit is made, Z contains those vectors already
% found. Z is a one-dimensional REAL array, dimensioned
% Z(NM,MM).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% 3*N+1 if M exceeds MM no eigenvalues or eigenvectors
% are computed,
% 4*N+J if the eigenvector corresponding to the J-th
% eigenvalue fails to converge in 5 iterations, then
% the eigenvalues and eigenvectors in W and Z should
% be correct for indices 1, 2, ..., J-1.
%
% RV1, RV2, RV3, RV4, RV5, and RV6 are temporary storage arrays,
% dimensioned RV1(N), RV2(N), RV3(N), RV4(N), RV5(N), and
% RV6(N).
%
% The ALGOL procedure STURMCNT contained in TRISTURM
% appears in TSTURM in-line.
%
% 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 TSTURM
%
persistent eps2 eps3 eps4 first firstCall group gt200 gt300 gt400 gt500 i ii ip isturm its j jj k m1 m2 machep norm p q r s s1 s2 t1 t2 u uk 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(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(ip), ip=0; end;
if isempty(jj), jj=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(its), its=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(group), group=0; end;
if isempty(isturm), isturm=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,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
rv1_shape=size(rv1);rv1=reshape(rv1,1,[]);
rv2_shape=size(rv2);rv2=reshape(rv2,1,[]);
rv3_shape=size(rv3);rv3=reshape(rv3,1,[]);
rv4_shape=size(rv4);rv4=reshape(rv4,1,[]);
rv5_shape=size(rv5);rv5=reshape(rv5,1,[]);
rv6_shape=size(rv6);rv6=reshape(rv6,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(uk), uk=0; end;
if isempty(xu), xu=0; end;
if isempty(x0), x0=0; end;
if isempty(x1), x1=0; end;
if isempty(eps2), eps2=0; end;
if isempty(eps3), eps3=0; end;
if isempty(eps4), eps4=0; end;
if isempty(norm), norm=0; end;
if isempty(machep), machep=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(first), first=false; end;
%
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT TSTURM
if( first )
[ machep ]=r1mach(4);
end;
first = false;
%
ierr = 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;
% .......... IN-LINE PROCEDURE FOR STURM SEQUENCE ..........
gt200=0;
gt300=0;
gt400=0;
gt500=0;
while (1);
if(gt500==0)
if(gt400==0)
if(gt300==0)
if(gt200==0)
s = fix(p - 1);
u = 1.0e0;
%
for i = p : q;
if( u==0.0e0 )
v = abs(e(i))./machep;
if( e2(i)==0.0e0 )
v = 0.0e0;
end;
else;
v = e2(i)./u;
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
rv2_shape=zeros(rv2_shape);rv2_shape(:)=rv2(1:numel(rv2_shape));rv2=rv2_shape;
rv3_shape=zeros(rv3_shape);rv3_shape(:)=rv3(1:numel(rv3_shape));rv3=rv3_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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
return;
else;
q = 0;
r = 0;
gt400=1;
break;
end;
elseif( isturm==3 ) ;
m1 = fix(s + 1);
x1 = ub;
isturm = 4;
continue;
elseif( isturm==4 ) ;
m2 = fix(s);
if( m1>m2 )
gt500=1;
break;
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);
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=1;
break;
end;
%to 200
end;
gt200=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 ..........
%to 300
end;
gt300=0;
x1 =(xu+x0).*0.5e0;
s1 = 2.0e0.*(abs(xu)+abs(x0)+abs(eps1));
s2 = s1 + abs(x0-xu);
if( s2~=s1 )
continue;
end;
% .......... K-TH EIGENVALUE FOUND ..........
rv5(k) = x1;
k = fix(k - 1);
if( k>=m1 )
gt200=1;
break;
end;
% .......... FIND VECTORS BY INVERSE ITERATION ..........
norm = abs(d(p));
ip = fix(p + 1);
%
for i = ip : q;
norm = max(norm,abs(d(i))+abs(e(i)));
end; i = fix(q+1);
% .......... EPS2 IS THE CRITERION FOR GROUPING,
% EPS3 REPLACES ZERO PIVOTS AND EQUAL
% ROOTS ARE MODIFIED BY EPS3,
% EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ..........
eps2 = 1.0e-3.*norm;
uk = sqrt(real(q-p+5));
eps3 = uk.*machep.*norm;
eps4 = uk.*eps3;
uk = eps4./sqrt(uk);
group = 0;
s = fix(p);
%
for k = m1 : m2;
r = fix(r + 1);
its = 1;
w(r) = rv5(k);
x1 = rv5(k);
% .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
if( k~=m1 )
if( x1-x0>=eps2 )
group = -1;
end;
group = fix(group + 1);
if( x1<=x0 )
x1 = x0 + eps3;
end;
end;
% .......... ELIMINATION WITH INTERCHANGES AND
% INITIALIZATION OF VECTOR ..........
v = 0.0e0;
%
for i = p : q;
rv6(i) = uk;
if( i~=p )
if( abs(e(i))<abs(u) )
xu = e(i)./u;
rv4(i) = xu;
rv1(i-1) = u;
rv2(i-1) = v;
rv3(i-1) = 0.0e0;
else;
xu = u./e(i);
rv4(i) = xu;
rv1(i-1) = e(i);
rv2(i-1) = d(i) - x1;
rv3(i-1) = 0.0e0;
if( i~=q )
rv3(i-1) = e(i+1);
end;
u = v - xu.*rv2(i-1);
v = -xu.*rv3(i-1);
continue;
end;
end;
u = d(i) - x1 - xu.*v;
if( i~=q )
v = e(i+1);
end;
end; i = fix(q+1);
%
if( u==0.0e0 )
u = eps3;
end;
rv1(q) = u;
rv2(q) = 0.0e0;
rv3(q) = 0.0e0;
% .......... BACK SUBSTITUTION
% FOR I=Q STEP -1 UNTIL P DO -- ..........
while( true );
for ii = p : q;
i = fix(p + q - ii);
rv6(i) =(rv6(i)-u.*rv2(i)-v.*rv3(i))./rv1(i);
v = u;
u = rv6(i);
end; ii = fix(q+1);
% .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
% MEMBERS OF GROUP ..........
if( group~=0 )
%
for jj = 1 : group;
j = fix(r - group - 1 + jj);
xu = 0.0e0;
%
for i = p : q;
xu = xu + rv6(i).*z(i,j);
end; i = fix(q+1);
%
for i = p : q;
%
rv6(i) = rv6(i) - xu.*z(i,j);
end; i = fix(q+1);
end; jj = fix(group+1);
end;
%
norm = 0.0e0;
%
for i = p : q;
norm = norm + abs(rv6(i));
end; i = fix(q+1);
%
if( norm>=1.0e0 )
break;
end;
if( its==5 )
% .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
ierr = fix(4.*n + r);
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
rv2_shape=zeros(rv2_shape);rv2_shape(:)=rv2(1:numel(rv2_shape));rv2=rv2_shape;
rv3_shape=zeros(rv3_shape);rv3_shape(:)=rv3(1:numel(rv3_shape));rv3=rv3_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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
return;
else;
if( norm==0.0e0 )
rv6(s) = eps4;
s = fix(s + 1);
if( s>q )
s = fix(p);
end;
else;
xu = eps4./norm;
%
for i = p : q;
rv6(i) = rv6(i).*xu;
end; i = fix(q+1);
end;
% .......... ELIMINATION OPERATIONS ON NEXT VECTOR
% ITERATE ..........
for i = ip : q;
u = rv6(i);
% .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
% WAS PERFORMED EARLIER IN THE
% TRIANGULARIZATION PROCESS ..........
if( rv1(i-1)==e(i) )
u = rv6(i-1);
rv6(i-1) = rv6(i);
end;
rv6(i) = u - rv4(i).*rv6(i-1);
end; i = fix(q+1);
%
its = fix(its + 1);
end;
end;
% .......... NORMALIZE SO THAT SUM OF SQUARES IS
% 1 AND EXPAND TO FULL ORDER ..........
u = 0.0e0;
%
for i = p : q;
u = u + rv6(i).^2;
end; i = fix(q+1);
%
xu = 1.0e0./sqrt(u);
%
for i = 1 : n;
z(i,r) = 0.0e0;
end; i = fix(n+1);
%
for i = p : q;
z(i,r) = rv6(i).*xu;
end; i = fix(q+1);
%
x0 = x1;
% .......... FORWARD SUBSTITUTION ..........
end;
gt500=1;
break;
% .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
% INTERVAL BY THE GERSCHGORIN BOUNDS ..........
%to 400
end;
gt400=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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
rv2_shape=zeros(rv2_shape);rv2_shape(:)=rv2(1:numel(rv2_shape));rv2=rv2_shape;
rv3_shape=zeros(rv3_shape);rv3_shape(:)=rv3(1:numel(rv3_shape));rv3=rv3_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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
return;
end;
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;
continue;
% .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL ..........
elseif( t1<=d(p) && d(p)<t2 ) ;
r = fix(r + 1);
%
for i = 1 : n;
z(i,r) = 0.0e0;
end; i = fix(n+1);
%
w(r) = d(p);
z(p,r) = 1.0e0;
end;
%
%to 500
end;
gt500=0;
if( q<n )
gt400=1;
break;
end;
break;
%100
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
rv2_shape=zeros(rv2_shape);rv2_shape(:)=rv2(1:numel(rv2_shape));rv2=rv2_shape;
rv3_shape=zeros(rv3_shape);rv3_shape(:)=rv3(1:numel(rv3_shape));rv3=rv3_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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
return;
end %subroutine tsturm
%DECK U11LS
|
|
