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.
Reference-LAPACK · Nov 10, 2023 · 5b0bc5e · 5b0bc5e
1 parent 7eb4057
commit 5b0bc5e
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 4 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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 )
 *