gcc/testsuite/gfortran.dg/norm2_3.f90

   1 ! { dg-do run }
   2 ! { dg-require-effective-target fortran_large_real }
   3 !
   4 !
   5 ! PR fortran/33197
   6 !
   7 ! Check implementation of L2 norm (Euclidean vector norm)
   8 !
   9 implicit none
  10
  11 integer,parameter :: qp = selected_real_kind (precision (0.0d0)+1)
  12
  13 real(qp) :: a(3) = [real(qp) :: 1, 2, huge(3.0_qp)]
  14 real(qp) :: b(3) = [real(qp) :: 1, 2, 3]
  15 real(qp) :: c(4) = [real(qp) :: 1, 2, 3, -1]
  16 real(qp) :: e(0) = [real(qp) :: ]
  17 real(qp) :: f(4) = [real(qp) :: 0, 0, 3, 0 ]
  18
  19 real(qp) :: d(4,1) = RESHAPE ([real(qp) :: 1, 2, 3, -1], [4,1])
  20 real(qp) :: g(4,1) = RESHAPE ([real(qp) :: 0, 0, 4, -1], [4,1])
  21
  22 ! Check compile-time version
  23
  24 if (abs (NORM2 ([real(qp) :: 1, 2, huge(3.0_qp)])   - huge(3.0_qp)) &
  25     > epsilon(0.0_qp)*huge(3.0_qp)) call abort()
  26
  27 if (abs (SNORM2([real(qp) :: 1, 2, huge(3.0_qp)],3) - huge(3.0_qp)) &
  28     > epsilon(0.0_qp)*huge(3.0_qp)) call abort()
  29
  30 if (abs (SNORM2([real(qp) :: 1, 2, 3],3) - NORM2([real(qp) :: 1, 2, 3])) &
  31     > epsilon(0.0_qp)*SNORM2([real(qp) :: 1, 2, 3],3)) call abort()
  32
  33 if (NORM2([real(qp) :: ]) /= 0.0_qp) call abort()
  34 if (abs (NORM2([real(qp) :: 0, 0, 3, 0]) - 3.0_qp) > epsilon(0.0_qp)) call abort()
  35
  36 ! Check TREE version
  37
  38 if (abs (NORM2 (a)   - huge(3.0_qp)) &
  39     > epsilon(0.0_qp)*huge(3.0_qp)) call abort()
  40
  41 if (abs (SNORM2(b,3) - NORM2(b)) &
  42     > epsilon(0.0_qp)*SNORM2(b,3)) call abort()
  43
  44 if (abs (SNORM2(c,4) - NORM2(c)) &
  45     > epsilon(0.0_qp)*SNORM2(c,4)) call abort()
  46
  47 if (ANY (abs (abs(d(:,1)) - NORM2(d, 2)) &
  48     > epsilon(0.0_qp))) call abort()
  49
  50 ! Check libgfortran version
  51
  52 if (ANY (abs (SNORM2(d,4) - NORM2(d, 1)) &
  53     > epsilon(0.0_qp)*SNORM2(d,4))) call abort()
  54
  55 if (abs (SNORM2(f,4) - NORM2(f, 1)) &
  56     > epsilon(0.0_qp)*SNORM2(d,4)) call abort()
  57
  58 if (ANY (abs (abs(g(:,1)) - NORM2(g, 2)) &
  59     > epsilon(0.0_qp))) call abort()
  60
  61 contains
  62    ! NORM2 algorithm based on BLAS, cf.
  63    ! http://www.netlib.org/blas/snrm2.f
  64    REAL(qp) FUNCTION SNORM2 (X,n)
  65       INTEGER, INTENT(IN) :: n
  66       REAL(qp), INTENT(IN) :: X(n)
  67
  68       REAL(qp) :: absXi, scale, SSQ
  69       INTEGER :: i
  70
  71       INTRINSIC :: ABS, SQRT
  72
  73       IF (N < 1) THEN
  74         snorm2 = 0.0_qp
  75       ELSE IF (N == 1) THEN
  76         snorm2 = ABS(X(1))
  77       ELSE
  78           scale = 0.0_qp
  79           SSQ = 1.0_qp
  80
  81           DO i = 1, N
  82               IF (X(i) /= 0.0_qp) THEN
  83                   absXi = ABS(X(i))
  84                   IF (scale < absXi) THEN
  85                       SSQ = 1.0_qp + SSQ * (scale/absXi)**2
  86                       scale = absXi
  87                   ELSE
  88                       SSQ = SSQ + (absXi/scale)**2
  89                   END IF
  90               END IF
  91           END DO
  92           snorm2 = scale * SQRT(SSQ)
  93       END IF
  94    END FUNCTION SNORM2
  95 end