[snrm2result,n,sx,incx]=snrm2(n,sx,incx);

function [snrm2result,n,sx,incx]=snrm2(n,sx,incx);
snrm2result=[];
persistent cuthi cutlo firstCall gt20 gt40 gt60 gt80 hitest i j next nn one summlv xmax zero ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(nn), nn=0; end;
%***BEGIN PROLOGUE  SNRM2
%***PURPOSE  Compute the Euclidean length (L2 norm) of a vector.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1A3B
%***TYPE      SINGLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
%***KEYWORDS  BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
%             LINEAR ALGEBRA, UNITARY, VECTOR
%***AUTHOR  Lawson, C. L., (JPL)
%           Hanson, R. J., (SNLA)
%           Kincaid, D. R., (U. of Texas)
%           Krogh, F. T., (JPL)
%***DESCRIPTION
%
%                B L A S  Subprogram
%    Description of Parameters
%
%     --Input--
%        N  number of elements in input vector(s)
%       SX  single precision vector with N elements
%     INCX  storage spacing between elements of SX
%
%     --Output--
%    SNRM2  single precision result (zero if N .LE. 0)
%
%     Euclidean norm of the N-vector stored in SX with storage
%     increment INCX .
%     If N .LE. 0, return with result = 0.
%     If N .GE. 1, then INCX must be .GE. 1
%
%     Four Phase Method using two built-in constants that are
%     hopefully applicable to all machines.
%         CUTLO = maximum of  SQRT(U/EPS)  over all known machines.
%         CUTHI = minimum of  SQRT(V)      over all known machines.
%     where
%         EPS = smallest no. such that EPS + 1. .GT. 1.
%         U   = smallest positive no.   (underflow limit)
%         V   = largest  no.            (overflow  limit)
%
%     Brief Outline of Algorithm.
%
%     Phase 1 scans zero components.
%     Move to phase 2 when a component is nonzero and .LE. CUTLO
%     Move to phase 3 when a component is .GT. CUTLO
%     Move to phase 4 when a component is .GE. CUTHI/M
%     where M = N for X() real and M = 2*N for complex.
%
%     Values for CUTLO and CUTHI.
%     From the environmental parameters listed in the IMSL converter
%     document the limiting values are as follows:
%     CUTLO, S.P.   U/EPS = 2**(-102) for  Honeywell.  Close seconds are
%                   Univac and DEC at 2**(-103)
%                   Thus CUTLO = 2**(-51) = 4.44089E-16
%     CUTHI, S.P.   V = 2**127 for Univac, Honeywell, and DEC.
%                   Thus CUTHI = 2**(63.5) = 1.30438E19
%     CUTLO, D.P.   U/EPS = 2**(-67) for Honeywell and DEC.
%                   Thus CUTLO = 2**(-33.5) = 8.23181D-11
%     CUTHI, D.P.   same as S.P.  CUTHI = 1.30438D19
%     DATA CUTLO, CUTHI /8.232D-11,  1.304D19/
%     DATA CUTLO, CUTHI /4.441E-16,  1.304E19/
%
%***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
%                 Krogh, Basic linear algebra subprograms for Fortran
%                 usage, Algorithm No. 539, Transactions on Mathematical
%                 Software 5, 3 (September 1979), pp. 308-323.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   791001  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  SNRM2
if isempty(next), next=0; end;
if isempty(gt20), gt20=0; end;
if isempty(gt40), gt40=0; end;
if isempty(gt60), gt60=0; end;
if isempty(gt80), gt80=0; end;
sx_shape=size(sx);sx=reshape(sx,1,[]);
if isempty(cutlo), cutlo=0; end;
if isempty(cuthi), cuthi=0; end;
if isempty(hitest), hitest=0; end;
if isempty(summlv), summlv=0; end;
if isempty(xmax), xmax=0; end;
if isempty(zero), zero=0; end;
if isempty(one), one=0; end;
if firstCall,   zero =[0.0e0];  end;
if firstCall,  one=[1.0e0];  end;
%
if firstCall,   cutlo =[4.441e-16];  end;
if firstCall,  cuthi=[1.304e19];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  SNRM2
if( n>0 )
%
next = 200;
summlv = zero;
nn = fix(n.*incx);
%
%                                                 BEGIN MAIN LOOP
%
i = 1;
gt20=0;
gt40=0;
gt60=0;
gt80=0;
while( true );
if(gt80==0)
if(gt60==0)
if(gt40==0)
if(gt20==0)
if( next~=200 )
if( next==300 )
gt40=1;
continue;
end;
if( next==400 )
%
%                   PHASE 2.  SUM IS SMALL.
%                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
%
if( abs(sx(i))>cutlo )
%
%                  PREPARE FOR PHASE 3.
%
summlv =(summlv.*xmax).*xmax;
gt60 =1;
continue;
end;
elseif( next~=500 ) ;
gt20 =1;
continue;
end;
%
%                     COMMON CODE FOR PHASES 2 AND 4.
%                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
%
if( abs(sx(i))<=xmax )
summlv = summlv +(sx(i)./xmax).^2;
else;
summlv = one + summlv.*(xmax./sx(i)).^2;
xmax = abs(sx(i));
end;
gt80 =1;
continue;
end;
end;
gt20=0;
if( abs(sx(i))>cutlo )
gt60 =1;
continue;
end;
next = 300;
xmax = zero;
end;
gt40=0;
%
%                        PHASE 1.  SUM IS ZERO
%
if( sx(i)==zero )
gt80 =1;
continue;
end;
if( abs(sx(i))<=cutlo )
%
%                                PREPARE FOR PHASE 2.
%
next = 400;
xmax = abs(sx(i));
%
summlv = summlv +(sx(i)./xmax).^2;
gt80 =1;
continue;
end;
end;
gt60=0;
%
%     FOR REAL OR D.P. SET HITEST = CUTHI/N
%     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
%
hitest = cuthi./n;
%
%                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
%
for j = i : incx: nn ;
if( abs(sx(j))>=hitest )
%
%                                PREPARE FOR PHASE 4.
%
i = fix(j);
next = 500;
summlv =(summlv./sx(i))./sx(i);
xmax = abs(sx(i));
summlv = summlv +(sx(i)./xmax).^2;
gt80 =1;
break;
else;
summlv = summlv + sx(j).^2;
end;
end;
if(gt80==1)
continue;
end;
break;
end;
gt80=0;
%
i = fix(i + incx);
if( i>nn )
%
%              end OF MAIN LOOP.
%
%              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
%
snrm2result = xmax.*sqrt(summlv);
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',incx); evalin('caller',[inputname(3),'=FUntemp;']); end
return;
end;
end;
snrm2result = sqrt(summlv);
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',incx); evalin('caller',[inputname(3),'=FUntemp;']); end
return;
else;
snrm2result = zero;
end;
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',incx); evalin('caller',[inputname(3),'=FUntemp;']); end
end %function snrm2
%DECK SNSQE