[denormresult,n,x]=denorm(n,x);

function [denormresult,n,x]=denorm(n,x);
denormresult=[];
persistent agiant firstCall floatn i one rdwarf rgiant s1 s2 s3 x1max x3max xabs zero ; if isempty(firstCall),firstCall=1;end; 

;
%***BEGIN PROLOGUE  DENORM
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DNSQ and DNSQE
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (ENORM-S, DENORM-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     Given an N-vector X, this function calculates the
%     Euclidean norm of X.
%
%     The Euclidean norm is computed by accumulating the sum of
%     squares in three different sums. The sums of squares for the
%     small and large components are scaled so that no overflows
%     occur. Non-destructive underflows are permitted. Underflows
%     and overflows do not occur in the computation of the unscaled
%     sum of squares for the intermediate components.
%     The definitions of small, intermediate and large components
%     depend on two constants, RDWARF and RGIANT. The main
%     restrictions on these constants are that RDWARF**2 not
%     underflow and RGIANT**2 not overflow. The constants
%     given here are suitable for every known computer.
%
%     The function statement is
%
%       doubleprecision function DENORM(N,X)
%
%     where
%
%       N is a positive integer input variable.
%
%       X is an input array of length N.
%
%***SEE ALSO  DNSQ, DNSQE
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800301  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DENORM
if isempty(i), i=0; end;
if isempty(agiant), agiant=0; end;
if isempty(floatn), floatn=0; end;
if isempty(one), one=0; end;
if isempty(rdwarf), rdwarf=0; end;
if isempty(rgiant), rgiant=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(s3), s3=0; end;
x_shape=size(x);x=reshape(x,1,[]);
if isempty(x1max), x1max=0; end;
if isempty(x3max), x3max=0; end;
if isempty(xabs), xabs=0; end;
if isempty(zero), zero=0; end;
if firstCall,   one =[1.0d0];  end;
if firstCall,  zero =[0.0d0];  end;
if firstCall,  rdwarf =[3.834d-20];  end;
if firstCall,  rgiant=[1.304d19];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DENORM
s1 = zero;
s2 = zero;
s3 = zero;
x1max = zero;
x3max = zero;
floatn = n;
agiant = rgiant./floatn;
for i = 1 : n;
xabs = abs(x(i));
if( xabs>rdwarf && xabs<agiant )
%
%           SUM FOR INTERMEDIATE COMPONENTS.
%
s2 = s2 + xabs.^2;
elseif( xabs<=rdwarf ) ;
%
%              SUM FOR SMALL COMPONENTS.
%
if( xabs<=x3max )
if( xabs~=zero )
s3 = s3 +(xabs./x3max).^2;
end;
else;
s3 = one + s3.*(x3max./xabs).^2;
x3max = xabs;
end;
%
%              SUM FOR LARGE COMPONENTS.
%
elseif( xabs<=x1max ) ;
s1 = s1 +(xabs./x1max).^2;
else;
s1 = one + s1.*(x1max./xabs).^2;
x1max = xabs;
end;
end; i = fix(n+1);
%
%     CALCULATION OF NORM.
%
if( s1~=zero )
denormresult = x1max.*sqrt(s1+(s2./x1max)./x1max);
elseif( s2==zero ) ;
denormresult = x3max.*sqrt(s3);
else;
if( s2>=x3max )
denormresult = sqrt(s2.*(one+(x3max./s2).*(x3max.*s3)));
end;
if( s2<x3max )
denormresult = sqrt(x3max.*((s2./x3max)+(x3max.*s3)));
end;
end;
%
%     LAST CARD OF FUNCTION DENORM.
%
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DERFC