Return to snrm2.f CVS log

Up to [Research Unix] / researchv10no / cmd / view2d / num

Annotation of researchv10no/cmd/view2d/num/snrm2.f, revision 1.1

1.1     ! root        1:       real function snrm2 ( n, sx, incx)
        !             2:       integer          next
        !             3:       real   sx(1),  cutlo, cuthi, hitest, sum, xmax, zero, one
        !             4:       data   zero, one /0.0e0, 1.0e0/
        !             5: c
        !             6: c     euclidean norm of the n-vector stored in sx() with storage
        !             7: c     increment incx .
        !             8: c     if    n .le. 0 return with result = 0.
        !             9: c     if n .ge. 1 then incx must be .ge. 1
        !            10: c
        !            11: c           c.l.lawson, 1978 jan 08
        !            12: c
        !            13: c     four phase method     using two built-in constants that are
        !            14: c     hopefully applicable to all machines.
        !            15: c         cutlo = maximum of  sqrt(u/eps)  over all known machines.
        !            16: c         cuthi = minimum of  sqrt(v)      over all known machines.
        !            17: c     where
        !            18: c         eps = smallest no. such that eps + 1. .gt. 1.
        !            19: c         u   = smallest positive no.   (underflow limit)
        !            20: c         v   = largest  no.            (overflow  limit)
        !            21: c
        !            22: c     brief outline of algorithm..
        !            23: c
        !            24: c     phase 1    scans zero components.
        !            25: c     move to phase 2 when a component is nonzero and .le. cutlo
        !            26: c     move to phase 3 when a component is .gt. cutlo
        !            27: c     move to phase 4 when a component is .ge. cuthi/m
        !            28: c     where m = n for x() real and m = 2*n for complex.
        !            29: c
        !            30: c     values for cutlo and cuthi..
        !            31: c     from the environmental parameters listed in the imsl converter
        !            32: c     document the limiting values are as follows..
        !            33: c     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
        !            34: c                   univac and dec at 2**(-103)
        !            35: c                   thus cutlo = 2**(-51) = 4.44089e-16
        !            36: c     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
        !            37: c                   thus cuthi = 2**(63.5) = 1.30438e19
        !            38: c     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
        !            39: c                   thus cutlo = 2**(-33.5) = 8.23181d-11
        !            40: c     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
        !            41: c     data cutlo, cuthi / 8.232d-11,  1.304d19 /
        !            42: c     data cutlo, cuthi / 4.441e-16,  1.304e19 /
        !            43:       data cutlo, cuthi / 4.441e-16,  1.304e19 /
        !            44: c
        !            45:       if(n .gt. 0) go to 10
        !            46:          snrm2  = zero
        !            47:          go to 300
        !            48: c
        !            49:    10 assign 30 to next
        !            50:       sum = zero
        !            51:       nn = n * incx
        !            52: c                                                 begin main loop
        !            53:       i = 1
        !            54:    20    go to next,(30, 50, 70, 110)
        !            55:    30 if( abs(sx(i)) .gt. cutlo) go to 85
        !            56:       assign 50 to next
        !            57:       xmax = zero
        !            58: c
        !            59: c                        phase 1.  sum is zero
        !            60: c
        !            61:    50 if( sx(i) .eq. zero) go to 200
        !            62:       if( abs(sx(i)) .gt. cutlo) go to 85
        !            63: c
        !            64: c                                prepare for phase 2.
        !            65:       assign 70 to next
        !            66:       go to 105
        !            67: c
        !            68: c                                prepare for phase 4.
        !            69: c
        !            70:   100 i = j
        !            71:       assign 110 to next
        !            72:       sum = (sum / sx(i)) / sx(i)
        !            73:   105 xmax = abs(sx(i))
        !            74:       go to 115
        !            75: c
        !            76: c                   phase 2.  sum is small.
        !            77: c                             scale to avoid destructive underflow.
        !            78: c
        !            79:    70 if( abs(sx(i)) .gt. cutlo ) go to 75
        !            80: c
        !            81: c                     common code for phases 2 and 4.
        !            82: c                     in phase 4 sum is large.  scale to avoid overflow.
        !            83: c
        !            84:   110 if( abs(sx(i)) .le. xmax ) go to 115
        !            85:          sum = one + sum * (xmax / sx(i))**2
        !            86:          xmax = abs(sx(i))
        !            87:          go to 200
        !            88: c
        !            89:   115 sum = sum + (sx(i)/xmax)**2
        !            90:       go to 200
        !            91: c
        !            92: c
        !            93: c                  prepare for phase 3.
        !            94: c
        !            95:    75 sum = (sum * xmax) * xmax
        !            96: c
        !            97: c
        !            98: c     for real or d.p. set hitest = cuthi/n
        !            99: c     for complex      set hitest = cuthi/(2*n)
        !           100: c
        !           101:    85 hitest = cuthi/float( n )
        !           102: c
        !           103: c                   phase 3.  sum is mid-range.  no scaling.
        !           104: c
        !           105:       do 95 j =i,nn,incx
        !           106:       if(abs(sx(j)) .ge. hitest) go to 100
        !           107:    95    sum = sum + sx(j)**2
        !           108:       snrm2 = sqrt( sum )
        !           109:       go to 300
        !           110: c
        !           111:   200 continue
        !           112:       i = i + incx
        !           113:       if ( i .le. nn ) go to 20
        !           114: c
        !           115: c              end of main loop.
        !           116: c
        !           117: c              compute square root and adjust for scaling.
        !           118: c
        !           119:       snrm2 = xmax * sqrt(sum)
        !           120:   300 continue
        !           121:       return
        !           122:       end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.