3047
3048
3049
3050
3051
3052
3053
3054 INTEGER INCX, N
3055 DOUBLE PRECISION ERRBND, PREC, USCLR
3056
3057
3058 DOUBLE PRECISION X( * )
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110 DOUBLE PRECISION ONE, TWO, ZERO
3111 parameter( one = 1.0d+0, two = 2.0d+0,
3112 $ zero = 0.0d+0 )
3113
3114
3115 INTEGER IX
3116 DOUBLE PRECISION ABSXI, ADDBND, FACT, SCALE, SSQ, SUMSCA, SUMSSQ
3117
3118
3119 INTRINSIC abs
3120
3121
3122
3123 usclr = zero
3124 sumssq = one
3125 sumsca = zero
3126 addbnd = two * two * two * prec
3127 fact = one + two * ( ( one + prec )**3 - one )
3128
3129 scale = zero
3130 ssq = one
3131 DO 10 ix = 1, 1 + ( n - 1 )*incx, incx
3132 IF( x( ix ).NE.zero ) THEN
3133 absxi = abs( x( ix ) )
3134 IF( scale.LT.absxi )THEN
3135 sumssq = one + ( ssq*( scale/absxi )**2 ) * fact
3136 errbnd = addbnd * sumssq
3137 sumssq = sumssq + errbnd
3138 ssq = one + ssq*( scale/absxi )**2
3139 sumsca = absxi
3140 scale = absxi
3141 ELSE
3142 sumssq = ssq + ( ( absxi/scale )**2 ) * fact
3143 errbnd = addbnd * sumssq
3144 sumssq = sumssq + errbnd
3145 ssq = ssq + ( absxi/scale )**2
3146 END IF
3147 END IF
3148 10 CONTINUE
3149
3150 usclr = scale * sqrt( ssq )
3151
3152
3153
3154 errbnd = sqrt( sumssq ) * ( one + two * ( 1.00001d+0 * prec ) )
3155
3156 errbnd = ( sumsca * errbnd ) - usclr
3157
3158 RETURN
3159
3160
3161