Skip to content

Scale the inputs of xGGLSE and xGGGLM into the safe range, as the other least-squares drivers do #1383

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rmlarsen wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Scale the inputs of xGGLSE and xGGGLM into the safe range, as the other least-squares drivers do #1383

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c7a7e38
Scale A, B and the right-hand sides of xGGLSE and xGGGLM independentl…
rmlarsen Sep 8, 2026
6b929b3
TESTING: Cover tiny inputs in xGGLSE and xGGGLM
rmlarsen Sep 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 76 additions & 5 deletions SRC/cggglm.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,8 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
* ===================================================================
*
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
COMPLEX CZERO, CONE
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
$ CONE = ( 1.0E+0, 0.0E+0 ) )
Expand All @@ -216,19 +218,24 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
LOGICAL LQUERY
INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
$ NB4, NP
INTEGER IA, IB, IBIG, ID, ISML, KA, KB, KD
REAL ANRM, BIGNUM, BNRM, DNRM, SMLNUM
* ..
* .. Local Arrays ..
REAL RWORK( 1 )
* ..
* .. External Subroutines ..
EXTERNAL CCOPY, CGEMV, CGGQRF, CTRTRS, CUNMQR,
$ CUNMRQ,
$ XERBLA
EXTERNAL CCOPY, CGEMV, CGGQRF, CLASCL, CSSCAL, CTRTRS,
$ CUNMQR, CUNMRQ, XERBLA
* ..
* .. External Functions ..
INTEGER ILAENV
REAL SROUNDUP_LWORK
EXTERNAL ILAENV, SROUNDUP_LWORK
REAL SLAMCH, CLANGE
EXTERNAL CLANGE, ILAENV, SLAMCH, SROUNDUP_LWORK
* ..
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN
INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -290,6 +297,62 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
RETURN
END IF
*
* Get machine parameters
*
SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
BIGNUM = ONE / SMLNUM
ISML = EXPONENT( SMLNUM )
IBIG = EXPONENT( BIGNUM ) - 1
*
* Scale A, B and d independently by powers of two 2**KA, 2**KB and
* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
* values whose EXPONENT lies in [ISML,IBIG]. x is then scaled by
* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is
* exact. A norm that is zero, infinite or NaN takes no part.
* CLASCL is called with both endpoints at or above one so that
* forming the factor raises no underflow.
*
ANRM = CLANGE( 'M', N, M, A, LDA, RWORK )
KA = 0
IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
IA = EXPONENT( ANRM )
IF( IA.LT.ISML ) THEN
KA = ISML - IA
ELSE IF( IA.GT.IBIG ) THEN
KA = IBIG - IA
END IF
END IF
IF( KA.NE.0 )
$ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
$ SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
*
BNRM = CLANGE( 'M', N, P, B, LDB, RWORK )
KB = 0
IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
IB = EXPONENT( BNRM )
IF( IB.LT.ISML ) THEN
KB = ISML - IB
ELSE IF( IB.GT.IBIG ) THEN
KB = IBIG - IB
END IF
END IF
IF( KB.NE.0 )
$ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
$ SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
*
DNRM = CLANGE( 'M', N, 1, D, N, RWORK )
KD = 0
IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
ID = EXPONENT( DNRM )
IF( ID.LT.ISML ) THEN
KD = ISML - ID
ELSE IF( ID.GT.IBIG ) THEN
KD = IBIG - ID
END IF
END IF
IF( KD.NE.0 )
$ CALL CSSCAL( N, SCALE( ONE, KD ), D, 1 )
*
* Compute the GQR factorization of matrices A and B:
*
* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
Expand Down Expand Up @@ -359,6 +422,14 @@ SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
CALL CUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
$ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
$ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
*
* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
*
IF( KA.NE.KD )
$ CALL CSSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
IF( KB.NE.KD )
$ CALL CSSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
*
WORK( 1 ) = CMPLX( M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) )
*
RETURN
Expand Down
92 changes: 87 additions & 5 deletions SRC/cgglse.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,26 +203,33 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
COMPLEX CONE
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
LOGICAL LQUERY
INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
$ NB4, NR
INTEGER IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT
REAL ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM
* ..
* .. Local Arrays ..
REAL RWORK( 1 )
* ..
* .. External Subroutines ..
EXTERNAL CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV,
$ CTRTRS,
$ CUNMQR, CUNMRQ, XERBLA
EXTERNAL CAXPY, CCOPY, CGEMV, CGGRQF, CLASCL, CSSCAL,
$ CTRMV, CTRTRS, CUNMQR, CUNMRQ, XERBLA
* ..
* .. External Functions ..
INTEGER ILAENV
REAL SROUNDUP_LWORK
EXTERNAL ILAENV, SROUNDUP_LWORK
REAL SLAMCH, CLANGE
EXTERNAL CLANGE, ILAENV, SLAMCH, SROUNDUP_LWORK
* ..
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN
INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -277,6 +284,72 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
IF( N.EQ.0 )
$ RETURN
*
* Get machine parameters
*
SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
BIGNUM = ONE / SMLNUM
ISML = EXPONENT( SMLNUM )
IBIG = EXPONENT( BIGNUM ) - 1
*
* Scale A, B, c and d by powers of two so that their largest
* entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies
* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and
* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which
* keeps the constraint and the residual consistent and scales x by
* 2**KT; KT brings the larger of the two right-hand sides into
* range. Scaling by a power of two is exact. A norm that is zero,
* infinite or NaN takes no part. CLASCL is called with both
* endpoints at or above one so that forming the factor raises no
* underflow.
*
ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
KA = 0
IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
IA = EXPONENT( ANRM )
IF( IA.LT.ISML ) THEN
KA = ISML - IA
ELSE IF( IA.GT.IBIG ) THEN
KA = IBIG - IA
END IF
END IF
IF( KA.NE.0 )
$ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
$ SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO )
*
BNRM = CLANGE( 'M', P, N, B, LDB, RWORK )
KB = 0
IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
IB = EXPONENT( BNRM )
IF( IB.LT.ISML ) THEN
KB = ISML - IB
ELSE IF( IB.GT.IBIG ) THEN
KB = IBIG - IB
END IF
END IF
IF( KB.NE.0 )
$ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
$ SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO )
*
CNRM = CLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK )
DNRM = CLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK )
IC = -HUGE( 0 )
IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) )
$ IC = EXPONENT( CNRM ) + KA
ID = -HUGE( 0 )
IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) )
$ ID = EXPONENT( DNRM ) + KB
IU = MAX( IC, ID )
KT = 0
IF( IU.GT.IBIG ) THEN
KT = IBIG - IU
ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN
KT = ISML - IU
END IF
IF( KA+KT.NE.0 )
$ CALL CSSCAL( M, SCALE( ONE, KA+KT ), C, 1 )
IF( KB+KT.NE.0 )
$ CALL CSSCAL( P, SCALE( ONE, KB+KT ), D, 1 )
*
* Compute the GRQ factorization of matrices B and A:
*
* B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P
Expand Down Expand Up @@ -357,6 +430,15 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK,
*
CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
$ WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
*
* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M)
* carries 2**(KA+KT)
*
IF( KT.NE.0 )
$ CALL CSSCAL( N, SCALE( ONE, -KT ), X, 1 )
IF( KA+KT.NE.0 .AND. M+P.GT.N )
$ CALL CSSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 )
*
WORK( 1 ) = CMPLX( P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) )
*
RETURN
Expand Down
79 changes: 74 additions & 5 deletions SRC/dggglm.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,18 +215,23 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
LOGICAL LQUERY
INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
$ NB4, NP
INTEGER IA, IB, IBIG, ID, ISML, KA, KB, KD
DOUBLE PRECISION ANRM, BIGNUM, BNRM, DNRM, SMLNUM
* ..
* .. Local Arrays ..
DOUBLE PRECISION RWORK( 1 )
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ,
$ DTRTRS,
$ XERBLA
EXTERNAL DCOPY, DGEMV, DGGQRF, DLASCL, DORMQR, DORMRQ,
$ DSCAL, DTRTRS, XERBLA
* ..
* .. External Functions ..
INTEGER ILAENV
EXTERNAL ILAENV
DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE, ILAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN
INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -288,6 +293,62 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
RETURN
END IF
*
* Get machine parameters
*
SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
BIGNUM = ONE / SMLNUM
ISML = EXPONENT( SMLNUM )
IBIG = EXPONENT( BIGNUM ) - 1
*
* Scale A, B and d independently by powers of two 2**KA, 2**KB and
* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the
* values whose EXPONENT lies in [ISML,IBIG]. x is then scaled by
* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is
* exact. A norm that is zero, infinite or NaN takes no part.
* DLASCL is called with both endpoints at or above one so that
* forming the factor raises no underflow.
*
ANRM = DLANGE( 'M', N, M, A, LDA, RWORK )
KA = 0
IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN
IA = EXPONENT( ANRM )
IF( IA.LT.ISML ) THEN
KA = ISML - IA
ELSE IF( IA.GT.IBIG ) THEN
KA = IBIG - IA
END IF
END IF
IF( KA.NE.0 )
$ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ),
$ SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO )
*
BNRM = DLANGE( 'M', N, P, B, LDB, RWORK )
KB = 0
IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN
IB = EXPONENT( BNRM )
IF( IB.LT.ISML ) THEN
KB = ISML - IB
ELSE IF( IB.GT.IBIG ) THEN
KB = IBIG - IB
END IF
END IF
IF( KB.NE.0 )
$ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ),
$ SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO )
*
DNRM = DLANGE( 'M', N, 1, D, N, RWORK )
KD = 0
IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN
ID = EXPONENT( DNRM )
IF( ID.LT.ISML ) THEN
KD = ISML - ID
ELSE IF( ID.GT.IBIG ) THEN
KD = IBIG - ID
END IF
END IF
IF( KD.NE.0 )
$ CALL DSCAL( N, SCALE( ONE, KD ), D, 1 )
*
* Compute the GQR factorization of matrices A and B:
*
* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M
Expand Down Expand Up @@ -355,6 +416,14 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK,
CALL DORMRQ( 'Left', 'Transpose', P, 1, NP,
$ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
$ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
*
* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB)
*
IF( KA.NE.KD )
$ CALL DSCAL( M, SCALE( ONE, KA-KD ), X, 1 )
IF( KB.NE.KD )
$ CALL DSCAL( P, SCALE( ONE, KB-KD ), Y, 1 )
*
WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
*
RETURN
Expand Down
Loading
Loading