| [nm,n,d,e,e2,m,w,ind,z,ierr,rv1,rv2,rv3,rv4,rv6]=tinvit(nm,n,d,e,e2,m,w,ind,z,ierr,rv1,rv2,rv3,rv4,rv6); |
function [nm,n,d,e,e2,m,w,ind,z,ierr,rv1,rv2,rv3,rv4,rv6]=tinvit(nm,n,d,e,e2,m,w,ind,z,ierr,rv1,rv2,rv3,rv4,rv6);
%***BEGIN PROLOGUE TINVIT
%***PURPOSE Compute the eigenvectors of symmetric tridiagonal matrix
% corresponding to specified eigenvalues, using inverse
% iteration.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C3
%***TYPE SINGLE PRECISION (TINVIT-S)
%***KEYWORDS EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the inverse iteration tech-
% nique in the ALGOL procedure TRISTURM by Peters and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
%
% This subroutine finds those eigenvectors of a TRIDIAGONAL
% SYMMETRIC matrix corresponding to specified eigenvalues,
% using 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.
%
% 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,
% with zeros corresponding to negligible elements of E.
% E(I) is considered negligible if it is not larger than
% the product of the relative machine precision and the sum
% of the magnitudes of D(I) and D(I-1). E2(1) must contain
% 0.0e0 if the eigenvalues are in ascending order, or 2.0e0
% if the eigenvalues are in descending order. If BISECT,
% TRIDIB, or IMTQLV has been used to find the eigenvalues,
% their output E2 array is exactly what is expected here.
% E2 is a one-dimensional REAL array, dimensioned E2(N).
%
% M is the number of specified eigenvalues for which eigenvectors
% are to be determined. M is an INTEGER variable.
%
% W contains the M eigenvalues in ascending or descending 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.
% If BISECT or TRIDIB has been used to determine the
% eigenvalues, their output IND array is suitable for input
% to TINVIT. IND is a one-dimensional INTEGER array,
% dimensioned IND(M).
%
% On Output
%
% ** All input arrays are unaltered.**
%
% Z contains the associated set of orthonormal eigenvectors.
% Any vector which fails to converge is set to zero.
% Z is a two-dimensional REAL array, dimensioned Z(NM,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 in 5 iterations.
%
% RV1, RV2 and RV3 are one-dimensional REAL arrays used for
% temporary storage. They are used to store the main diagonal
% and the two adjacent diagonals of the triangular matrix
% produced in the inverse iteration process. RV1, RV2 and
% RV3 are dimensioned RV1(N), RV2(N) and RV3(N).
%
% RV4 and RV6 are one-dimensional REAL arrays used for temporary
% storage. RV4 holds the multipliers of the Gaussian
% elimination process. RV6 holds the approximate eigenvectors
% in this process. RV4 and RV6 are dimensioned RV4(N) and
% 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 TINVIT
%
persistent eps2 eps3 eps4 group i igo ii ip its j jj norm order p q r s tag u uk v x0 x1 xu ;
if isempty(i), i=0; end;
if isempty(j), j=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(its), its=0; end;
if isempty(tag), tag=0; end;
if isempty(group), group=0; end;
if isempty(igo), igo=0; end;
ind_shape=size(ind);ind=reshape(ind,1,[]);
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,[]);
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;
%
%***FIRST EXECUTABLE STATEMENT TINVIT
ierr = 0;
if( m==0 )
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
return;
end;
tag = 0;
order = 1.0e0 - e2(1);
q = 0;
% .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
while( true );
p = fix(q + 1);
%
for q = p : n;
if( q==n )
break;
end;
if( e2(q+1)==0.0e0 )
break;
end;
end;
% .......... FIND VECTORS BY INVERSE ITERATION ..........
tag = fix(tag + 1);
s = 0;
%
igo=0;
for r = 1 : m;
if( ind(r)==tag )
its = 1;
x1 = w(r);
while (1);
if( s==0 )
% .......... CHECK FOR ISOLATED ROOT ..........
xu = 1.0e0;
if( p==q )
% .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
rv6(p) = 1.0e0;
break;
else;
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;
eps3 = norm;
while( true );
eps3 = 0.5e0.*eps3;
if( norm+eps3<=norm )
break;
end;
end;
uk = sqrt(real(q-p+5));
eps3 = uk.*eps3;
eps4 = uk.*eps3;
uk = eps4./uk;
s = fix(p);
group = 0;
end;
elseif( abs(x1-x0)>=eps2 ) ;
group = 0;
else;
group = fix(group + 1);
if( order.*(x1-x0)<=0.0e0 )
x1 = x0 + order.*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;
% .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
% E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY ..........
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 )
j = fix(r);
%
for jj = 1 : group;
while( true );
j = fix(j - 1);
if( ind(j)==tag )
break;
end;
end;
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;
end;
%
norm = 0.0e0;
%
for i = p : q;
norm = norm + abs(rv6(i));
end; i = fix(q+1);
%
if( norm>=1.0e0 )
% .......... 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);
igo=1;
break;
else;
if( its==5 )
break;
end;
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;
if(igo==1)
break;
end;
% .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
ierr = fix(-r);
xu = 0.0e0;
% .......... FORWARD SUBSTITUTION ..........
%
break;
end;
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;
end;
end;
%
if( q>=n )
break;
end;
end;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
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;
rv6_shape=zeros(rv6_shape);rv6_shape(:)=rv6(1:numel(rv6_shape));rv6=rv6_shape;
end %subroutine tinvit
%DECK TQL1
|
|
