db1nrm2()

subroutine db1nrm2	(	integer	n,
		integer	incx,
		double precision	thresh )

Definition at line 1137 of file dblat1.f.

*     Compare NRM2 with a reference computation using combinations
*     of the following values:
*
*     0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN
*
*     one of these values is used to initialize x(1) and x(2:N) is
*     filled with random values from [-1,1] scaled by another of
*     these values.
*
*     This routine is adapted from the test suite provided by
*     Anderson E. (2017)
*     Algorithm 978: Safe Scaling in the Level 1 BLAS
*     ACM Trans Math Softw 44:1--28
*     https://doi.org/10.1145/3061665
*
*     .. Scalar Arguments ..
      INTEGER           INCX, N
      DOUBLE PRECISION  THRESH
*
*  =====================================================================
*     .. Parameters ..
      INTEGER           NMAX, NOUT, NV
      parameter(nmax=20, nout=6, nv=10)
      DOUBLE PRECISION  HALF, ONE, TWO, ZERO
      parameter(half=0.5d+0, one=1.0d+0, two= 2.0d+0,
     &                  zero=0.0d+0)
*     .. External Functions ..
      DOUBLE PRECISION  DNRM2
      EXTERNAL          dnrm2
*     .. Intrinsic Functions ..
      INTRINSIC         abs, dble, max, min, sqrt
*     .. Model parameters ..
      DOUBLE PRECISION  BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP
      parameter(bignum=0.99792015476735990583d+292,
     &                  safmax=0.44942328371557897693d+308,
     &                  safmin=0.22250738585072013831d-307,
     &                  smlnum=0.10020841800044863890d-291,
     &                  ulp=0.22204460492503130808d-015)
*     .. Local Scalars ..
      DOUBLE PRECISION  ROGUE, SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2,
     &                  YMAX, YMIN, YNRM, ZNRM
      INTEGER           I, IV, IW, IX
      LOGICAL           FIRST
*     .. Local Arrays ..
      DOUBLE PRECISION  VALUES(NV), WORK(NMAX), X(NMAX), Z(NMAX)
*     .. Executable Statements ..
      values(1) = zero
      values(2) = two*safmin
      values(3) = smlnum
      values(4) = ulp
      values(5) = one
      values(6) = one / ulp
      values(7) = bignum
      values(8) = safmax
      values(9) = dxvals(v0,2)
      values(10) = dxvals(v0,3)
      rogue = -1234.5678d+0
      first = .true.
*
*     Check that the arrays are large enough
*
      IF (n*abs(incx).GT.nmax) THEN
         WRITE (nout,99) "DNRM2", nmax, incx, n, n*abs(incx)
         RETURN
      END IF
*
*     Zero-sized inputs are tested in STEST1.
      IF (n.LE.0) THEN
         RETURN
      END IF
*
*     Generate (N-1) values in (-1,1).
*
      DO i = 2, n
         CALL random_number(work(i))
         work(i) = one - two*work(i)
      END DO
*
*     Compute the sum of squares of the random values
*     by an unscaled algorithm.
*
      workssq = zero
      DO i = 2, n
         workssq = workssq + work(i)*work(i)
      END DO
*
*     Construct the test vector with one known value
*     and the rest from the random work array multiplied
*     by a scaling factor.
*
      DO iv = 1, nv
         v0 = values(iv)
         IF (abs(v0).GT.one) THEN
            v0 = v0*half
         END IF
         z(1) = v0
         DO iw = 1, nv
            v1 = values(iw)
            IF (abs(v1).GT.one) THEN
               v1 = (v1*half) / sqrt(dble(n))
            END IF
            DO i = 2, n
               z(i) = v1*work(i)
            END DO
*
*           Compute the expected value of the 2-norm
*
            y1 = abs(v0)
            IF (n.GT.1) THEN
               y2 = abs(v1)*sqrt(workssq)
            ELSE
               y2 = zero
            END IF
            ymin = min(y1, y2)
            ymax = max(y1, y2)
*
*           Expected value is NaN if either is NaN. The test
*           for YMIN == YMAX avoids further computation if both
*           are infinity.
*
            IF ((y1.NE.y1).OR.(y2.NE.y2)) THEN
*              Add to propagate NaN
               ynrm = y1 + y2
            ELSE IF (ymax == zero) THEN
               ynrm = zero
            ELSE IF (ymin == ymax) THEN
               ynrm = sqrt(two)*ymax
            ELSE
               ynrm = ymax*sqrt(one + (ymin / ymax)**2)
            END IF
*
*           Fill the input array to DNRM2 with steps of incx
*
            DO i = 1, n
               x(i) = rogue
            END DO
            ix = 1
            IF (incx.LT.0) ix = 1 - (n-1)*incx
            DO i = 1, n
               x(ix) = z(i)
               ix = ix + incx
            END DO
*
*           Call DNRM2 to compute the 2-norm
*
            snrm = dnrm2(n,x,incx)
*
*           Compare SNRM and ZNRM.  Roundoff error grows like O(n)
*           in this implementation so we scale the test ratio accordingly.
*
            IF (incx.EQ.0) THEN
               znrm = sqrt(dble(n))*abs(x(1))
            ELSE
               znrm = ynrm
            END IF
*
*           The tests for NaN rely on the compiler not being overly
*           aggressive and removing the statements altogether.
            IF ((snrm.NE.snrm).OR.(znrm.NE.znrm)) THEN
               IF ((snrm.NE.snrm).NEQV.(znrm.NE.znrm)) THEN
                  trat = one / ulp
               ELSE
                  trat = zero
               END IF
            ELSE IF (snrm == znrm) THEN
               trat = zero
            ELSE IF (znrm == zero) THEN
               trat = snrm / ulp
            ELSE
               trat = (abs(snrm-znrm) / znrm) / (dble(n)*ulp)
            END IF
            IF ((trat.NE.trat).OR.(trat.GE.thresh)) THEN
               IF (first) THEN
                  first = .false.
                  WRITE(nout,99999)
               END IF
               WRITE (nout,98) "DNRM2", n, incx, iv, iw, trat
            END IF
         END DO
      END DO
99999 FORMAT ('                                       FAIL')
   99 FORMAT ( ' Not enough space to test ', a6, ': NMAX = ',i6,
     + ', INCX = ',i6,/,'   N = ',i6,', must be at least ',i6 )
   98 FORMAT( 1x, a6, ': N=', i6,', INCX=', i4, ', IV=', i2, ', IW=',
     +  i2, ', test=', e15.8 )
      RETURN
      CONTAINS
      DOUBLE PRECISION FUNCTION dxvals(XX,K)
*     .. Scalar Arguments ..
      DOUBLE PRECISION  XX
      INTEGER           K
*     .. Parameters ..
      DOUBLE PRECISION  ZERO
      parameter(zero=0.0d+0)
*     .. Local Scalars ..
      DOUBLE PRECISION  X, Y, Z
*     .. Intrinsic Functions ..
      INTRINSIC         huge
*     .. Executable Statements ..
      x = zero
      y = huge(xx)
      z = y*y
      IF (k.EQ.1) THEN
         x = -z
      ELSE IF (k.EQ.2) THEN
         x = z
      ELSE IF (k.EQ.3) THEN
         x = z / z
      END IF
      dxvals = x
      RETURN
      END

Here is the caller graph for this function: