[nm,n,low,igh,a,ort]=orthes(nm,n,low,igh,a,ort);

function [nm,n,low,igh,a,ort]=orthes(nm,n,low,igh,a,ort);
%***BEGIN PROLOGUE  ORTHES
%***PURPOSE  Reduce a real general matrix to upper Hessenberg form
%            using orthogonal similarity transformations.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C1B2
%***TYPE      SINGLE PRECISION (ORTHES-S, CORTH-C)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure ORTHES,
%     NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
%
%     Given a REAL GENERAL matrix, this subroutine
%     reduces a submatrix situated in rows and columns
%     LOW through IGH to upper Hessenberg form by
%     orthogonal similarity transformations.
%
%     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.
%
%        LOW and IGH are two INTEGER variables determined by the
%          balancing subroutine  BALANC.  If  BALANC  has not been
%          used, set LOW=1 and IGH equal to the order of the matrix, N.
%
%        A contains the general matrix to be reduced to upper
%          Hessenberg form.  A is a two-dimensional REAL array,
%          dimensioned A(NM,N).
%
%     On OUTPUT
%
%        A contains the upper Hessenberg matrix.  Some information about
%          the orthogonal transformations used in the reduction
%          is stored in the remaining triangle under the Hessenberg
%          matrix.
%
%        ORT contains further information about the orthogonal trans-
%          formations used in the reduction.  Only elements LOW+1
%          through IGH are used.  ORT is a one-dimensional REAL array,
%          dimensioned ORT(IGH).
%
%     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  ORTHES
%
persistent f g h i ii j jj kp1 la m mp scalemlv ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(m), m=0; end;
if isempty(ii), ii=0; end;
if isempty(jj), jj=0; end;
if isempty(la), la=0; end;
if isempty(mp), mp=0; end;
if isempty(kp1), kp1=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
ort_shape=size(ort);ort=reshape(ort,1,[]);
if isempty(f), f=0; end;
if isempty(g), g=0; end;
if isempty(h), h=0; end;
if isempty(scalemlv), scalemlv=0; end;
%
%***FIRST EXECUTABLE STATEMENT  ORTHES
la = fix(igh - 1);
kp1 = fix(low + 1);
if( la>=kp1 )
%
for m = kp1 : la;
h = 0.0e0;
ort(m) = 0.0e0;
scalemlv = 0.0e0;
%     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
for i = m : igh;
scalemlv = scalemlv + abs(a(i,m-1));
end; i = fix(igh+1);
%
if( scalemlv~=0.0e0 )
mp = fix(m + igh);
%     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
for ii = m : igh;
i = fix(mp - ii);
ort(i) = a(i,m-1)./scalemlv;
h = h + ort(i).*ort(i);
end; ii = fix(igh+1);
%
g = -(abs(sqrt(h)).*sign(ort(m)));
h = h - ort(m).*g;
ort(m) = ort(m) - g;
%     .......... FORM (I-(U*UT)/H) * A ..........
for j = m : n;
f = 0.0e0;
%     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
for ii = m : igh;
i = fix(mp - ii);
f = f + ort(i).*a(i,j);
end; ii = fix(igh+1);
%
f = f./h;
%
for i = m : igh;
a(i,j) = a(i,j) - f.*ort(i);
end; i = fix(igh+1);
%
end; j = fix(n+1);
%     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
for i = 1 : igh;
f = 0.0e0;
%     .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
for jj = m : igh;
j = fix(mp - jj);
f = f + ort(j).*a(i,j);
end; jj = fix(igh+1);
%
f = f./h;
%
for j = m : igh;
a(i,j) = a(i,j) - f.*ort(j);
end; j = fix(igh+1);
%
end; i = fix(igh+1);
%
ort(m) = scalemlv.*ort(m);
a(m,m-1) = scalemlv.*g;
end;
end; m = fix(la+1);
end;
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ort_shape=zeros(ort_shape);ort_shape(:)=ort(1:numel(ort_shape));ort=ort_shape;
end
%DECK ORTHO4