LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zb1nrm2()

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

Definition at line 790 of file zblat1.f.

791* Compare NRM2 with a reference computation using combinations
792* of the following values:
793*
794* 0, very small, small, ulp, 1, 1/ulp, big, very big, infinity, NaN
795*
796* one of these values is used to initialize x(1) and x(2:N) is
797* filled with random values from [-1,1] scaled by another of
798* these values.
799*
800* This routine is adapted from the test suite provided by
801* Anderson E. (2017)
802* Algorithm 978: Safe Scaling in the Level 1 BLAS
803* ACM Trans Math Softw 44:1--28
804* https://doi.org/10.1145/3061665
805*
806* .. Scalar Arguments ..
807 INTEGER INCX, N
808 DOUBLE PRECISION THRESH
809*
810* =====================================================================
811* .. Parameters ..
812 INTEGER NMAX, NOUT, NV
813 parameter(nmax=20, nout=6, nv=10)
814 DOUBLE PRECISION HALF, ONE, THREE, TWO, ZERO
815 parameter(half=0.5d+0, one=1.0d+0, two= 2.0d+0,
816 & three=3.0d+0, zero=0.0d+0)
817* .. External Functions ..
818 DOUBLE PRECISION DZNRM2
819 EXTERNAL dznrm2
820* .. Intrinsic Functions ..
821 INTRINSIC aimag, abs, dcmplx, dble, max, min, sqrt
822* .. Model parameters ..
823 DOUBLE PRECISION BIGNUM, SAFMAX, SAFMIN, SMLNUM, ULP
824 parameter(bignum=0.99792015476735990583d+292,
825 & safmax=0.44942328371557897693d+308,
826 & safmin=0.22250738585072013831d-307,
827 & smlnum=0.10020841800044863890d-291,
828 & ulp=0.22204460492503130808d-015)
829* .. Local Scalars ..
830 COMPLEX*16 ROGUE
831 DOUBLE PRECISION SNRM, TRAT, V0, V1, WORKSSQ, Y1, Y2,
832 & YMAX, YMIN, YNRM, ZNRM
833 INTEGER I, IV, IW, IX, KS
834 LOGICAL FIRST
835* .. Local Arrays ..
836 COMPLEX*16 X(NMAX), Z(NMAX)
837 DOUBLE PRECISION VALUES(NV), WORK(NMAX)
838* .. Executable Statements ..
839 values(1) = zero
840 values(2) = two*safmin
841 values(3) = smlnum
842 values(4) = ulp
843 values(5) = one
844 values(6) = one / ulp
845 values(7) = bignum
846 values(8) = safmax
847 values(9) = dxvals(v0,2)
848 values(10) = dxvals(v0,3)
849 rogue = dcmplx(1234.5678d+0,-1234.5678d+0)
850 first = .true.
851*
852* Check that the arrays are large enough
853*
854 IF (n*abs(incx).GT.nmax) THEN
855 WRITE (nout,99) "DZNRM2", nmax, incx, n, n*abs(incx)
856 RETURN
857 END IF
858*
859* Zero-sized inputs are tested in STEST1.
860 IF (n.LE.0) THEN
861 RETURN
862 END IF
863*
864* Generate 2*(N-1) values in (-1,1).
865*
866 ks = 2*(n-1)
867 DO i = 1, ks
868 CALL random_number(work(i))
869 work(i) = one - two*work(i)
870 END DO
871*
872* Compute the sum of squares of the random values
873* by an unscaled algorithm.
874*
875 workssq = zero
876 DO i = 1, ks
877 workssq = workssq + work(i)*work(i)
878 END DO
879*
880* Construct the test vector with one known value
881* and the rest from the random work array multiplied
882* by a scaling factor.
883*
884 DO iv = 1, nv
885 v0 = values(iv)
886 IF (abs(v0).GT.one) THEN
887 v0 = v0*half*half
888 END IF
889 z(1) = dcmplx(v0,-three*v0)
890 DO iw = 1, nv
891 v1 = values(iw)
892 IF (abs(v1).GT.one) THEN
893 v1 = (v1*half) / sqrt(dble(ks+1))
894 END IF
895 DO i = 1, n-1
896 z(i+1) = dcmplx(v1*work(2*i-1),v1*work(2*i))
897 END DO
898*
899* Compute the expected value of the 2-norm
900*
901 y1 = abs(v0) * sqrt(10.0d0)
902 IF (n.GT.1) THEN
903 y2 = abs(v1)*sqrt(workssq)
904 ELSE
905 y2 = zero
906 END IF
907 ymin = min(y1, y2)
908 ymax = max(y1, y2)
909*
910* Expected value is NaN if either is NaN. The test
911* for YMIN == YMAX avoids further computation if both
912* are infinity.
913*
914 IF ((y1.NE.y1).OR.(y2.NE.y2)) THEN
915* add to propagate NaN
916 ynrm = y1 + y2
917 ELSE IF (ymin == ymax) THEN
918 ynrm = sqrt(two)*ymax
919 ELSE IF (ymax == zero) THEN
920 ynrm = zero
921 ELSE
922 ynrm = ymax*sqrt(one + (ymin / ymax)**2)
923 END IF
924*
925* Fill the input array to DZNRM2 with steps of incx
926*
927 DO i = 1, n
928 x(i) = rogue
929 END DO
930 ix = 1
931 IF (incx.LT.0) ix = 1 - (n-1)*incx
932 DO i = 1, n
933 x(ix) = z(i)
934 ix = ix + incx
935 END DO
936*
937* Call DZNRM2 to compute the 2-norm
938*
939 snrm = dznrm2(n,x,incx)
940*
941* Compare SNRM and ZNRM. Roundoff error grows like O(n)
942* in this implementation so we scale the test ratio accordingly.
943*
944 IF (incx.EQ.0) THEN
945 y1 = abs(dble(x(1)))
946 y2 = abs(aimag(x(1)))
947 ymin = min(y1, y2)
948 ymax = max(y1, y2)
949 IF ((y1.NE.y1).OR.(y2.NE.y2)) THEN
950* add to propagate NaN
951 znrm = y1 + y2
952 ELSE IF (ymin == ymax) THEN
953 znrm = sqrt(two)*ymax
954 ELSE IF (ymax == zero) THEN
955 znrm = zero
956 ELSE
957 znrm = ymax * sqrt(one + (ymin / ymax)**2)
958 END IF
959 znrm = sqrt(dble(n)) * znrm
960 ELSE
961 znrm = ynrm
962 END IF
963*
964* The tests for NaN rely on the compiler not being overly
965* aggressive and removing the statements altogether.
966 IF ((snrm.NE.snrm).OR.(znrm.NE.znrm)) THEN
967 IF ((snrm.NE.snrm).NEQV.(znrm.NE.znrm)) THEN
968 trat = one / ulp
969 ELSE
970 trat = zero
971 END IF
972 ELSE IF (znrm == zero) THEN
973 trat = snrm / ulp
974 ELSE
975 trat = (abs(snrm-znrm) / znrm) / (two*dble(n)*ulp)
976 END IF
977 IF ((trat.NE.trat).OR.(trat.GE.thresh)) THEN
978 IF (first) THEN
979 first = .false.
980 WRITE(nout,99999)
981 END IF
982 WRITE (nout,98) "DZNRM2", n, incx, iv, iw, trat
983 END IF
984 END DO
985 END DO
98699999 FORMAT (' FAIL')
987 99 FORMAT ( ' Not enough space to test ', a6, ': NMAX = ',i6,
988 + ', INCX = ',i6,/,' N = ',i6,', must be at least ',i6 )
989 98 FORMAT( 1x, a6, ': N=', i6,', INCX=', i4, ', IV=', i2, ', IW=',
990 + i2, ', test=', e15.8 )
991 RETURN
992 CONTAINS
993 DOUBLE PRECISION FUNCTION dxvals(XX,K)
994* .. Scalar Arguments ..
995 DOUBLE PRECISION XX
996 INTEGER K
997* .. Local Scalars ..
998 DOUBLE PRECISION X, Y, YY, Z
999* .. Intrinsic Functions ..
1000 INTRINSIC huge
1001* .. Executable Statements ..
1002 y = huge(xx)
1003 z = yy
1004 IF (k.EQ.1) THEN
1005 x = -z
1006 ELSE IF (k.EQ.2) THEN
1007 x = z
1008 ELSE IF (k.EQ.3) THEN
1009 x = z / z
1010 END IF
1011 dxvals = x
1012 RETURN
1013 END
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
Here is the caller graph for this function: