[nm,n,a,b,dl,ierr]=reduc(nm,n,a,b,dl,ierr);

function [nm,n,a,b,dl,ierr]=reduc(nm,n,a,b,dl,ierr);
%***BEGIN PROLOGUE  REDUC
%***PURPOSE  Reduce a generalized symmetric eigenproblem to a standard
%            symmetric eigenproblem using Cholesky factorization.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C1C
%***TYPE      SINGLE PRECISION (REDUC-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure REDUC1,
%     NUM. MATH. 11, 99-110(1968) by Martin and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
%
%     This subroutine reduces the generalized SYMMETRIC eigenproblem
%     Ax=(LAMBDA)Bx, where B is POSITIVE DEFINITE, to the standard
%     symmetric eigenproblem using the Cholesky factorization of B.
%
%     On Input
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameters, A and B, as declared in the calling
%          program dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrices A and B.  If the Cholesky
%          factor L of B is already available, N should be prefixed
%          with a minus sign.  N is an INTEGER variable.
%
%        A and B contain the real symmetric input matrices.  Only
%          the full upper triangles of the matrices need be supplied.
%          If N is negative, the strict lower triangle of B contains,
%          instead, the strict lower triangle of its Cholesky factor L.
%          A and B are two-dimensional REAL arrays, dimensioned A(NM,N)
%          and B(NM,N).
%
%       DL contains, if N is negative, the diagonal elements of L.
%          DL is a one-dimensional REAL array, dimensioned DL(N).
%
%     On Output
%
%        A contains in its full lower triangle the full lower triangle
%          of the symmetric matrix derived from the reduction to the
%          standard form.  The strict upper triangle of A is unaltered.
%
%        B contains in its strict lower triangle the strict lower
%          triangle of its Cholesky factor L.  The full upper triangle
%          of B is unaltered.
%
%        DL contains the diagonal elements of L.
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          7*N+1      if B is not positive definite.
%
%     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  REDUC
%
persistent i i1 j j1 k nn x y ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(i1), i1=0; end;
if isempty(j1), j1=0; end;
if isempty(nn), nn=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([nm])).*prod([nm])-numel(b))],nm,[]);
dl_shape=size(dl);dl=reshape(dl,1,[]);
if isempty(x), x=0; end;
if isempty(y), y=0; end;
%
%***FIRST EXECUTABLE STATEMENT  REDUC
ierr = 0;
nn = fix(abs(n));
if( n>=0 )
%     .......... FORM L IN THE ARRAYS B AND DL ..........
for i = 1 : n;
i1 = fix(i - 1);
%
for j = i : n;
x = b(i,j);
if( i~=1 )
%
for k = 1 : i1;
x = x - b(i,k).*b(j,k);
end; k = fix(i1+1);
end;
%
if( j~=i )
b(j,i) = x./y;
elseif( x<=0.0e0 ) ;
%     .......... SET ERROR -- B IS NOT POSITIVE DEFINITE ..........
ierr = fix(7.*n + 1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
dl_shape=zeros(dl_shape);dl_shape(:)=dl(1:numel(dl_shape));dl=dl_shape;
return;
else;
y = sqrt(x);
dl(i) = y;
end;
end; j = fix(n+1);
end; i = fix(n+1);
end;
%     .......... FORM THE TRANSPOSE OF THE UPPER TRIANGLE OF INV(L)*A
%                IN THE LOWER TRIANGLE OF THE ARRAY A ..........
for i = 1 : nn;
i1 = fix(i - 1);
y = dl(i);
%
for j = i : nn;
x = a(i,j);
if( i~=1 )
%
for k = 1 : i1;
x = x - b(i,k).*a(j,k);
end; k = fix(i1+1);
end;
%
a(j,i) = x./y;
end; j = fix(nn+1);
end; i = fix(nn+1);
%     .......... PRE-MULTIPLY BY INV(L) AND OVERWRITE ..........
for j = 1 : nn;
j1 = fix(j - 1);
%
for i = j : nn;
x = a(i,j);
if( i~=j )
i1 = fix(i - 1);
%
for k = j : i1;
x = x - a(k,j).*b(i,k);
end; k = fix(i1+1);
end;
%
if( j~=1 )
%
for k = 1 : j1;
x = x - a(j,k).*b(i,k);
end; k = fix(j1+1);
end;
%
a(i,j) = x./dl(i);
end; i = fix(nn+1);
%
end; j = fix(nn+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
dl_shape=zeros(dl_shape);dl_shape(:)=dl(1:numel(dl_shape));dl=dl_shape;
end
%DECK REORT