Parent Directory | Revision Log

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

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

File size: 3986 byte(s)

File size: 3986 byte(s)

Setting up web subdirectory in repository

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 |