| [nm,n,mb,a,t,r,ierr,nv,rv]=bqr(nm,n,mb,a,t,r,ierr,nv,rv); |
function [nm,n,mb,a,t,r,ierr,nv,rv]=bqr(nm,n,mb,a,t,r,ierr,nv,rv);
%***BEGIN PROLOGUE BQR
%***PURPOSE Compute some of the eigenvalues of a real symmetric
% matrix using the QR method with shifts of origin.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4A6
%***TYPE SINGLE PRECISION (BQR-S)
%***KEYWORDS EIGENVALUES, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure BQR,
% NUM. MATH. 16, 85-92(1970) by Martin, Reinsch, and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971).
%
% This subroutine finds the eigenvalue of smallest (usually)
% magnitude of a REAL SYMMETRIC BAND matrix using the
% QR algorithm with shifts of origin. Consecutive calls
% can be made to find further eigenvalues.
%
% On INPUT
%
% NM must be set to the row dimension of the two-dimensional
% array parameter, A, as declared in the calling program
% dimension statement. NM is an INTEGER variable.
%
% N is the order of the matrix A. N is an INTEGER variable.
% N must be less than or equal to NM.
%
% MB is the (half) band width of the matrix, defined as the
% number of adjacent diagonals, including the principal
% diagonal, required to specify the non-zero portion of the
% lower triangle of the matrix. MB is an INTEGER variable.
% MB must be less than or equal to N on first call.
%
% A contains the lower triangle of the symmetric band input
% matrix stored as an N by MB array. Its lowest subdiagonal
% is stored in the last N+1-MB positions of the first column,
% its next subdiagonal in the last N+2-MB positions of the
% second column, further subdiagonals similarly, and finally
% its principal diagonal in the N positions of the last column.
% Contents of storages not part of the matrix are arbitrary.
% On a subsequent call, its output contents from the previous
% call should be passed. A is a two-dimensional REAL array,
% dimensioned A(NM,MB).
%
% T specifies the shift (of eigenvalues) applied to the diagonal
% of A in forming the input matrix. What is actually determined
% is the eigenvalue of A+TI (I is the identity matrix) nearest
% to T. On a subsequent call, the output value of T from the
% previous call should be passed if the next nearest eigenvalue
% is sought. T is a REAL variable.
%
% R should be specified as zero on the first call, and as its
% output value from the previous call on a subsequent call.
% It is used to determine when the last row and column of
% the transformed band matrix can be regarded as negligible.
% R is a REAL variable.
%
% NV must be set to the dimension of the array parameter RV
% as declared in the calling program dimension statement.
% NV is an INTEGER variable.
%
% On OUTPUT
%
% A contains the transformed band matrix. The matrix A+TI
% derived from the output parameters is similar to the
% input A+TI to within rounding errors. Its last row and
% column are null (if IERR is zero).
%
% T contains the computed eigenvalue of A+TI (if IERR is zero),
% where I is the identity matrix.
%
% R contains the maximum of its input value and the norm of the
% last column of the input matrix A.
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% J if the J-th eigenvalue has not been
% determined after a total of 30 iterations.
%
% RV is a one-dimensional REAL array of dimension NV which is
% at least (2*MB**2+4*MB-3), used for temporary storage. The
% first (3*MB-2) locations correspond to the ALGOL array B,
% the next (2*MB-1) locations correspond to the ALGOL array H,
% and the final (2*MB**2-MB) locations correspond to the MB
% by (2*MB-1) ALGOL array U.
%
% NOTE. For a subsequent call, N should be replaced by N-1, but
% MB should not be altered even when it exceeds the current N.
%
% Calls PYTHAG(A,B) for SQRT(A**2 + B**2).
%
% 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 PYTHAG
%***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 BQR
%
persistent f g i igo ii ik imult its j jk jm k kj kj1 kk km l ll m m1 m2 m21 m3 m31 m4 mk mn mz ni q s scalemlv ;
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(m), m=0; end;
if isempty(ii), ii=0; end;
if isempty(ik), ik=0; end;
if isempty(jk), jk=0; end;
if isempty(jm), jm=0; end;
if isempty(kj), kj=0; end;
if isempty(kk), kk=0; end;
if isempty(km), km=0; end;
if isempty(ll), ll=0; end;
if isempty(mk), mk=0; end;
if isempty(mn), mn=0; end;
if isempty(mz), mz=0; end;
if isempty(igo), igo=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(m3), m3=0; end;
if isempty(m4), m4=0; end;
if isempty(ni), ni=0; end;
if isempty(its), its=0; end;
if isempty(kj1), kj1=0; end;
if isempty(m21), m21=0; end;
if isempty(m31), m31=0; end;
if isempty(imult), imult=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
rv_shape=size(rv);rv=reshape(rv,1,[]);
if isempty(f), f=0; end;
if isempty(g), g=0; end;
if isempty(q), q=0; end;
if isempty(s), s=0; end;
if isempty(scalemlv), scalemlv=0; end;
%
%***FIRST EXECUTABLE STATEMENT BQR
ierr = 0;
m1 = fix(min(mb,n));
m = fix(m1 - 1);
m2 = fix(m + m);
m21 = fix(m2 + 1);
m3 = fix(m21 + m);
m31 = fix(m3 + 1);
m4 = fix(m31 + m2);
mn = fix(m + n);
mz = fix(mb - m1);
its = 0;
% .......... TEST FOR CONVERGENCE ..........
while (1);
g = a(n,mb);
if( m~=0 )
f = 0.0e0;
%
for k = 1 : m;
mk = fix(k + mz);
f = f + abs(a(n,mk));
end; k = fix(m+1);
%
if( its==0 && f>r )
r = f;
end;
if( r+f>r )
if( its==30 )
% .......... SET ERROR -- NO CONVERGENCE TO
% EIGENVALUE AFTER 30 ITERATIONS ..........
ierr = fix(n);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
rv_shape=zeros(rv_shape);rv_shape(:)=rv(1:numel(rv_shape));rv=rv_shape;
return;
else;
its = fix(its + 1);
% .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR ..........
if( f<=0.25e0.*r || its>=5 )
f = a(n,mb-1);
if( f~=0.0e0 )
q =(a(n-1,mb)-g)./(2.0e0.*f);
[s ,q]=pythag(q,1.0e0);
g = g - f./(q+(abs(s).*sign(q)));
end;
t = t + g;
%
for i = 1 : n;
a(i,mb) = a(i,mb) - g;
end; i = fix(n+1);
end;
%
for k = m31 : m4;
rv(k) = 0.0e0;
end; k = fix(m4+1);
%
for ii = 1 : mn;
i = fix(ii - m);
ni = fix(n - ii);
igo=1;
if( ni<0 )
igo=0;
else;
% .......... FORM COLUMN OF SHIFTED MATRIX A-G*I ..........
l = fix(max(1,2-i));
%
for k = 1 : m3;
rv(k) = 0.0e0;
end; k = fix(m3+1);
%
for k = l : m1;
km = fix(k + m);
mk = fix(k + mz);
rv(km) = a(ii,mk);
end; k = fix(m1+1);
%
ll = fix(min(m,ni));
if( ll~=0 )
%
for k = 1 : ll;
km = fix(k + m21);
ik = fix(ii + k);
mk = fix(mb - k);
rv(km) = a(ik,mk);
end; k = fix(ll+1);
end;
% .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
ll = fix(m2);
imult = 0;
end;
% .......... MULTIPLICATION PROCEDURE ..........
while (1);
if(igo==1)
kj = fix(m4 - m1);
%
for j = 1 : ll;
kj = fix(kj + m1);
jm = fix(j + m3);
if( rv(jm)~=0.0e0 )
f = 0.0e0;
%
for k = 1 : m1;
kj = fix(kj + 1);
jk = fix(j + k - 1);
f = f + rv(kj).*rv(jk);
end; k = fix(m1+1);
%
f = f./rv(jm);
kj = fix(kj - m1);
%
for k = 1 : m1;
kj = fix(kj + 1);
jk = fix(j + k - 1);
rv(jk) = rv(jk) - rv(kj).*f;
end; k = fix(m1+1);
%
kj = fix(kj - m1);
end;
end; j = fix(ll+1);
%
if( imult~=0 )
% .......... STORE COLUMN OF NEW A MATRIX ..........
for k = l : m1;
mk = fix(k + mz);
a(i,mk) = rv(k);
end; k = fix(m1+1);
%to 115
break;
else;
% .......... HOUSEHOLDER REFLECTION ..........
f = rv(m21);
s = 0.0e0;
rv(m4) = 0.0e0;
scalemlv = 0.0e0;
%
for k = m21 : m3;
scalemlv = scalemlv + abs(rv(k));
end; k = fix(m3+1);
%
if( scalemlv~=0.0e0 )
%
for k = m21 : m3;
s = s +(rv(k)./scalemlv).^2;
end; k = fix(m3+1);
%
s = scalemlv.*scalemlv.*s;
g = -(abs(sqrt(s)).*sign(f));
rv(m21) = g;
rv(m4) = s - f.*g;
kj = fix(m4 + m2.*m1 + 1);
rv(kj) = f - g;
%
for k = 2 : m1;
kj = fix(kj + 1);
km = fix(k + m2);
rv(kj) = rv(km);
end; k = fix(m1+1);
end;
% .......... SAVE COLUMN OF TRIANGULAR FACTOR R ..........
for k = l : m1;
km = fix(k + m);
mk = fix(k + mz);
a(ii,mk) = rv(km);
end; k = fix(m1+1);
end;
%to 110
end;
%
l = fix(max(1,m1+1-i));
if( i>0 )
% .......... PERFORM ADDITIONAL STEPS ..........
for k = 1 : m21;
rv(k) = 0.0e0;
end; k = fix(m21+1);
%
ll = fix(min(m1,ni+m1));
% .......... GET ROW OF TRIANGULAR FACTOR R ..........
for kk = 1 : ll;
k = fix(kk - 1);
km = fix(k + m1);
ik = fix(i + k);
mk = fix(mb - k);
rv(km) = a(ik,mk);
end; kk = fix(ll+1);
% .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
ll = fix(m1);
imult = 1;
%back to 105
continue;
end;
%out of 105
break;
end;
% .......... UPDATE HOUSEHOLDER REFLECTIONS ..........
if( l>1 )
l = fix(l - 1);
end;
kj1 = fix(m4 + l.*m1);
%
for j = l : m2;
jm = fix(j + m3);
rv(jm) = rv(jm+1);
%
for k = 1 : m1;
kj1 = fix(kj1 + 1);
kj = fix(kj1 - m1);
rv(kj) = rv(kj1);
end; k = fix(m1+1);
end; j = fix(m2+1);
%
end;
%
%back to 100
continue;
end;
end;
end;
break;
end;
% .......... CONVERGENCE ..........
t = t + g;
%
for i = 1 : n;
a(i,mb) = a(i,mb) - g;
end; i = fix(n+1);
%
for k = 1 : m1;
mk = fix(k + mz);
a(n,mk) = 0.0e0;
%
end; k = fix(m1+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
rv_shape=zeros(rv_shape);rv_shape(:)=rv(1:numel(rv_shape));rv=rv_shape;
end %subroutine bqr
%DECK BSGQ8
|
|
