| [n,eps1,d,e,e2,m,w,ind,bd,type,idef,ierr]=ratqr(n,eps1,d,e,e2,m,w,ind,bd,type,idef,ierr); |
function [n,eps1,d,e,e2,m,w,ind,bd,type,idef,ierr]=ratqr(n,eps1,d,e,e2,m,w,ind,bd,type,idef,ierr);
%***BEGIN PROLOGUE RATQR
%***PURPOSE Compute the largest or smallest eigenvalues of a symmetric
% tridiagonal matrix using the rational QR method with Newton
% correction.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4A5, D4C2A
%***TYPE SINGLE PRECISION (RATQR-S)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure RATQR,
% NUM. MATH. 11, 264-272(1968) by REINSCH and BAUER.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 257-265(1971).
%
% This subroutine finds the algebraically smallest or largest
% eigenvalues of a SYMMETRIC TRIDIAGONAL matrix by the
% rational QR method with Newton corrections.
%
% On Input
%
% N is the order of the matrix. N is an INTEGER variable.
%
% EPS1 is a theoretical absolute error tolerance for the
% computed eigenvalues. If the input EPS1 is non-positive, or
% indeed smaller than its default value, it is reset at each
% iteration to the respective default value, namely, the
% product of the relative machine precision and the magnitude
% of the current eigenvalue iterate. The theoretical absolute
% error in the K-th eigenvalue is usually not greater than
% K times EPS1. 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 in
% its last N-1 positions. E2(1) is arbitrary. E2 is a one-
% dimensional REAL array, dimensioned E2(N).
%
% M is the number of eigenvalues to be found. M is an INTEGER
% variable.
%
% IDEF should be set to 1 if the input matrix is known to be
% positive definite, to -1 if the input matrix is known to
% be negative definite, and to 0 otherwise. IDEF is an
% INTEGER variable.
%
% TYPE should be set to true if the smallest eigenvalues are
% to be found, and to false if the largest eigenvalues are
% to be found. TYPE is a LOGICAL variable.
%
% On Output
%
% EPS1 is unaltered unless it has been reset to its
% (last) default value.
%
% D and E are unaltered (unless W overwrites D).
%
% 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 set
% to 0.0e0 if the smallest eigenvalues have been found, and
% to 2.0e0 if the largest eigenvalues have been found. E2
% is otherwise unaltered (unless overwritten by BD).
%
% W contains the M algebraically smallest eigenvalues in
% ascending order, or the M largest eigenvalues in descending
% order. If an error exit is made because of an incorrect
% specification of IDEF, no eigenvalues are found. If the
% Newton iterates for a particular eigenvalue are not monotone,
% the best estimate obtained is returned and IERR is set.
% W is a one-dimensional REAL array, dimensioned W(N). W need
% not be distinct from D.
%
% 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(N).
%
% BD contains refined bounds for the theoretical errors of the
% corresponding eigenvalues in W. These bounds are usually
% within the tolerance specified by EPS1. BD is a one-
% dimensional REAL array, dimensioned BD(N). BD need not be
% distinct from E2.
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% 6*N+1 if IDEF is set to 1 and TYPE to true
% when the matrix is NOT positive definite, or
% if IDEF is set to -1 and TYPE to false
% when the matrix is NOT negative definite,
% no eigenvalues are computed, or
% M is greater than N,
% 5*N+K if successive iterates to the K-th eigenvalue
% are NOT monotone increasing, where K refers
% to the last such occurrence.
%
% Note that subroutine TRIDIB is generally faster and more
% accurate than RATQR if the eigenvalues are clustered.
%
% 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 RATQR
%
persistent delta ep err f first firstCall gt i ii j jdef jj k k1 machep p q qp r s tot ; 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(ii), ii=0; end;
if isempty(jj), jj=0; end;
if isempty(k1), k1=0; end;
if isempty(jdef), jdef=0; end;
if isempty(gt), gt=zeros(1,4); 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,[]);
bd_shape=size(bd);bd=reshape(bd,1,[]);
if isempty(f), f=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(ep), ep=0; end;
if isempty(qp), qp=0; end;
if isempty(err), err=0; end;
if isempty(tot), tot=0; end;
if isempty(delta), delta=0; end;
if isempty(machep), machep=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 RATQR
if( first )
[ machep ]=r1mach(4);
end;
first = false;
%
ierr = 0;
jdef = fix(idef);
% .......... COPY D ARRAY INTO W ..........
for i = 1 : n;
w(i) = d(i);
end; i = fix(n+1);
%
gt(:)=0;
while (1);
if(gt(4)==0)
if(gt(1)==0)
if( ~(type) )
j = 1;
gt(4)=1;
continue;
end;
end;
gt(1)=0;
err = 0.0e0;
s = 0.0e0;
% .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DEFINE
% INITIAL SHIFT FROM LOWER GERSCHGORIN BOUND.
% COPY E2 ARRAY INTO BD ..........
tot = w(1);
q = 0.0e0;
j = 0;
%
for i = 1 : n;
p = q;
if( i==1 )
e2(i) = 0.0e0;
elseif( p<=machep.*(abs(d(i))+abs(d(i-1))) ) ;
e2(i) = 0.0e0;
end;
bd(i) = e2(i);
% .......... COUNT ALSO IF ELEMENT OF E2 HAS UNDERFLOWED ..........
if( e2(i)==0.0e0 )
j = fix(j + 1);
end;
ind(i) = fix(j);
q = 0.0e0;
if( i~=n )
q = abs(e(i+1));
end;
tot = min(w(i)-p-q,tot);
end; i = fix(n+1);
%
if( jdef==1 && tot<0.0e0 )
tot = 0.0e0;
else;
%
for i = 1 : n;
w(i) = w(i) - tot;
%
end; i = fix(n+1);
end;
%
for k = 1 : m;
gt(3)=0;
% .......... NEXT QR TRANSFORMATION ..........
while (1);
tot = tot + s;
delta = w(n) - s;
i = fix(n);
f = abs(machep.*tot);
if( eps1<f )
eps1 = f;
end;
if( delta>eps1 )
% .......... REPLACE SMALL SUB-DIAGONAL SQUARES BY ZERO
% TO REDUCE THE INCIDENCE OF UNDERFLOWS ..........
if( k~=n )
k1 = fix(k + 1);
for j = k1 : n;
if( bd(j)<=(machep.*(w(j)+w(j-1))).^2 )
bd(j) = 0.0e0;
end;
end; j = fix(n+1);
end;
%
f = bd(n)./delta;
qp = delta + f;
p = 1.0e0;
if( k~=n )
k1 = fix(n - k);
% .......... FOR I=N-1 STEP -1 UNTIL K DO -- ..........
for ii = 1 : k1;
i = fix(n - ii);
q = w(i) - s - f;
r = q./qp;
p = p.*r + 1.0e0;
ep = f.*r;
w(i+1) = qp + ep;
delta = q - ep;
if( delta>eps1 )
f = bd(i)./q;
qp = delta + f;
bd(i+1) = qp.*ep;
else;
if( delta>=(-eps1) )
gt(3)=1;
break;
end;
% .......... SET ERROR -- IDEF SPECIFIED INCORRECTLY ..........
ierr = fix(6.*n + 1);
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;
bd_shape=zeros(bd_shape);bd_shape(:)=bd(1:numel(bd_shape));bd=bd_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
end;
if(gt(3)==1)
break;
end;
end;
%
w(k) = qp;
s = qp./p;
if( tot+s>tot )
continue;
end;
% .......... SET ERROR -- IRREGULAR END OF ITERATION.
% DEFLATE MINIMUM DIAGONAL ELEMENT ..........
ierr = fix(5.*n + k);
s = 0.0e0;
delta = qp;
%
for j = k : n;
if( w(j)<=delta )
i = fix(j);
delta = w(j);
end;
end; j = fix(n+1);
elseif( delta<(-eps1) ) ;
ierr = fix(6.*n + 1);
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;
bd_shape=zeros(bd_shape);bd_shape(:)=bd(1:numel(bd_shape));bd=bd_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
break;
end;
gt(3)=0;
% .......... CONVERGENCE ..........
if( i<n )
bd(i+1) = bd(i).*f./qp;
end;
ii = fix(ind(i));
if( i~=k )
k1 = fix(i - k);
% .......... FOR J=I-1 STEP -1 UNTIL K DO -- ..........
for jj = 1 : k1;
j = fix(i - jj);
w(j+1) = w(j) - s;
bd(j+1) = bd(j);
ind(j+1) = fix(ind(j));
end; jj = fix(k1+1);
end;
%
w(k) = tot;
err = err + abs(delta);
bd(k) = err;
ind(k) = fix(ii);
end;
%
if( type )
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;
bd_shape=zeros(bd_shape);bd_shape(:)=bd(1:numel(bd_shape));bd=bd_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
end;
f = bd(1);
e2(1) = 2.0e0;
bd(1) = f;
j = 2;
end;
gt(4)=0;
% .......... NEGATE ELEMENTS OF W FOR LARGEST VALUES ..........
for i = 1 : n;
w(i) = -w(i);
end; i = fix(n+1);
%
jdef = fix(-jdef);
if( j==1 )
gt(1)=1;
continue;
end;
if( j~=2 )
ierr = fix(6.*n + 1);
end;
break;
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;
bd_shape=zeros(bd_shape);bd_shape(:)=bd(1:numel(bd_shape));bd=bd_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
end %subroutine ratqr
%DECK RC3JJ
|
|
