From 5b0bc5e555a74c35782a7665bec212172bd42cb3 Mon Sep 17 00:00:00 2001 From: Christoph Conrads Date: Fri, 10 Nov 2023 16:29:41 +0100 Subject: [PATCH] xORBDB5/xUNBDB5: ensure xORBDB6 input has unit norm Call sites may run into problems when this subroutine computes nonzero vectors that are very small in norm. --- SRC/cunbdb5.f | 9 ++++++++- SRC/dorbdb5.f | 9 ++++++++- SRC/sorbdb5.f | 9 ++++++++- SRC/zunbdb5.f | 9 ++++++++- 4 files changed, 32 insertions(+), 4 deletions(-) diff --git a/SRC/cunbdb5.f b/SRC/cunbdb5.f index 184e445c8e..22513cf8b1 100644 --- a/SRC/cunbdb5.f +++ b/SRC/cunbdb5.f @@ -179,7 +179,7 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, REAL EPS, NORM, SCL, SSQ * .. * .. External Subroutines .. - EXTERNAL CLASSQ, CUNBDB6, XERBLA + EXTERNAL CLASSQ, CUNBDB6, CSCAL, XERBLA * .. * .. External Functions .. REAL SLAMCH, SCNRM2 @@ -227,6 +227,13 @@ SUBROUTINE CUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, NORM = SCL * SQRT( SSQ ) * IF( NORM .GT. N * EPS ) THEN +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL CSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL CSCAL( M2, ONE / NORM, X2, INCX2 ) CALL CUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) * diff --git a/SRC/dorbdb5.f b/SRC/dorbdb5.f index 5a52107b77..cbd58ae547 100644 --- a/SRC/dorbdb5.f +++ b/SRC/dorbdb5.f @@ -179,7 +179,7 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, DOUBLE PRECISION EPS, NORM, SCL, SSQ * .. * .. External Subroutines .. - EXTERNAL DLASSQ, DORBDB6, XERBLA + EXTERNAL DLASSQ, DORBDB6, DSCAL, XERBLA * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DNRM2 @@ -227,6 +227,13 @@ SUBROUTINE DORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, NORM = SCL * SQRT( SSQ ) * IF( NORM .GT. N * EPS ) THEN +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL DSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL DSCAL( M2, ONE / NORM, X2, INCX2 ) CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) * diff --git a/SRC/sorbdb5.f b/SRC/sorbdb5.f index 6aa166e7b1..8fb88876fe 100644 --- a/SRC/sorbdb5.f +++ b/SRC/sorbdb5.f @@ -179,7 +179,7 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, REAL EPS, NORM, SCL, SSQ * .. * .. External Subroutines .. - EXTERNAL SLASSQ, SORBDB6, XERBLA + EXTERNAL SLASSQ, SORBDB6, SSCAL, XERBLA * .. * .. External Functions .. REAL SLAMCH, SNRM2 @@ -227,6 +227,13 @@ SUBROUTINE SORBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, NORM = SCL * SQRT( SSQ ) * IF( NORM .GT. N * EPS ) THEN +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL SSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL SSCAL( M2, ONE / NORM, X2, INCX2 ) CALL SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) * diff --git a/SRC/zunbdb5.f b/SRC/zunbdb5.f index b867ec14b4..c451ae921a 100644 --- a/SRC/zunbdb5.f +++ b/SRC/zunbdb5.f @@ -179,7 +179,7 @@ SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, DOUBLE PRECISION EPS, NORM, SCL, SSQ * .. * .. External Subroutines .. - EXTERNAL ZLASSQ, ZUNBDB6, XERBLA + EXTERNAL ZLASSQ, ZUNBDB6, ZSCAL, XERBLA * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DZNRM2 @@ -227,6 +227,13 @@ SUBROUTINE ZUNBDB5( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, NORM = SCL * SQRT( SSQ ) * IF( NORM .GT. N * EPS ) THEN +* Scale vector to unit norm to avoid problems in the caller code. +* Computing the reciprocal is undesirable but +* * xLASCL cannot be used because of the vector increments and +* * the round-off error has a negligible impact on +* orthogonalization. + CALL ZSCAL( M1, ONE / NORM, X1, INCX1 ) + CALL ZSCAL( M2, ONE / NORM, X2, INCX2 ) CALL ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, $ LDQ2, WORK, LWORK, CHILDINFO ) *