[n,nv,a,d,e,e2]=tred3(n,nv,a,d,e,e2);

function [n,nv,a,d,e,e2]=tred3(n,nv,a,d,e,e2);
%***BEGIN PROLOGUE  TRED3
%***PURPOSE  Reduce a real symmetric matrix stored in packed form to
%            symmetric tridiagonal matrix using orthogonal
%            transformations.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C1B1
%***TYPE      SINGLE PRECISION (TRED3-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure TRED3,
%     NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
%
%     This subroutine reduces a REAL SYMMETRIC matrix, stored as
%     a one-dimensional array, to a symmetric tridiagonal matrix
%     using orthogonal similarity transformations.
%
%     On Input
%
%        N is the order of the matrix A.  N is an INTEGER variable.
%
%        NV is an INTEGER variable set equal to the dimension of the
%          array A as specified in the calling program.  NV must not
%          be less than  N*(N+1)/2.
%
%        A contains the lower triangle, stored row-wise, of the real
%          symmetric packed matrix.  A is a one-dimensional REAL
%          array, dimensioned A(NV).
%
%     On Output
%
%        A contains information about the orthogonal transformations
%          used in the reduction in its first N*(N+1)/2 positions.
%
%        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 set
%          to zero.  E is a one-dimensional REAL array, dimensioned
%          E(N).
%
%        E2 contains the squares of the corresponding elements of E.
%          E2 may coincide with E if the squares are not needed.
%          E2 is a one-dimensional REAL array, dimensioned E2(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
%   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  TRED3
%
persistent f g h hh i ii iz j jk k l 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(ii), ii=0; end;
if isempty(iz), iz=0; end;
if isempty(jk), jk=0; end;
a_shape=size(a);a=reshape(a,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,[]);
if isempty(f), f=0; end;
if isempty(g), g=0; end;
if isempty(h), h=0; end;
if isempty(hh), hh=0; end;
if isempty(scalemlv), scalemlv=0; end;
%
%     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
%***FIRST EXECUTABLE STATEMENT  TRED3
for ii = 1 : n;
i = fix(n + 1 - ii);
l = fix(i - 1);
iz =fix(fix((i.*l)./2));
h = 0.0e0;
scalemlv = 0.0e0;
if( l>=1 )
%     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
for k = 1 : l;
iz = fix(iz + 1);
d(k) = a(iz);
scalemlv = scalemlv + abs(d(k));
end; k = fix(l+1);
%
if( scalemlv~=0.0e0 )
%
for k = 1 : l;
d(k) = d(k)./scalemlv;
h = h + d(k).*d(k);
end; k = fix(l+1);
%
e2(i) = scalemlv.*scalemlv.*h;
f = d(l);
g = -(abs(sqrt(h)).*sign(f));
e(i) = scalemlv.*g;
h = h - f.*g;
d(l) = f - g;
a(iz) = scalemlv.*d(l);
if( l~=1 )
f = 0.0e0;
%
for j = 1 : l;
g = 0.0e0;
jk =fix(fix((j.*(j-1))./2));
%     .......... FORM ELEMENT OF A*U ..........
for k = 1 : l;
jk = fix(jk + 1);
if( k>j )
jk = fix(jk + k - 2);
end;
g = g + a(jk).*d(k);
end; k = fix(l+1);
%     .......... FORM ELEMENT OF P ..........
e(j) = g./h;
f = f + e(j).*d(j);
end; j = fix(l+1);
%
hh = f./(h+h);
jk = 0;
%     .......... FORM REDUCED A ..........
for j = 1 : l;
f = d(j);
g = e(j) - hh.*f;
e(j) = g;
%
for k = 1 : j;
jk = fix(jk + 1);
a(jk) = a(jk) - f.*e(k) - g.*d(k);
end; k = fix(j+1);
end; j = fix(l+1);
end;
d(i) = a(iz+1);
a(iz+1) = scalemlv.*sqrt(h);
continue;
end;
end;
e(i) = 0.0e0;
e2(i) = 0.0e0;
%
d(i) = a(iz+1);
a(iz+1) = scalemlv.*sqrt(h);
end; ii = fix(n+1);
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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;
end %subroutine tred3
%DECK TRI3