Parent Directory | Revision Log

Revision **1** -
(**hide annotations**)
(**download**)

*Fri Oct 29 20:54:12 2004 UTC*
(19 years, 5 months ago)
by *aw0a*

File size: 3986 byte(s)

File size: 3986 byte(s)

Setting up web subdirectory in repository

1 | aw0a | 1 | C dnrm2.f |

2 | C is freely available from netlib. It is not subject to any GNU License | ||

3 | C set by the authors of the ASCEND math programming system. | ||

4 | C $Date: 1996/04/30 18:14:07 $ $Revision: 1.1.1.1 $ | ||

5 | C | ||

6 | double precision function dnrm2 ( n, dx, incx) | ||

7 | integer next | ||

8 | double precision dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one | ||

9 | data zero, one /0.0d0, 1.0d0/ | ||

10 | c | ||

11 | c euclidean norm of the n-vector stored in dx() with storage | ||

12 | c increment incx . | ||

13 | c if n .le. 0 return with result = 0. | ||

14 | c if n .ge. 1 then incx must be .ge. 1 | ||

15 | c | ||

16 | c c.l.lawson, 1978 jan 08 | ||

17 | c | ||

18 | c four phase method using two built-in constants that are | ||

19 | c hopefully applicable to all machines. | ||

20 | c cutlo = maximum of dsqrt(u/eps) over all known machines. | ||

21 | c cuthi = minimum of dsqrt(v) over all known machines. | ||

22 | c where | ||

23 | c eps = smallest no. such that eps + 1. .gt. 1. | ||

24 | c u = smallest positive no. (underflow limit) | ||

25 | c v = largest no. (overflow limit) | ||

26 | c | ||

27 | c brief outline of algorithm.. | ||

28 | c | ||

29 | c phase 1 scans zero components. | ||

30 | c move to phase 2 when a component is nonzero and .le. cutlo | ||

31 | c move to phase 3 when a component is .gt. cutlo | ||

32 | c move to phase 4 when a component is .ge. cuthi/m | ||

33 | c where m = n for x() real and m = 2*n for complex. | ||

34 | c | ||

35 | c values for cutlo and cuthi.. | ||

36 | c from the environmental parameters listed in the imsl converter | ||

37 | c document the limiting values are as follows.. | ||

38 | c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are | ||

39 | c univac and dec at 2**(-103) | ||

40 | c thus cutlo = 2**(-51) = 4.44089e-16 | ||

41 | c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. | ||

42 | c thus cuthi = 2**(63.5) = 1.30438e19 | ||

43 | c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. | ||

44 | c thus cutlo = 2**(-33.5) = 8.23181d-11 | ||

45 | c cuthi, d.p. same as s.p. cuthi = 1.30438d19 | ||

46 | c data cutlo, cuthi / 8.232d-11, 1.304d19 / | ||

47 | c data cutlo, cuthi / 4.441e-16, 1.304e19 / | ||

48 | data cutlo, cuthi / 8.232d-11, 1.304d19 / | ||

49 | c | ||

50 | if(n .gt. 0) go to 10 | ||

51 | dnrm2 = zero | ||

52 | go to 300 | ||

53 | c | ||

54 | 10 assign 30 to next | ||

55 | sum = zero | ||

56 | nn = n * incx | ||

57 | c begin main loop | ||

58 | i = 1 | ||

59 | 20 go to next,(30, 50, 70, 110) | ||

60 | 30 if( dabs(dx(i)) .gt. cutlo) go to 85 | ||

61 | assign 50 to next | ||

62 | xmax = zero | ||

63 | c | ||

64 | c phase 1. sum is zero | ||

65 | c | ||

66 | 50 if( dx(i) .eq. zero) go to 200 | ||

67 | if( dabs(dx(i)) .gt. cutlo) go to 85 | ||

68 | c | ||

69 | c prepare for phase 2. | ||

70 | assign 70 to next | ||

71 | go to 105 | ||

72 | c | ||

73 | c prepare for phase 4. | ||

74 | c | ||

75 | 100 i = j | ||

76 | assign 110 to next | ||

77 | sum = (sum / dx(i)) / dx(i) | ||

78 | 105 xmax = dabs(dx(i)) | ||

79 | go to 115 | ||

80 | c | ||

81 | c phase 2. sum is small. | ||

82 | c scale to avoid destructive underflow. | ||

83 | c | ||

84 | 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 | ||

85 | c | ||

86 | c common code for phases 2 and 4. | ||

87 | c in phase 4 sum is large. scale to avoid overflow. | ||

88 | c | ||

89 | 110 if( dabs(dx(i)) .le. xmax ) go to 115 | ||

90 | sum = one + sum * (xmax / dx(i))**2 | ||

91 | xmax = dabs(dx(i)) | ||

92 | go to 200 | ||

93 | c | ||

94 | 115 sum = sum + (dx(i)/xmax)**2 | ||

95 | go to 200 | ||

96 | c | ||

97 | c | ||

98 | c prepare for phase 3. | ||

99 | c | ||

100 | 75 sum = (sum * xmax) * xmax | ||

101 | c | ||

102 | c | ||

103 | c for real or d.p. set hitest = cuthi/n | ||

104 | c for complex set hitest = cuthi/(2*n) | ||

105 | c | ||

106 | 85 hitest = cuthi/float( n ) | ||

107 | c | ||

108 | c phase 3. sum is mid-range. no scaling. | ||

109 | c | ||

110 | do 95 j =i,nn,incx | ||

111 | if(dabs(dx(j)) .ge. hitest) go to 100 | ||

112 | 95 sum = sum + dx(j)**2 | ||

113 | dnrm2 = dsqrt( sum ) | ||

114 | go to 300 | ||

115 | c | ||

116 | 200 continue | ||

117 | i = i + incx | ||

118 | if ( i .le. nn ) go to 20 | ||

119 | c | ||

120 | c end of main loop. | ||

121 | c | ||

122 | c compute square root and adjust for scaling. | ||

123 | c | ||

124 | dnrm2 = xmax * dsqrt(sum) | ||

125 | 300 continue | ||

126 | return | ||

127 | end | ||

128 |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |