| [nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6]=bandv(nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6); |
function [nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6]=bandv(nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6);
%***BEGIN PROLOGUE BANDV
%***PURPOSE Form the eigenvectors of a real symmetric band matrix
% associated with a set of ordered approximate eigenvalues
% by inverse iteration.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C3
%***TYPE SINGLE PRECISION (BANDV-S)
%***KEYWORDS EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine finds those eigenvectors of a REAL SYMMETRIC
% BAND matrix corresponding to specified eigenvalues, using inverse
% iteration. The subroutine may also be used to solve systems
% of linear equations with a symmetric or non-symmetric band
% coefficient matrix.
%
% On INPUT
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, A and Z, 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.
%
% MBW is the number of columns of the array A used to store the
% band matrix. If the matrix is symmetric, MBW is its (half)
% band width, denoted MB and 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. If the subroutine is being used to solve systems
% of linear equations and the coefficient matrix is not
% symmetric, it must however have the same number of adjacent
% diagonals above the main diagonal as below, and in this
% case, MBW=2*MB-1. MBW is an INTEGER variable. MB must not
% be greater than N.
%
% 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 column MB.
% If the subroutine is being used to solve systems of linear
% equations and the coefficient matrix is not symmetric, A is
% N by 2*MB-1 instead with lower triangle as above and with
% its first superdiagonal stored in the first N-1 positions of
% column MB+1, its second superdiagonal in the first N-2
% positions of column MB+2, further superdiagonals similarly,
% and finally its highest superdiagonal in the first N+1-MB
% positions of the last column. Contents of storage locations
% not part of the matrix are arbitrary. A is a two-dimensional
% REAL array, dimensioned A(NM,MBW).
%
% E21 specifies the ordering of the eigenvalues and contains
% 0.0E0 if the eigenvalues are in ascending order, or
% 2.0E0 if the eigenvalues are in descending order.
% If the subroutine is being used to solve systems of linear
% equations, E21 should be set to 1.0E0 if the coefficient
% matrix is symmetric and to -1.0E0 if not. E21 is a REAL
% variable.
%
% M is the number of specified eigenvalues or the number of
% systems of linear equations. M is an INTEGER variable.
%
% W contains the M eigenvalues in ascending or descending order.
% If the subroutine is being used to solve systems of linear
% equations (A-W(J)*I)*X(J)=B(J), where I is the identity
% matrix, W(J) should be set accordingly, for J=1,2,...,M.
% W is a one-dimensional REAL array, dimensioned W(M).
%
% Z contains the constant matrix columns (B(J),J=1,2,...,M), if
% the subroutine is used to solve systems of linear equations.
% Z is a two-dimensional REAL array, dimensioned Z(NM,M).
%
% 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 and W are unaltered.
%
% Z contains the associated set of orthogonal eigenvectors.
% Any vector which fails to converge is set to zero. If the
% subroutine is used to solve systems of linear equations,
% Z contains the solution matrix columns (X(J),J=1,2,...,M).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% -J if the eigenvector corresponding to the J-th
% eigenvalue fails to converge, or if the J-th
% system of linear equations is nearly singular.
%
% RV and RV6 are temporary storage arrays. If the subroutine
% is being used to solve systems of linear equations, the
% determinant (up to sign) of A-W(M)*I is available, upon
% return, as the product of the first N elements of RV.
% RV and RV6 are one-dimensional REAL arrays. Note that RV
% is dimensioned RV(NV), where NV must be at least N*(2*MB-1).
% RV6 is dimensioned RV6(N).
%
% 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 (NONE)
%***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 BANDV
%
persistent eps2 eps3 eps4 group i igo ii ij ij1 its j jj k kj kj1 m1 m21 maxj maxk mb norm order r s u uk v x0 x1 xu ;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(r), r=0; end;
if isempty(ii), ii=0; end;
if isempty(ij), ij=0; end;
if isempty(jj), jj=0; end;
if isempty(kj), kj=0; end;
if isempty(mb), mb=0; end;
if isempty(m1), m1=0; end;
if isempty(ij1), ij1=0; end;
if isempty(its), its=0; end;
if isempty(kj1), kj1=0; end;
if isempty(m21), m21=0; end;
if isempty(igo), igo=0; end;
if isempty(maxj), maxj=0; end;
if isempty(maxk), maxk=0; end;
if isempty(group), group=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
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,[]);
rv_shape=size(rv);rv=reshape(rv,1,[]);
rv6_shape=size(rv6);rv6=reshape(rv6,1,[]);
if isempty(u), u=0; end;
if isempty(v), v=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(order), order=0; end;
if isempty(s), s=0; end;
%
%***FIRST EXECUTABLE STATEMENT BANDV
ierr = 0;
if( m~=0 )
mb = fix(mbw);
if( e21<0.0e0 )
mb =fix(fix((mbw+1)./2));
end;
m1 = fix(mb - 1);
m21 = fix(m1 + mb);
order = 1.0e0 - abs(e21);
% .......... FIND VECTORS BY INVERSE ITERATION ..........
for r = 1 : m;
its = 1;
x1 = w(r);
if( r==1 )
% .......... COMPUTE NORM OF MATRIX ..........
norm = 0.0e0;
%
for j = 1 : mb;
jj = fix(mb + 1 - j);
kj = fix(jj + m1);
ij = 1;
s = 0.0e0;
%
for i = jj : n;
s = s + abs(a(i,j));
if( e21<0.0e0 )
s = s + abs(a(ij,kj));
ij = fix(ij + 1);
end;
end; i = fix(n+1);
%
norm = max(norm,s);
end; j = fix(mb+1);
%
if( e21<0.0e0 )
norm = 0.5e0.*norm;
end;
% .......... 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 ..........
if( norm==0.0e0 )
norm = 1.0e0;
end;
eps2 = 1.0e-3.*norm.*abs(order);
eps3 = norm;
while (1);
eps3 = 0.5e0.*eps3;
if( norm+eps3<=norm )
break;
end;
end;
uk = sqrt(real(n));
eps3 = uk.*eps3;
eps4 = uk.*eps3;
group = 0;
% .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
elseif( abs(x1-x0)>=eps2 ) ;
group = 0;
else;
group = fix(group + 1);
if( order.*(x1-x0)<=0.0e0 )
x1 = x0 + order.*eps3;
end;
end;
% .......... EXPAND MATRIX, SUBTRACT EIGENVALUE,
% AND INITIALIZE VECTOR ..........
for i = 1 : n;
ij = fix(i + min(0,i-m1).*n);
kj = fix(ij + mb.*n);
ij1 = fix(kj + m1.*n);
if( m1~=0 )
%
for j = 1 : m1;
if( ij>m1 )
rv(ij) = a(i,j);
elseif( ij<=0 ) ;
rv(ij1) = 0.0e0;
ij1 = fix(ij1 + n);
end;
ij = fix(ij + n);
ii = fix(i + j);
if( ii<=n )
jj = fix(mb - j);
if( e21<0.0e0 )
ii = fix(i);
jj = fix(mb + j);
end;
rv(kj) = a(ii,jj);
kj = fix(kj + n);
end;
end; j = fix(m1+1);
end;
%
rv(ij) = a(i,mb) - x1;
rv6(i) = eps4;
if( order==0.0e0 )
rv6(i) = z(i,r);
end;
end; i = fix(n+1);
%
if( m1~=0 )
% .......... ELIMINATION WITH INTERCHANGES ..........
for i = 1 : n;
ii = fix(i + 1);
maxk = fix(min(i+m1-1,n));
maxj = fix(min(n-i,m21-2).*n);
%
for k = i : maxk;
kj1 = fix(k);
j = fix(kj1 + n);
jj = fix(j + maxj);
%
for kj = j : n: jj ;
rv(kj1) = rv(kj);
kj1 = fix(kj);
end; kj = fix(jj +1);
%
rv(kj1) = 0.0e0;
end; k = fix(maxk+1);
%
if( i~=n )
u = 0.0e0;
maxk = fix(min(i+m1,n));
maxj = fix(min(n-ii,m21-2).*n);
%
for j = i : maxk;
if( abs(rv(j))>=abs(u) )
u = rv(j);
k = fix(j);
end;
end; j = fix(maxk+1);
%
j = fix(i + n);
jj = fix(j + maxj);
if( k~=i )
kj = fix(k);
%
for ij = i : n: jj ;
v = rv(ij);
rv(ij) = rv(kj);
rv(kj) = v;
kj = fix(kj + n);
end; ij = fix(jj +1);
%
if( order==0.0e0 )
v = rv6(i);
rv6(i) = rv6(k);
rv6(k) = v;
end;
end;
if( u~=0.0e0 )
%
for k = ii : maxk;
v = rv(k)./u;
kj = fix(k);
%
for ij = j : n: jj ;
kj = fix(kj + n);
rv(kj) = rv(kj) - v.*rv(ij);
end; ij = fix(jj +1);
%
if( order==0.0e0 )
rv6(k) = rv6(k) - v.*rv6(i);
end;
end; k = fix(maxk+1);
end;
end;
%
end; i = fix(n+1);
end;
% .......... BACK SUBSTITUTION
% FOR I=N STEP -1 UNTIL 1 DO -- ..........
igo=1;
while(igo==1);
igo=0;
for ii = 1 : n;
i = fix(n + 1 - ii);
maxj = fix(min(ii,m21));
if( maxj~=1 )
ij1 = fix(i);
j = fix(ij1 + n);
jj = fix(j +(maxj-2).*n);
%
for ij = j : n: jj ;
ij1 = fix(ij1 + 1);
rv6(i) = rv6(i) - rv(ij).*rv6(ij1);
end; ij = fix(jj +1);
end;
%
v = rv(i);
if( abs(v)<eps3 )
% .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM ..........
if( order==0.0e0 )
ierr = fix(-r);
end;
v = (abs(eps3).*sign(v));
end;
rv6(i) = rv6(i)./v;
end; ii = fix(n+1);
%
xu = 1.0e0;
if( order~=0.0e0 )
% .......... 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 = 1 : n;
xu = xu + rv6(i).*z(i,j);
end; i = fix(n+1);
%
for i = 1 : n;
rv6(i) = rv6(i) - xu.*z(i,j);
end; i = fix(n+1);
%
end; jj = fix(group+1);
end;
%
norm = 0.0e0;
%
for i = 1 : n;
norm = norm + abs(rv6(i));
end; i = fix(n+1);
%
if( norm>=0.1e0 )
% .......... NORMALIZE SO THAT SUM OF SQUARES IS
% 1 AND EXPAND TO FULL ORDER ..........
u = 0.0e0;
%
for i = 1 : n;
u = u + rv6(i).^2;
end; i = fix(n+1);
%
xu = 1.0e0./sqrt(u);
% .......... IN-LINE PROCEDURE FOR CHOOSING
% A NEW STARTING VECTOR ..........
elseif( its>=n ) ;
% .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
ierr = fix(-r);
xu = 0.0e0;
else;
its = fix(its + 1);
xu = eps4./(uk+1.0e0);
rv6(1) = eps4;
%
for i = 2 : n;
rv6(i) = xu;
end; i = fix(n+1);
%
rv6(its) = rv6(its) - eps4.*uk;
igo=1;
continue;
end;
end;
end;
for i = 1 : n;
z(i,r) = rv6(i).*xu;
end; i = fix(n+1);
%
x0 = x1;
end;
end;
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
rv_shape=zeros(rv_shape);rv_shape(:)=rv(1:numel(rv_shape));rv=rv_shape;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
end
%DECK BCRH
|
|
