| [scnrm2result,n,cx,incx]=scnrm2(n,cx,incx); |
function [scnrm2result,n,cx,incx]=scnrm2(n,cx,incx);
scnrm2result=[];
persistent absx cuthi cutlo firstCall gt hitest i igo imagmlv next nn one scalemlv summlv xmax zero ; if isempty(firstCall),firstCall=1;end;
;
if isempty(i), i=0; end;
if isempty(nn), nn=0; end;
%***BEGIN PROLOGUE SCNRM2
%***PURPOSE Compute the unitary norm of a complex vector.
%***LIBRARY SLATEC (BLAS)
%***CATEGORY D1A3B
%***TYPE COMPLEX (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)
% CX complex vector with N elements
% INCX storage spacing between elements of CX
%
% --Output--
% SCNRM2 single precision result (zero if N .LE. 0)
%
% Unitary norm of the complex N-vector stored in CX 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 SCNRM2
if isempty(imagmlv), imagmlv=false; end;
if isempty(scalemlv), scalemlv=false; end;
if isempty(next), next=0; end;
if isempty(igo), igo=0; end;
if isempty(gt), gt=zeros(1,40); end;
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(absx), absx=0; end;
if isempty(zero), zero=0; end;
if isempty(one), one=0; end;
cx_shape=size(cx);cx=reshape(cx,1,[]);
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 SCNRM2
if( n>0 )
%
next = 20;
summlv = zero;
nn = fix(n.*incx);
else;
scnrm2result = zero;
cx_shape=zeros(cx_shape);cx_shape(:)=cx(1:numel(cx_shape));cx=cx_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
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
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',cx); evalin('caller',[inputname(2),'=FUntemp;']); end
return;
end;
%
% BEGIN MAIN LOOP
%
gt(:)=0;
for i = 1 : incx: nn ;
while (1);
if(gt(40)==0)
if(gt(35)==0)
if(gt(30)==0)
if(gt(25)==0)
if(gt(20)==0)
if(gt(15)==0)
if(gt(10)==0)
if(gt(5)==0)
absx = abs(real(cx(i)));
imagmlv = false;
if( next~=20 )
if( next==40 )
gt(10)=1;
continue;
end;
if( next==80 )
gt(25)=1;
continue;
end;
if( next==140 )
gt(40)=1;
continue;
end;
if( next==100 )
gt(35)=1;
continue;
end;
end;
end;
gt(5)=0;
if( absx>cutlo )
gt(30)=1;
continue;
end;
next = 40;
scalemlv = false;
%
% PHASE 1. SUM IS ZERO
%
end;
gt(10)=0;
if( absx==zero )
gt(20)=1;
continue;
end;
if( absx>cutlo )
gt(30)=1;
continue;
end;
%
% PREPARE FOR PHASE 2.
%
next = 80;
end;
gt(15)=0;
scalemlv = true;
xmax = absx;
%
summlv = summlv +(absx./xmax).^2;
%
% CONTROL SELECTION OF REAL AND IMAGINARY PARTS.
%
end;
gt(20)=0;
if( imagmlv )
break;
end;
%
absx = abs(imag(cx(i)));
imagmlv = true;
if( next==20 )
gt(5)=1;
continue;
end;
if( next==40 )
gt(10)=1;
continue;
end;
if( next~=80 )
if( next==140 )
gt(40)=1;
continue;
end;
if( next~=100 )
break;
end;
gt(35 )=1;
continue;
end;
%
% PHASE 2. SUM IS SMALL.
% SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
%
end;
gt(25)=0;
if( absx<=cutlo )
gt(35)=1;
continue;
end;
%
% PREPARE FOR PHASE 3.
%
summlv =(summlv.*xmax).*xmax;
%
end;
gt(30)=0;
next = 140;
scalemlv = false;
%
% FOR REAL OR D.P. SET HITEST = CUTHI/N
% FOR COMPLEX SET HITEST = CUTHI/(2*N)
%
hitest = cuthi./n;
gt(40 )=1;
continue;
%
% COMMON CODE FOR PHASES 2 AND 4.
% IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
%
end;
gt(35)=0;
if( absx<=xmax )
summlv = summlv +(absx./xmax).^2;
else;
summlv = one + summlv.*(xmax./absx).^2;
xmax = absx;
end;
gt(20 )=1;
continue;
%
% PHASE 3. SUM IS MID-RANGE. NO SCALING.
%
end;
gt(40)=0;
if( absx<hitest )
summlv = summlv + absx.^2;
gt(20 )=1;
continue;
end;
%
% PREPARE FOR PHASE 4.
%
next = 100;
summlv =(summlv./absx)./absx;
gt(15 )=1;
continue;
break;
end;
end;
%
% end OF MAIN LOOP.
% COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
%
scnrm2result = sqrt(summlv);
if( scalemlv )
scnrm2result = scnrm2result.*xmax;
end;
cx_shape=zeros(cx_shape);cx_shape(:)=cx(1:numel(cx_shape));cx=cx_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
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
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',cx); evalin('caller',[inputname(2),'=FUntemp;']); end
end %function scnrm2
%!!!DECK SCNRM2
%!!REAL function SCNRM2(N,Cx,Incx)
%!!IMPLICIT NONE
%!!!*--SCNRM25
%!!INTEGER i , Incx , N , nn
%!!!***BEGIN PROLOGUE SCNRM2
%!!!***PURPOSE Compute the unitary norm of a complex vector.
%!!!***LIBRARY SLATEC (BLAS)
%!!!***CATEGORY D1A3B
%!!!***TYPE COMPLEX (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)
%!!! CX complex vector with N elements
%!!! INCX storage spacing between elements of CX
%!!!
%!!! --Output--
%!!! SCNRM2 single precision result (zero if N .LE. 0)
%!!!
%!!! Unitary norm of the complex N-vector stored in CX 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 SCNRM2
%!!LOGICAL imagmlv , scale
%!!INTEGER next , igo,gt(40)
%!!REAL cutlo , cuthi , hitest , sum , xmax , absx , zero , one
%!!COMPLEX Cx(*)
%!!SAVE cutlo , cuthi , zero , one
%!!DATA zero , one/0.0E0 , 1.0E0/
%!!!
%!!DATA cutlo , cuthi/4.441E-16 , 1.304E19/
%!!!***FIRST EXECUTABLE STATEMENT SCNRM2
%!!IF ( N>0 ) THEN
%!! !
%!! next = 20
%!! sum = zero
%!! nn = N*Incx
%!!ELSE
%!! SCNRM2 = zero
%!! RETURN
%!!ENDIF
%!!!
%!!! BEGIN MAIN LOOP
%!!!
%!!gt=0
%!!DO i = 1 , nn , Incx
%!! do
%!! if (gt(40)==0) then
%!! if (gt(35)==0) then
%!! if (gt(30)==0) then
%!! if (gt(25)==0) then
%!! if (gt(20)==0) then
%!! if (gt(15)==0) then
%!! if (gt(10)==0) then
%!! if (gt(5)==0) then
%!! absx = ABS(REAL(Cx(i)))
%!! imagmlv = false
%!! IF ( next~=20 ) THEN
%!! IF ( next==40 ) GOTO 100
%!! IF ( next==80 ) GOTO 250
%!! IF ( next==140 ) GOTO 400
%!! IF ( next==100 ) GOTO 350
%!! ENDIF
%!! endif
%!! gt(5)=0
%!!50 IF ( absx>cutlo ) GOTO 300
%!! next = 40
%!! scale = false
%!! !
%!! ! PHASE 1. SUM IS ZERO
%!! !
%!! endif
%!! gt(10)=0
%!!100 IF ( absx==zero ) GOTO 200
%!! IF ( absx>cutlo ) GOTO 300
%!! !
%!! ! PREPARE FOR PHASE 2.
%!! !
%!! next = 80
%!! endif
%!! gt(15)=0
%!!150 scale = true
%!! xmax = absx
%!! !
%!! sum = sum + (absx/xmax)**2
%!! !
%!! ! CONTROL SELECTION OF REAL AND IMAGINARY PARTS.
%!! !
%!! endif
%!! gt(20)=0
%!!200 IF ( imagmlv ) exit
%!! !
%!! absx = ABS(AIMAG(Cx(i)))
%!! imagmlv = true
%!! IF ( next==20 ) GOTO 50
%!! IF ( next==40 ) GOTO 100
%!! IF ( next~=80 ) THEN
%!! IF ( next==140 ) GOTO 400
%!! IF ( next~=100 ) exit
%!! gt(35 )=1
%!! cycle
%!! ENDIF
%!! !
%!! ! PHASE 2. SUM IS SMALL.
%!! ! SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
%!! !
%!! endif
%!! gt(25)=0
%!!250 IF ( absx<=cutlo ) GOTO 350
%!! !
%!! ! PREPARE FOR PHASE 3.
%!! !
%!! sum = (sum*xmax)*xmax
%!! !
%!! endif
%!! gt(30)=0
%!!300 next = 140
%!! scale = false
%!! !
%!! ! FOR REAL OR D.P. SET HITEST = CUTHI/N
%!! ! FOR COMPLEX SET HITEST = CUTHI/(2*N)
%!! !
%!! hitest = cuthi/N
%!! GOTO 400
%!! !
%!! ! COMMON CODE FOR PHASES 2 AND 4.
%!! ! IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
%!! !
%!! endif
%!! gt(35)=0
%!!350 IF ( absx<=xmax ) THEN
%!! sum = sum + (absx/xmax)**2
%!! ELSE
%!! sum = one + sum*(xmax/absx)**2
%!! xmax = absx
%!! ENDIF
%!! GOTO 200
%!! !
%!! ! PHASE 3. SUM IS MID-RANGE. NO SCALING.
%!! !
%!! endif
%!! gt(40)=0
%!!400 IF ( absx<hitest ) THEN
%!! sum = sum + absx**2
%!! GOTO 200
%!! ENDIF
%!! !
%!! ! PREPARE FOR PHASE 4.
%!! !
%!! next = 100
%!! sum = (sum/absx)/absx
%!! GOTO 150
%!! exit
%!! enddo
%!!ENDDO
%!!!
%!!! end OF MAIN LOOP.
%!!! COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
%!!!
%!!SCNRM2 = SQRT(sum)
%!!IF ( scale ) SCNRM2 = SCNRM2*xmax
%!!end function SCNRM2
%DECK SCOEF
|
|
