From b0b131e0870c0a8866decbcb95fb10dd455f307d Mon Sep 17 00:00:00 2001 From: Mark Gates Date: Fri, 2 Nov 2018 15:48:02 -0400 Subject: [PATCH 1/2] Update Frobenius norms for better accuracy. --- SRC/Makefile | 4 +-- SRC/clangb.f | 23 +++++++++---- SRC/clange.f | 22 ++++++++---- SRC/clanhb.f | 48 +++++++++++++++++++------- SRC/clanhe.f | 45 +++++++++++++++++------- SRC/clanhp.f | 46 ++++++++++++++++++------- SRC/clanhs.f | 23 +++++++++---- SRC/clansb.f | 42 +++++++++++++++++------ SRC/clansp.f | 54 ++++++++++++++++++++--------- SRC/clansy.f | 40 +++++++++++++++++----- SRC/clantb.f | 55 ++++++++++++++++++++---------- SRC/clantp.f | 53 +++++++++++++++++++++-------- SRC/clantr.f | 54 ++++++++++++++++++++--------- SRC/dcombssq.f | 92 ++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/dlangb.f | 26 +++++++++----- SRC/dlange.f | 22 ++++++++---- SRC/dlanhs.f | 25 ++++++++++---- SRC/dlansb.f | 44 ++++++++++++++++++------ SRC/dlansp.f | 48 +++++++++++++++++++------- SRC/dlansy.f | 42 +++++++++++++++++------ SRC/dlantb.f | 57 ++++++++++++++++++++----------- SRC/dlantp.f | 55 +++++++++++++++++++++--------- SRC/dlantr.f | 56 ++++++++++++++++++++---------- SRC/scombssq.f | 92 ++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/slangb.f | 26 +++++++++----- SRC/slange.f | 22 ++++++++---- SRC/slanhs.f | 25 ++++++++++---- SRC/slansb.f | 44 ++++++++++++++++++------ SRC/slansp.f | 48 +++++++++++++++++++------- SRC/slansy.f | 42 +++++++++++++++++------ SRC/slantb.f | 57 ++++++++++++++++++++----------- SRC/slantp.f | 55 +++++++++++++++++++++--------- SRC/slantr.f | 56 ++++++++++++++++++++---------- SRC/zlangb.f | 23 +++++++++---- SRC/zlange.f | 22 ++++++++---- SRC/zlanhb.f | 48 +++++++++++++++++++------- SRC/zlanhe.f | 45 +++++++++++++++++------- SRC/zlanhp.f | 46 ++++++++++++++++++------- SRC/zlanhs.f | 23 +++++++++---- SRC/zlansb.f | 42 +++++++++++++++++------ SRC/zlansp.f | 54 ++++++++++++++++++++--------- SRC/zlansy.f | 40 +++++++++++++++++----- SRC/zlantb.f | 55 ++++++++++++++++++++---------- SRC/zlantp.f | 53 +++++++++++++++++++++-------- SRC/zlantr.f | 54 ++++++++++++++++++++--------- 45 files changed, 1454 insertions(+), 494 deletions(-) create mode 100644 SRC/dcombssq.f create mode 100644 SRC/scombssq.f diff --git a/SRC/Makefile b/SRC/Makefile index 1814521999..f583115ca7 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -174,7 +174,7 @@ SLASRC = \ ssytrd_2stage.o ssytrd_sy2sb.o ssytrd_sb2st.o ssb2st_kernels.o \ ssyevd_2stage.o ssyev_2stage.o ssyevx_2stage.o ssyevr_2stage.o \ ssbev_2stage.o ssbevx_2stage.o ssbevd_2stage.o ssygv_2stage.o \ - sgesvdq.o + sgesvdq.o scombssq.o DSLASRC = spotrs.o sgetrs.o spotrf.o sgetrf.o @@ -373,7 +373,7 @@ DLASRC = \ dsytrd_2stage.o dsytrd_sy2sb.o dsytrd_sb2st.o dsb2st_kernels.o \ dsyevd_2stage.o dsyev_2stage.o dsyevx_2stage.o dsyevr_2stage.o \ dsbev_2stage.o dsbevx_2stage.o dsbevd_2stage.o dsygv_2stage.o \ - dgesvdq.o + dgesvdq.o dcombssq.o ifdef USEXBLAS DXLASRC = dgesvxx.o dgerfsx.o dla_gerfsx_extended.o dla_geamv.o \ diff --git a/SRC/clangb.f b/SRC/clangb.f index 14a163ea7d..9818360fe4 100644 --- a/SRC/clangb.f +++ b/SRC/clangb.f @@ -130,6 +130,7 @@ REAL FUNCTION CLANGB( NORM, N, KL, KU, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER KL, KU, LDAB, N @@ -147,14 +148,17 @@ REAL FUNCTION CLANGB( NORM, N, KL, KU, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, K, L - REAL SCALE, SUM, VALUE, TEMP + REAL SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -207,15 +211,22 @@ REAL FUNCTION CLANGB( NORM, N, KL, KU, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N L = MAX( 1, J-KU ) K = KU + 1 - J + L - CALL CLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANGB = VALUE diff --git a/SRC/clange.f b/SRC/clange.f index 50f705a18d..00895c8bc7 100644 --- a/SRC/clange.f +++ b/SRC/clange.f @@ -120,6 +120,7 @@ REAL FUNCTION CLANGE( NORM, M, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, M, N @@ -137,14 +138,17 @@ REAL FUNCTION CLANGE( NORM, M, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL SCALE, SUM, VALUE, TEMP + REAL SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -196,13 +200,19 @@ REAL FUNCTION CLANGE( NORM, M, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL CLASSQ( M, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANGE = VALUE diff --git a/SRC/clanhb.f b/SRC/clanhb.f index 2b034b19b5..f78de23df5 100644 --- a/SRC/clanhb.f +++ b/SRC/clanhb.f @@ -137,6 +137,7 @@ REAL FUNCTION CLANHB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -154,14 +155,17 @@ REAL FUNCTION CLANHB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, REAL, SQRT @@ -233,39 +237,57 @@ REAL FUNCTION CLANHB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL CLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 J = 1, N IF( REAL( AB( L, J ) ).NE.ZERO ) THEN ABSA = ABS( REAL( AB( L, J ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANHB = VALUE diff --git a/SRC/clanhe.f b/SRC/clanhe.f index 101d778eb0..33d6c8b014 100644 --- a/SRC/clanhe.f +++ b/SRC/clanhe.f @@ -129,6 +129,7 @@ REAL FUNCTION CLANHE( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -146,14 +147,17 @@ REAL FUNCTION CLANHE( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, REAL, SQRT @@ -223,31 +227,48 @@ REAL FUNCTION CLANHE( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL CLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J-1, A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J, A( J+1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* DO 130 I = 1, N IF( REAL( A( I, I ) ).NE.ZERO ) THEN ABSA = ABS( REAL( A( I, I ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( SSQ( 1 ).LT.ABSA ) THEN + SSQ( 2 ) = ONE + SSQ( 2 )*( SSQ( 1 ) / ABSA )**2 + SSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + SSQ( 2 ) = SSQ( 2 ) + ( ABSA / SSQ( 1 ) )**2 END IF END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANHE = VALUE diff --git a/SRC/clanhp.f b/SRC/clanhp.f index c8927d5033..e0e23abc7e 100644 --- a/SRC/clanhp.f +++ b/SRC/clanhp.f @@ -122,6 +122,7 @@ REAL FUNCTION CLANHP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -139,14 +140,17 @@ REAL FUNCTION CLANHP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, REAL, SQRT @@ -225,31 +229,48 @@ REAL FUNCTION CLANHP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL CLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( REAL( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( REAL( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -258,7 +279,8 @@ REAL FUNCTION CLANHP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANHP = VALUE diff --git a/SRC/clanhs.f b/SRC/clanhs.f index 35623b73d7..661b4f9013 100644 --- a/SRC/clanhs.f +++ b/SRC/clanhs.f @@ -114,6 +114,7 @@ REAL FUNCTION CLANHS( NORM, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, N @@ -131,14 +132,17 @@ REAL FUNCTION CLANHS( NORM, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL SCALE, SUM, VALUE + REAL SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -190,13 +194,20 @@ REAL FUNCTION CLANHS( NORM, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL CLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N, J+1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANHS = VALUE diff --git a/SRC/clansb.f b/SRC/clansb.f index fbc50674c0..1085fd8805 100644 --- a/SRC/clansb.f +++ b/SRC/clansb.f @@ -135,6 +135,7 @@ REAL FUNCTION CLANSB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -152,14 +153,17 @@ REAL FUNCTION CLANSB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -227,29 +231,47 @@ REAL FUNCTION CLANSB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL CLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF - CALL CLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANSB = VALUE diff --git a/SRC/clansp.f b/SRC/clansp.f index fd64366c6e..628dc0a75f 100644 --- a/SRC/clansp.f +++ b/SRC/clansp.f @@ -120,6 +120,7 @@ REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -137,14 +138,17 @@ REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, AIMAG, REAL, SQRT @@ -219,40 +223,57 @@ REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL CLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( REAL( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( REAL( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( AIMAG( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( AIMAG( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -261,7 +282,8 @@ REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANSP = VALUE diff --git a/SRC/clansy.f b/SRC/clansy.f index 3aa787410c..537fb7ba9d 100644 --- a/SRC/clansy.f +++ b/SRC/clansy.f @@ -128,6 +128,7 @@ REAL FUNCTION CLANSY( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -145,14 +146,17 @@ REAL FUNCTION CLANSY( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT @@ -218,21 +222,39 @@ REAL FUNCTION CLANSY( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL CLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL CLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM - CALL CLASSQ( N, A, LDA+1, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANSY = VALUE diff --git a/SRC/clantb.f b/SRC/clantb.f index 4b4361c796..8066d0ef69 100644 --- a/SRC/clantb.f +++ b/SRC/clantb.f @@ -146,6 +146,7 @@ REAL FUNCTION CLANTB( NORM, UPLO, DIAG, N, K, AB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER K, LDAB, N @@ -164,14 +165,17 @@ REAL FUNCTION CLANTB( NORM, UPLO, DIAG, N, K, AB, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, L - REAL SCALE, SUM, VALUE + REAL SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -313,46 +317,61 @@ REAL FUNCTION CLANTB( NORM, UPLO, DIAG, N, K, AB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 280 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL CLASSQ( MIN( J-1, K ), - $ AB( MAX( K+2-J, 1 ), J ), 1, SCALE, - $ SUM ) + $ AB( MAX( K+2-J, 1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 280 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 290 J = 1, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL CLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 300 J = 1, N - 1 - CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 310 J = 1, N - CALL CLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANTB = VALUE diff --git a/SRC/clantp.f b/SRC/clantp.f index 148ac5436b..b0c48eb468 100644 --- a/SRC/clantp.f +++ b/SRC/clantp.f @@ -130,6 +130,7 @@ REAL FUNCTION CLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER N @@ -148,14 +149,17 @@ REAL FUNCTION CLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, K - REAL SCALE, SUM, VALUE + REAL SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT @@ -308,45 +312,64 @@ REAL FUNCTION CLANTP( NORM, UPLO, DIAG, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 280 J = 2, N - CALL CLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J-1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 280 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 290 J = 1, N - CALL CLASSQ( J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 300 J = 1, N - 1 - CALL CLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 300 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 310 J = 1, N - CALL CLASSQ( N-J+1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( N-J+1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANTP = VALUE diff --git a/SRC/clantr.f b/SRC/clantr.f index 4e1843d3dc..30ed649250 100644 --- a/SRC/clantr.f +++ b/SRC/clantr.f @@ -147,6 +147,7 @@ REAL FUNCTION CLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER LDA, M, N @@ -165,14 +166,17 @@ REAL FUNCTION CLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J - REAL SCALE, SUM, VALUE + REAL SUM, VALUE +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. * .. External Subroutines .. - EXTERNAL CLASSQ + EXTERNAL CLASSQ, SCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -313,38 +317,56 @@ REAL FUNCTION CLANTR( NORM, UPLO, DIAG, M, N, A, LDA, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 290 J = 2, N - CALL CLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( M, J-1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 300 J = 1, N - CALL CLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( MIN( M, J ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 310 J = 1, N - CALL CLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 320 J = 1, N - CALL CLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL CLASSQ( M-J+1, A( J, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 320 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * CLANTR = VALUE diff --git a/SRC/dcombssq.f b/SRC/dcombssq.f new file mode 100644 index 0000000000..79f6d95c92 --- /dev/null +++ b/SRC/dcombssq.f @@ -0,0 +1,92 @@ +*> \brief \b DCOMBSSQ adds two scaled sum of squares quantities. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* +* Definition: +* =========== +* +* SUBROUTINE DCOMBSSQ( V1, V2 ) +* +* .. Array Arguments .. +* DOUBLE PRECISION V1( 2 ), V2( 2 ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2. +*> That is, +*> +*> V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq +*> + V2_scale**2 * V2_sumsq +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in,out] V1 +*> \verbatim +*> V1 is DOUBLE PRECISION array, dimension (2). +*> The first scaled sum. +*> V1(1) = V1_scale, V1(2) = V1_sumsq. +*> \endverbatim +*> +*> \param[in] V2 +*> \verbatim +*> V2 is DOUBLE PRECISION array, dimension (2). +*> The second scaled sum. +*> V2(1) = V2_scale, V2(2) = V2_sumsq. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2018 +* +*> \ingroup OTHERauxiliary +* +* ===================================================================== + SUBROUTINE DCOMBSSQ( V1, V2 ) +* +* -- LAPACK auxiliary routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2018 +* +* .. Array Arguments .. + DOUBLE PRECISION V1( 2 ), V2( 2 ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) +* .. +* .. Executable Statements .. +* + IF( V1( 1 ).GE.V2( 1 ) ) THEN + IF( V1( 1 ).NE.ZERO ) THEN + V1( 2 ) = V1( 2 ) + ( V2( 1 ) / V1( 1 ) )**2 * V2( 2 ) + END IF + ELSE + V1( 2 ) = V2( 2 ) + ( V1( 1 ) / V2( 1 ) )**2 * V1( 2 ) + V1( 1 ) = V2( 1 ) + END IF + RETURN +* +* End of DCOMBSSQ +* + END diff --git a/SRC/dlangb.f b/SRC/dlangb.f index 078573b87a..0c4f938f7f 100644 --- a/SRC/dlangb.f +++ b/SRC/dlangb.f @@ -129,6 +129,7 @@ DOUBLE PRECISION FUNCTION DLANGB( NORM, N, KL, KU, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER KL, KU, LDAB, N @@ -139,22 +140,24 @@ DOUBLE PRECISION FUNCTION DLANGB( NORM, N, KL, KU, AB, LDAB, * * ===================================================================== * -* * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I, J, K, L - DOUBLE PRECISION SCALE, SUM, VALUE, TEMP + DOUBLE PRECISION SUM, VALUE, TEMP * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -206,15 +209,22 @@ DOUBLE PRECISION FUNCTION DLANGB( NORM, N, KL, KU, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N L = MAX( 1, J-KU ) K = KU + 1 - J + L - CALL DLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANGB = VALUE diff --git a/SRC/dlange.f b/SRC/dlange.f index 9dbf45e818..6b32fbefd7 100644 --- a/SRC/dlange.f +++ b/SRC/dlange.f @@ -119,6 +119,7 @@ DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, M, N @@ -135,10 +136,13 @@ DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE, TEMP + DOUBLE PRECISION SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Subroutines .. - EXTERNAL DLASSQ + EXTERNAL DLASSQ, DCOMBSSQ * .. * .. External Functions .. LOGICAL LSAME, DISNAN @@ -194,13 +198,19 @@ DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL DLASSQ( M, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANGE = VALUE diff --git a/SRC/dlanhs.f b/SRC/dlanhs.f index 691dbc21ec..a859d22164 100644 --- a/SRC/dlanhs.f +++ b/SRC/dlanhs.f @@ -113,6 +113,7 @@ DOUBLE PRECISION FUNCTION DLANHS( NORM, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, N @@ -129,15 +130,18 @@ DOUBLE PRECISION FUNCTION DLANHS( NORM, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT * .. @@ -188,13 +192,20 @@ DOUBLE PRECISION FUNCTION DLANHS( NORM, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL DLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( N, J+1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANHS = VALUE diff --git a/SRC/dlansb.f b/SRC/dlansb.f index 4ccf5f27e1..a82dc41b1f 100644 --- a/SRC/dlansb.f +++ b/SRC/dlansb.f @@ -134,6 +134,7 @@ DOUBLE PRECISION FUNCTION DLANSB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -150,15 +151,18 @@ DOUBLE PRECISION FUNCTION DLANSB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -225,29 +229,47 @@ DOUBLE PRECISION FUNCTION DLANSB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL DLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL DLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF - CALL DLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANSB = VALUE diff --git a/SRC/dlansp.f b/SRC/dlansp.f index a1829db75c..b6ad1ffcf0 100644 --- a/SRC/dlansp.f +++ b/SRC/dlansp.f @@ -119,6 +119,7 @@ DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -135,15 +136,18 @@ DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -217,31 +221,48 @@ DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL DLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL DLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( AP( K ).NE.ZERO ) THEN ABSA = ABS( AP( K ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -250,7 +271,8 @@ DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANSP = VALUE diff --git a/SRC/dlansy.f b/SRC/dlansy.f index 2372fce0a8..87d514c118 100644 --- a/SRC/dlansy.f +++ b/SRC/dlansy.f @@ -127,6 +127,7 @@ DOUBLE PRECISION FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -143,15 +144,18 @@ DOUBLE PRECISION FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -216,21 +220,39 @@ DOUBLE PRECISION FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL DLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL DLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM - CALL DLASSQ( N, A, LDA+1, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANSY = VALUE diff --git a/SRC/dlantb.f b/SRC/dlantb.f index 3d2bfe7e4b..0d46f6cc89 100644 --- a/SRC/dlantb.f +++ b/SRC/dlantb.f @@ -145,6 +145,7 @@ DOUBLE PRECISION FUNCTION DLANTB( NORM, UPLO, DIAG, N, K, AB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER K, LDAB, N @@ -162,15 +163,18 @@ DOUBLE PRECISION FUNCTION DLANTB( NORM, UPLO, DIAG, N, K, AB, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, L - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -311,46 +315,61 @@ DOUBLE PRECISION FUNCTION DLANTB( NORM, UPLO, DIAG, N, K, AB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 280 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL DLASSQ( MIN( J-1, K ), - $ AB( MAX( K+2-J, 1 ), J ), 1, SCALE, - $ SUM ) + $ AB( MAX( K+2-J, 1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 280 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 290 J = 1, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL DLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 300 J = 1, N - 1 - CALL DLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 310 J = 1, N - CALL DLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANTB = VALUE diff --git a/SRC/dlantp.f b/SRC/dlantp.f index f84a9e9d7d..a7b89dec75 100644 --- a/SRC/dlantp.f +++ b/SRC/dlantp.f @@ -129,6 +129,7 @@ DOUBLE PRECISION FUNCTION DLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER N @@ -146,15 +147,18 @@ DOUBLE PRECISION FUNCTION DLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, K - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -306,45 +310,64 @@ DOUBLE PRECISION FUNCTION DLANTP( NORM, UPLO, DIAG, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 280 J = 2, N - CALL DLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( J-1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 280 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 290 J = 1, N - CALL DLASSQ( J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 300 J = 1, N - 1 - CALL DLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N-J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 300 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 310 J = 1, N - CALL DLASSQ( N-J+1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( N-J+1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANTP = VALUE diff --git a/SRC/dlantr.f b/SRC/dlantr.f index 8585b2f689..b0f85ef4ed 100644 --- a/SRC/dlantr.f +++ b/SRC/dlantr.f @@ -146,6 +146,7 @@ DOUBLE PRECISION FUNCTION DLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER LDA, M, N @@ -163,15 +164,18 @@ DOUBLE PRECISION FUNCTION DLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL DLASSQ +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. +* .. External Subroutines .. + EXTERNAL DLASSQ, DCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT * .. @@ -311,38 +315,56 @@ DOUBLE PRECISION FUNCTION DLANTR( NORM, UPLO, DIAG, M, N, A, LDA, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 290 J = 2, N - CALL DLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( M, J-1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 300 J = 1, N - CALL DLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( MIN( M, J ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 310 J = 1, N - CALL DLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 320 J = 1, N - CALL DLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL DLASSQ( M-J+1, A( J, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 320 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * DLANTR = VALUE diff --git a/SRC/scombssq.f b/SRC/scombssq.f new file mode 100644 index 0000000000..76bc0e320e --- /dev/null +++ b/SRC/scombssq.f @@ -0,0 +1,92 @@ +*> \brief \b SCOMBSSQ adds two scaled sum of squares quantities +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* +* Definition: +* =========== +* +* SUBROUTINE SCOMBSSQ( V1, V2 ) +* +* .. Array Arguments .. +* REAL V1( 2 ), V2( 2 ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2. +*> That is, +*> +*> V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq +*> + V2_scale**2 * V2_sumsq +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in,out] V1 +*> \verbatim +*> V1 is REAL array, dimension (2). +*> The first scaled sum. +*> V1(1) = V1_scale, V1(2) = V1_sumsq. +*> \endverbatim +*> +*> \param[in] V2 +*> \verbatim +*> V2 is REAL array, dimension (2). +*> The second scaled sum. +*> V2(1) = V2_scale, V2(2) = V2_sumsq. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2018 +* +*> \ingroup OTHERauxiliary +* +* ===================================================================== + SUBROUTINE SCOMBSSQ( V1, V2 ) +* +* -- LAPACK auxiliary routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2018 +* +* .. Array Arguments .. + REAL V1( 2 ), V2( 2 ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0D+0 ) +* .. +* .. Executable Statements .. +* + IF( V1( 1 ).GE.V2( 1 ) ) THEN + IF( V1( 1 ).NE.ZERO ) THEN + V1( 2 ) = V1( 2 ) + ( V2( 1 ) / V1( 1 ) )**2 * V2( 2 ) + END IF + ELSE + V1( 2 ) = V2( 2 ) + ( V1( 1 ) / V2( 1 ) )**2 * V1( 2 ) + V1( 1 ) = V2( 1 ) + END IF + RETURN +* +* End of SCOMBSSQ +* + END diff --git a/SRC/slangb.f b/SRC/slangb.f index fd538b1b70..706e07501c 100644 --- a/SRC/slangb.f +++ b/SRC/slangb.f @@ -129,6 +129,7 @@ REAL FUNCTION SLANGB( NORM, N, KL, KU, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER KL, KU, LDAB, N @@ -139,22 +140,24 @@ REAL FUNCTION SLANGB( NORM, N, KL, KU, AB, LDAB, * * ===================================================================== * -* * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. INTEGER I, J, K, L - REAL SCALE, SUM, VALUE, TEMP + REAL SUM, VALUE, TEMP * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -206,15 +209,22 @@ REAL FUNCTION SLANGB( NORM, N, KL, KU, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N L = MAX( 1, J-KU ) K = KU + 1 - J + L - CALL SLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANGB = VALUE diff --git a/SRC/slange.f b/SRC/slange.f index 2eb8d7d140..0c80f1d40a 100644 --- a/SRC/slange.f +++ b/SRC/slange.f @@ -119,6 +119,7 @@ REAL FUNCTION SLANGE( NORM, M, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, M, N @@ -135,10 +136,13 @@ REAL FUNCTION SLANGE( NORM, M, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL SCALE, SUM, VALUE, TEMP + REAL SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Subroutines .. - EXTERNAL SLASSQ + EXTERNAL SLASSQ, SCOMBSSQ * .. * .. External Functions .. LOGICAL LSAME, SISNAN @@ -194,13 +198,19 @@ REAL FUNCTION SLANGE( NORM, M, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL SLASSQ( M, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANGE = VALUE diff --git a/SRC/slanhs.f b/SRC/slanhs.f index c5a077fbf1..8913031a21 100644 --- a/SRC/slanhs.f +++ b/SRC/slanhs.f @@ -113,6 +113,7 @@ REAL FUNCTION SLANHS( NORM, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, N @@ -129,15 +130,18 @@ REAL FUNCTION SLANHS( NORM, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL SCALE, SUM, VALUE + REAL SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT * .. @@ -188,13 +192,20 @@ REAL FUNCTION SLANHS( NORM, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL SLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( N, J+1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANHS = VALUE diff --git a/SRC/slansb.f b/SRC/slansb.f index 8f3fe1eb95..23519025de 100644 --- a/SRC/slansb.f +++ b/SRC/slansb.f @@ -134,6 +134,7 @@ REAL FUNCTION SLANSB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -150,15 +151,18 @@ REAL FUNCTION SLANSB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -225,29 +229,47 @@ REAL FUNCTION SLANSB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL SLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF - CALL SLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANSB = VALUE diff --git a/SRC/slansp.f b/SRC/slansp.f index 35390cd1ca..7e29d778be 100644 --- a/SRC/slansp.f +++ b/SRC/slansp.f @@ -119,6 +119,7 @@ REAL FUNCTION SLANSP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -135,15 +136,18 @@ REAL FUNCTION SLANSP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -217,31 +221,48 @@ REAL FUNCTION SLANSP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL SLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL SLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( AP( K ).NE.ZERO ) THEN ABSA = ABS( AP( K ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -250,7 +271,8 @@ REAL FUNCTION SLANSP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANSP = VALUE diff --git a/SRC/slansy.f b/SRC/slansy.f index c8400e5308..66ff1c5c76 100644 --- a/SRC/slansy.f +++ b/SRC/slansy.f @@ -127,6 +127,7 @@ REAL FUNCTION SLANSY( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -143,15 +144,18 @@ REAL FUNCTION SLANSY( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - REAL ABSA, SCALE, SUM, VALUE + REAL ABSA, SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -216,21 +220,39 @@ REAL FUNCTION SLANSY( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL SLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL SLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM - CALL SLASSQ( N, A, LDA+1, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANSY = VALUE diff --git a/SRC/slantb.f b/SRC/slantb.f index 3588779cb6..5b94618e19 100644 --- a/SRC/slantb.f +++ b/SRC/slantb.f @@ -145,6 +145,7 @@ REAL FUNCTION SLANTB( NORM, UPLO, DIAG, N, K, AB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER K, LDAB, N @@ -162,15 +163,18 @@ REAL FUNCTION SLANTB( NORM, UPLO, DIAG, N, K, AB, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, L - REAL SCALE, SUM, VALUE + REAL SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. @@ -311,46 +315,61 @@ REAL FUNCTION SLANTB( NORM, UPLO, DIAG, N, K, AB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 280 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL SLASSQ( MIN( J-1, K ), - $ AB( MAX( K+2-J, 1 ), J ), 1, SCALE, - $ SUM ) + $ AB( MAX( K+2-J, 1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 280 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 290 J = 1, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL SLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 300 J = 1, N - 1 - CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 310 J = 1, N - CALL SLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANTB = VALUE diff --git a/SRC/slantp.f b/SRC/slantp.f index 1423f5ca36..ab781deac3 100644 --- a/SRC/slantp.f +++ b/SRC/slantp.f @@ -129,6 +129,7 @@ REAL FUNCTION SLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER N @@ -146,15 +147,18 @@ REAL FUNCTION SLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, K - REAL SCALE, SUM, VALUE + REAL SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. @@ -306,45 +310,64 @@ REAL FUNCTION SLANTP( NORM, UPLO, DIAG, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 280 J = 2, N - CALL SLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( J-1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 280 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 290 J = 1, N - CALL SLASSQ( J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + J 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 300 J = 1, N - 1 - CALL SLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N-J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 300 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 310 J = 1, N - CALL SLASSQ( N-J+1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( N-J+1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANTP = VALUE diff --git a/SRC/slantr.f b/SRC/slantr.f index 63b8558926..8168a598cb 100644 --- a/SRC/slantr.f +++ b/SRC/slantr.f @@ -146,6 +146,7 @@ REAL FUNCTION SLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER LDA, M, N @@ -163,15 +164,18 @@ REAL FUNCTION SLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J - REAL SCALE, SUM, VALUE + REAL SUM, VALUE * .. -* .. External Subroutines .. - EXTERNAL SLASSQ +* .. Local Arrays .. + REAL SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, SISNAN EXTERNAL LSAME, SISNAN * .. +* .. External Subroutines .. + EXTERNAL SLASSQ, SCOMBSSQ +* .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT * .. @@ -311,38 +315,56 @@ REAL FUNCTION SLANTR( NORM, UPLO, DIAG, M, N, A, LDA, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 290 J = 2, N - CALL SLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( M, J-1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 300 J = 1, N - CALL SLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( MIN( M, J ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 310 J = 1, N - CALL SLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 320 J = 1, N - CALL SLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL SLASSQ( M-J+1, A( J, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL SCOMBSSQ( SSQ, COLSSQ ) 320 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * SLANTR = VALUE diff --git a/SRC/zlangb.f b/SRC/zlangb.f index 949bb2c01d..e40a470fde 100644 --- a/SRC/zlangb.f +++ b/SRC/zlangb.f @@ -130,6 +130,7 @@ DOUBLE PRECISION FUNCTION ZLANGB( NORM, N, KL, KU, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER KL, KU, LDAB, N @@ -147,14 +148,17 @@ DOUBLE PRECISION FUNCTION ZLANGB( NORM, N, KL, KU, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, K, L - DOUBLE PRECISION SCALE, SUM, VALUE, TEMP + DOUBLE PRECISION SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -207,15 +211,22 @@ DOUBLE PRECISION FUNCTION ZLANGB( NORM, N, KL, KU, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N L = MAX( 1, J-KU ) K = KU + 1 - J + L - CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANGB = VALUE diff --git a/SRC/zlange.f b/SRC/zlange.f index 5407decef9..8162786fb7 100644 --- a/SRC/zlange.f +++ b/SRC/zlange.f @@ -120,6 +120,7 @@ DOUBLE PRECISION FUNCTION ZLANGE( NORM, M, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, M, N @@ -137,14 +138,17 @@ DOUBLE PRECISION FUNCTION ZLANGE( NORM, M, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE, TEMP + DOUBLE PRECISION SUM, VALUE, TEMP +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -196,13 +200,19 @@ DOUBLE PRECISION FUNCTION ZLANGE( NORM, M, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL ZLASSQ( M, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANGE = VALUE diff --git a/SRC/zlanhb.f b/SRC/zlanhb.f index b3717804f9..16b5c117c5 100644 --- a/SRC/zlanhb.f +++ b/SRC/zlanhb.f @@ -137,6 +137,7 @@ DOUBLE PRECISION FUNCTION ZLANHB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -154,14 +155,17 @@ DOUBLE PRECISION FUNCTION ZLANHB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, MAX, MIN, SQRT @@ -233,39 +237,57 @@ DOUBLE PRECISION FUNCTION ZLANHB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 J = 1, N IF( DBLE( AB( L, J ) ).NE.ZERO ) THEN ABSA = ABS( DBLE( AB( L, J ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANHB = VALUE diff --git a/SRC/zlanhe.f b/SRC/zlanhe.f index 7c7f7f3be4..5aef9a7569 100644 --- a/SRC/zlanhe.f +++ b/SRC/zlanhe.f @@ -129,6 +129,7 @@ DOUBLE PRECISION FUNCTION ZLANHE( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -146,14 +147,17 @@ DOUBLE PRECISION FUNCTION ZLANHE( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, SQRT @@ -223,31 +227,48 @@ DOUBLE PRECISION FUNCTION ZLANHE( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J-1, A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J, A( J+1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* DO 130 I = 1, N IF( DBLE( A( I, I ) ).NE.ZERO ) THEN ABSA = ABS( DBLE( A( I, I ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( SSQ( 1 ).LT.ABSA ) THEN + SSQ( 2 ) = ONE + SSQ( 2 )*( SSQ( 1 ) / ABSA )**2 + SSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + SSQ( 2 ) = SSQ( 2 ) + ( ABSA / SSQ( 1 ) )**2 END IF END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANHE = VALUE diff --git a/SRC/zlanhp.f b/SRC/zlanhp.f index 9ded607460..d795aeca9d 100644 --- a/SRC/zlanhp.f +++ b/SRC/zlanhp.f @@ -122,6 +122,7 @@ DOUBLE PRECISION FUNCTION ZLANHP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -139,14 +140,17 @@ DOUBLE PRECISION FUNCTION ZLANHP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, SQRT @@ -225,31 +229,48 @@ DOUBLE PRECISION FUNCTION ZLANHP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( DBLE( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( DBLE( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -258,7 +279,8 @@ DOUBLE PRECISION FUNCTION ZLANHP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANHP = VALUE diff --git a/SRC/zlanhs.f b/SRC/zlanhs.f index f2d36b3042..bd8e86be9d 100644 --- a/SRC/zlanhs.f +++ b/SRC/zlanhs.f @@ -114,6 +114,7 @@ DOUBLE PRECISION FUNCTION ZLANHS( NORM, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM INTEGER LDA, N @@ -131,14 +132,17 @@ DOUBLE PRECISION FUNCTION ZLANHS( NORM, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -190,13 +194,20 @@ DOUBLE PRECISION FUNCTION ZLANHS( NORM, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 90 J = 1, N - CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 90 CONTINUE - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANHS = VALUE diff --git a/SRC/zlansb.f b/SRC/zlansb.f index 3468c49b3a..245dcaf4b7 100644 --- a/SRC/zlansb.f +++ b/SRC/zlansb.f @@ -135,6 +135,7 @@ DOUBLE PRECISION FUNCTION ZLANSB( NORM, UPLO, N, K, AB, LDAB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER K, LDAB, N @@ -152,14 +153,17 @@ DOUBLE PRECISION FUNCTION ZLANSB( NORM, UPLO, N, K, AB, LDAB, * .. * .. Local Scalars .. INTEGER I, J, L - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -227,29 +231,47 @@ DOUBLE PRECISION FUNCTION ZLANSB( NORM, UPLO, N, K, AB, LDAB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( K.GT.0 ) THEN IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE L = K + 1 ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE L = 1 END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) ELSE L = 1 END IF - CALL ZLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANSB = VALUE diff --git a/SRC/zlansp.f b/SRC/zlansp.f index 84fb972bbe..fa92204874 100644 --- a/SRC/zlansp.f +++ b/SRC/zlansp.f @@ -120,6 +120,7 @@ DOUBLE PRECISION FUNCTION ZLANSP( NORM, UPLO, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER N @@ -137,14 +138,17 @@ DOUBLE PRECISION FUNCTION ZLANSP( NORM, UPLO, N, AP, WORK ) * .. * .. Local Scalars .. INTEGER I, J, K - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DIMAG, SQRT @@ -219,40 +223,57 @@ DOUBLE PRECISION FUNCTION ZLANSP( NORM, UPLO, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE K = 2 IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 120 CONTINUE END IF - SUM = 2*SUM + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* K = 1 + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE DO 130 I = 1, N IF( DBLE( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( DBLE( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( DIMAG( AP( K ) ).NE.ZERO ) THEN ABSA = ABS( DIMAG( AP( K ) ) ) - IF( SCALE.LT.ABSA ) THEN - SUM = ONE + SUM*( SCALE / ABSA )**2 - SCALE = ABSA + IF( COLSSQ( 1 ).LT.ABSA ) THEN + COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2 + COLSSQ( 1 ) = ABSA ELSE - SUM = SUM + ( ABSA / SCALE )**2 + COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2 END IF END IF IF( LSAME( UPLO, 'U' ) ) THEN @@ -261,7 +282,8 @@ DOUBLE PRECISION FUNCTION ZLANSP( NORM, UPLO, N, AP, WORK ) K = K + N - I + 1 END IF 130 CONTINUE - VALUE = SCALE*SQRT( SUM ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANSP = VALUE diff --git a/SRC/zlansy.f b/SRC/zlansy.f index 58269a911f..e022f85e1e 100644 --- a/SRC/zlansy.f +++ b/SRC/zlansy.f @@ -128,6 +128,7 @@ DOUBLE PRECISION FUNCTION ZLANSY( NORM, UPLO, N, A, LDA, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER NORM, UPLO INTEGER LDA, N @@ -145,14 +146,17 @@ DOUBLE PRECISION FUNCTION ZLANSY( NORM, UPLO, N, A, LDA, WORK ) * .. * .. Local Scalars .. INTEGER I, J - DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + DOUBLE PRECISION ABSA, SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT @@ -218,21 +222,39 @@ DOUBLE PRECISION FUNCTION ZLANSY( NORM, UPLO, N, A, LDA, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. +* + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE +* +* Sum off-diagonals * - SCALE = ZERO - SUM = ONE IF( LSAME( UPLO, 'U' ) ) THEN DO 110 J = 2, N - CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 110 CONTINUE ELSE DO 120 J = 1, N - 1 - CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 120 CONTINUE END IF - SUM = 2*SUM - CALL ZLASSQ( N, A, LDA+1, SCALE, SUM ) - VALUE = SCALE*SQRT( SUM ) + SSQ( 2 ) = 2*SSQ( 2 ) +* +* Sum diagonal +* + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANSY = VALUE diff --git a/SRC/zlantb.f b/SRC/zlantb.f index 3077ba1515..f02509223d 100644 --- a/SRC/zlantb.f +++ b/SRC/zlantb.f @@ -146,6 +146,7 @@ DOUBLE PRECISION FUNCTION ZLANTB( NORM, UPLO, DIAG, N, K, AB, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER K, LDAB, N @@ -164,14 +165,17 @@ DOUBLE PRECISION FUNCTION ZLANTB( NORM, UPLO, DIAG, N, K, AB, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, L - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -313,46 +317,61 @@ DOUBLE PRECISION FUNCTION ZLANTB( NORM, UPLO, DIAG, N, K, AB, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 280 J = 2, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL ZLASSQ( MIN( J-1, K ), - $ AB( MAX( K+2-J, 1 ), J ), 1, SCALE, - $ SUM ) + $ AB( MAX( K+2-J, 1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 280 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 290 J = 1, N + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE CALL ZLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ), - $ 1, SCALE, SUM ) + $ 1, COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N IF( K.GT.0 ) THEN DO 300 J = 1, N - 1 - CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 310 J = 1, N - CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANTB = VALUE diff --git a/SRC/zlantp.f b/SRC/zlantp.f index 69dbaa5bca..d32a00f138 100644 --- a/SRC/zlantp.f +++ b/SRC/zlantp.f @@ -130,6 +130,7 @@ DOUBLE PRECISION FUNCTION ZLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER N @@ -148,14 +149,17 @@ DOUBLE PRECISION FUNCTION ZLANTP( NORM, UPLO, DIAG, N, AP, WORK ) * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J, K - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT @@ -308,45 +312,64 @@ DOUBLE PRECISION FUNCTION ZLANTP( NORM, UPLO, DIAG, N, AP, WORK ) ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 280 J = 2, N - CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J-1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 280 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 290 J = 1, N - CALL ZLASSQ( J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + J 290 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = N + SSQ( 1 ) = ONE + SSQ( 2 ) = N K = 2 DO 300 J = 1, N - 1 - CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 300 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE K = 1 DO 310 J = 1, N - CALL ZLASSQ( N-J+1, AP( K ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( N-J+1, AP( K ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) K = K + N - J + 1 310 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANTP = VALUE diff --git a/SRC/zlantr.f b/SRC/zlantr.f index 04ee482f7d..e763fe5f70 100644 --- a/SRC/zlantr.f +++ b/SRC/zlantr.f @@ -147,6 +147,7 @@ DOUBLE PRECISION FUNCTION ZLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * December 2016 * + IMPLICIT NONE * .. Scalar Arguments .. CHARACTER DIAG, NORM, UPLO INTEGER LDA, M, N @@ -165,14 +166,17 @@ DOUBLE PRECISION FUNCTION ZLANTR( NORM, UPLO, DIAG, M, N, A, LDA, * .. Local Scalars .. LOGICAL UDIAG INTEGER I, J - DOUBLE PRECISION SCALE, SUM, VALUE + DOUBLE PRECISION SUM, VALUE +* .. +* .. Local Arrays .. + DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 ) * .. * .. External Functions .. LOGICAL LSAME, DISNAN EXTERNAL LSAME, DISNAN * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, DCOMBSSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MIN, SQRT @@ -313,38 +317,56 @@ DOUBLE PRECISION FUNCTION ZLANTR( NORM, UPLO, DIAG, M, N, A, LDA, ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). +* SSQ(1) is scale +* SSQ(2) is sum-of-squares +* For better accuracy, sum each column separately. * IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 290 J = 2, N - CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 290 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 300 J = 1, N - CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 300 CONTINUE END IF ELSE IF( LSAME( DIAG, 'U' ) ) THEN - SCALE = ONE - SUM = MIN( M, N ) + SSQ( 1 ) = ONE + SSQ( 2 ) = MIN( M, N ) DO 310 J = 1, N - CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, - $ SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 310 CONTINUE ELSE - SCALE = ZERO - SUM = ONE + SSQ( 1 ) = ZERO + SSQ( 2 ) = ONE DO 320 J = 1, N - CALL ZLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) + COLSSQ( 1 ) = ZERO + COLSSQ( 2 ) = ONE + CALL ZLASSQ( M-J+1, A( J, J ), 1, + $ COLSSQ( 1 ), COLSSQ( 2 ) ) + CALL DCOMBSSQ( SSQ, COLSSQ ) 320 CONTINUE END IF END IF - VALUE = SCALE*SQRT( SUM ) + VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) ) END IF * ZLANTR = VALUE From 0d9289b1ce52778fc7798e4c3991e7732df3e7bb Mon Sep 17 00:00:00 2001 From: Mark Gates Date: Mon, 5 Nov 2018 21:38:22 -0500 Subject: [PATCH 2/2] add combssq to CMake --- SRC/CMakeLists.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index efe405e163..e6f1809327 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -150,7 +150,8 @@ set(SLASRC stplqt.f stplqt2.f stpmlqt.f ssytrd_2stage.f ssytrd_sy2sb.f ssytrd_sb2st.F ssb2st_kernels.f ssyevd_2stage.f ssyev_2stage.f ssyevx_2stage.f ssyevr_2stage.f - ssbev_2stage.f ssbevx_2stage.f ssbevd_2stage.f ssygv_2stage.f) + ssbev_2stage.f ssbevx_2stage.f ssbevd_2stage.f ssygv_2stage.f + scombssq.f) set(DSLASRC spotrs.f sgetrs.f spotrf.f sgetrf.f) @@ -341,7 +342,8 @@ set(DLASRC dtplqt.f dtplqt2.f dtpmlqt.f dsytrd_2stage.f dsytrd_sy2sb.f dsytrd_sb2st.F dsb2st_kernels.f dsyevd_2stage.f dsyev_2stage.f dsyevx_2stage.f dsyevr_2stage.f - dsbev_2stage.f dsbevx_2stage.f dsbevd_2stage.f dsygv_2stage.f) + dsbev_2stage.f dsbevx_2stage.f dsbevd_2stage.f dsygv_2stage.f + dcombssq.f) set(DXLASRC dgesvxx.f dgerfsx.f dla_gerfsx_extended.f dla_geamv.f dla_gercond.f dla_gerpvgrw.f dsysvxx.f dsyrfsx.f