[nm,n,mb,a,w,matz,z,fv1,fv2,ierr]=rsb(nm,n,mb,a,w,matz,z,fv1,fv2,ierr);

function [nm,n,mb,a,w,matz,z,fv1,fv2,ierr]=rsb(nm,n,mb,a,w,matz,z,fv1,fv2,ierr);
%***BEGIN PROLOGUE  RSB
%***PURPOSE  Compute the eigenvalues and, optionally, the eigenvectors
%            of a symmetric band matrix.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4A6
%***TYPE      SINGLE PRECISION (RSB-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine calls the recommended sequence of
%     subroutines from the eigensystem subroutine package (EISPACK)
%     to find the eigenvalues and eigenvectors (if desired)
%     of a REAL SYMMETRIC BAND 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.
%
%        MB is the half band width of the matrix, 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.  MB must be less than or
%          equal to N.  MB is an INTEGER variable.
%
%        A contains the lower triangle of the real symmetric band
%          matrix.  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 the last column.  Contents of storage
%          locations not part of the matrix are arbitrary.  A is a
%          two-dimensional REAL array, dimensioned A(NM,MB).
%
%        MATZ is an INTEGER variable set equal to zero if only
%          eigenvalues are desired.  Otherwise, it is set to any
%          non-zero integer for both eigenvalues and eigenvectors.
%
%     On Output
%
%        A has been destroyed.
%
%        W contains the eigenvalues in ascending order.  W is a one-
%          dimensional REAL array, dimensioned W(N).
%
%        Z contains the eigenvectors if MATZ is not zero.  The
%          eigenvectors are orthonormal.  Z is a two-dimensional
%          REAL array, dimensioned Z(NM,N).
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          10*N       if N is greater than NM,
%          12*N       if MB is either non-positive or greater than N,
%          J          if the J-th eigenvalue has not been
%                     determined after 30 iterations.
%                     The eigenvalues and eigenvectors, if requested,
%                     should be correct for indices 1, 2, ..., IERR-1.
%
%        FV1 and FV2 are one-dimensional REAL arrays used for temporary
%          storage, dimensioned FV1(N) and FV2(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  BANDR, TQL2, TQLRAT
%***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  RSB
%
persistent tf ; 

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,[]);
fv1_shape=size(fv1);fv1=reshape(fv1,1,[]);
fv2_shape=size(fv2);fv2=reshape(fv2,1,[]);
if isempty(tf), tf=false; end;
%
%***FIRST EXECUTABLE STATEMENT  RSB
if( n>nm )
ierr = fix(10.*n);
elseif( mb<=0 ) ;
ierr = fix(12.*n);
elseif( mb>n ) ;
ierr = fix(12.*n);
%
elseif( matz~=0 ) ;
%     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
tf = true;
fv1_orig=fv1;    [nm,n,mb,a,w,fv1,dumvar7,tf,z]=bandr(nm,n,mb,a,w,fv1,fv1,tf,z);    fv1(dumvar7~=fv1_orig)=dumvar7(dumvar7~=fv1_orig);
[nm,n,w,fv1,z,ierr]=tql2(nm,n,w,fv1,z,ierr);
else;
%     .......... FIND EIGENVALUES ONLY ..........
tf = false;
[nm,n,mb,a,w,fv1,fv2,tf,z]=bandr(nm,n,mb,a,w,fv1,fv2,tf,z);
[n,w,fv2,ierr]=tqlrat(n,w,fv2,ierr);
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;
fv1_shape=zeros(fv1_shape);fv1_shape(:)=fv1(1:numel(fv1_shape));fv1=fv1_shape;
fv2_shape=zeros(fv2_shape);fv2_shape(:)=fv2(1:numel(fv2_shape));fv2=fv2_shape;
end
%DECK RSCO