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!
9implicit none
10
11integer,parameter :: qp = selected_real_kind (precision (0.0d0)+1)
12
13real(qp) :: a(3) = [real(qp) :: 1, 2, huge(3.0_qp)]
14real(qp) :: b(3) = [real(qp) :: 1, 2, 3]
15real(qp) :: c(4) = [real(qp) :: 1, 2, 3, -1]
16real(qp) :: e(0) = [real(qp) :: ]
17real(qp) :: f(4) = [real(qp) :: 0, 0, 3, 0 ]
18
19real(qp) :: d(4,1) = RESHAPE ([real(qp) :: 1, 2, 3, -1], [4,1])
20real(qp) :: g(4,1) = RESHAPE ([real(qp) :: 0, 0, 4, -1], [4,1])
21
22! Check compile-time version
23
24if (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
27if (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
30if (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
33if (NORM2([real(qp) :: ]) /= 0.0_qp) call abort()
34if (abs (NORM2([real(qp) :: 0, 0, 3, 0]) - 3.0_qp) > epsilon(0.0_qp)) call abort()
35
36! Check TREE version
37
38if (abs (NORM2 (a)   - huge(3.0_qp)) &
39    > epsilon(0.0_qp)*huge(3.0_qp)) call abort()
40
41if (abs (SNORM2(b,3) - NORM2(b)) &
42    > epsilon(0.0_qp)*SNORM2(b,3)) call abort()
43
44if (abs (SNORM2(c,4) - NORM2(c)) &
45    > epsilon(0.0_qp)*SNORM2(c,4)) call abort()
46
47if (ANY (abs (abs(d(:,1)) - NORM2(d, 2)) &
48    > epsilon(0.0_qp))) call abort()
49
50! Check libgfortran version
51
52if (ANY (abs (SNORM2(d,4) - NORM2(d, 1)) &
53    > epsilon(0.0_qp)*SNORM2(d,4))) call abort()
54
55if (abs (SNORM2(f,4) - NORM2(f, 1)) &
56    > epsilon(0.0_qp)*SNORM2(d,4)) call abort()
57
58if (ANY (abs (abs(g(:,1)) - NORM2(g, 2)) &
59    > epsilon(0.0_qp))) call abort()
60
61contains
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
95end
96