Mercurial > octave-nkf
view libcruft/odessa/odessa_vnorm.f @ 4720:e759d01692db ss-2-1-53
[project @ 2004-01-23 04:13:37 by jwe]
author | jwe |
---|---|
date | Fri, 23 Jan 2004 04:13:37 +0000 |
parents | 258c1d15ad78 |
children |
line wrap: on
line source
DOUBLE PRECISION FUNCTION ODESSA_VNORM (N, V, W) C----------------------------------------------------------------------- C THIS FUNCTION ROUTINE COMPUTES THE WEIGHTED ROOT-MEAN-SQUARE NORM C OF THE VECTOR OF LENGTH N CONTAINED IN THE ARRAY V, WITH WEIGHTS C CONTAINED IN THE ARRAY W OF LENGTH N.. C ODESSA_VNORM = SQRT( (1/N) * SUM( V(I)*W(I) )**2 ) C PROTECTION FOR UNDERFLOW/OVERFLOW IS ACCOMPLISHED USING TWO C CONSTANTS WHICH ARE HOPEFULLY APPLICABLE FOR ALL MACHINES. C THESE ARE: C CUTLO = maximum of SQRT(U/EPS) over all known machines C CUTHI = minimum of SQRT(Z) over all known machines C WHERE C EPS = smallest number s.t. EPS + 1 .GT. 1 C U = smallest positive number (underflow limit) C Z = largest number (overflow limit) C C DETAILS OF THE ALGORITHM AND OF VALUES OF CUTLO AND CUTHI ARE C FOUND IN THE BLAS ROUTINE SNRM2 (SEE ALSO ALGORITHM 539, TRANS. C MATH. SOFTWARE, VOL. 5 NO. 3, 1979, 308-323. C FOR SINGLE PRECISION, THE FOLLOWING VALUES SHOULD BE UNIVERSAL: C DATA CUTLO,CUTHI /4.441E-16,1.304E19/ C FOR DOUBLE PRECISION, USE C DATA CUTLO,CUTHI /8.232D-11,1.304D19/ C C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEXT,I,J,N DIMENSION V(N),W(N) DATA CUTLO,CUTHI /8.232D-11,1.304D19/ DATA ZERO,ONE/0.0D0,1.0D0/ C BLAS ALGORITHM NEXT = 1 SUM = ZERO I = 1 20 SX = V(I)*W(I) GO TO (30,40,70,80),NEXT 30 IF (DABS(SX).GT.CUTLO) GO TO 110 NEXT = 2 XMAX = ZERO 40 IF (SX.EQ.ZERO) GO TO 130 IF (DABS(SX).GT.CUTLO) GO TO 110 NEXT = 3 GO TO 60 50 I=J NEXT = 4 SUM = (SUM/SX)/SX 60 XMAX = DABS(SX) GO TO 90 70 IF(DABS(SX).GT.CUTLO) GO TO 100 80 IF(DABS(SX).LE.XMAX) GO TO 90 SUM = ONE + SUM * (XMAX/SX)**2 XMAX = DABS(SX) GO TO 130 90 SUM = SUM + (SX/XMAX)**2 GO TO 130 100 SUM = (SUM*XMAX)*XMAX 110 HITEST = CUTHI/DBLE(N) DO 120 J = I,N SX = V(J)*W(J) IF(DABS(SX).GE.HITEST) GO TO 50 SUM = SUM + SX**2 120 CONTINUE ODESSA_VNORM = DSQRT(SUM) GO TO 140 130 CONTINUE I = I + 1 IF (I.LE.N) GO TO 20 ODESSA_VNORM = XMAX * DSQRT(SUM) 140 CONTINUE RETURN C----------------------- END OF FUNCTION ODESSA_VNORM ------------------------- END