Fix potential overflow in the calculation of MAXIT

This commit is contained in:
Martin Kroeker 2023-11-06 21:22:26 +01:00 committed by GitHub
parent cd8eb83bae
commit cf8295da5c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 50 additions and 18 deletions

View File

@ -204,6 +204,17 @@
*> algorithm through its inner loop. The algorithms stops *> algorithm through its inner loop. The algorithms stops
*> (and so fails to converge) if the number of passes *> (and so fails to converge) if the number of passes
*> through the inner loop exceeds MAXITR*N**2. *> through the inner loop exceeds MAXITR*N**2.
*>
*> \endverbatim
*
*> \par Note:
* ===========
*>
*> \verbatim
*> Bug report from Cezary Dendek.
*> On November 3rd 2023, the INTEGER variable MAXIT = MAXITR*N**2 is
*> removed since it can overflow pretty easily (for N larger or equal
*> than 18,919). We instead use MAXITDIVN = MAXITR*N.
*> \endverbatim *> \endverbatim
* *
* Authors: * Authors:
@ -214,7 +225,7 @@
*> \author Univ. of Colorado Denver *> \author Univ. of Colorado Denver
*> \author NAG Ltd. *> \author NAG Ltd.
* *
*> \ingroup complexOTHERcomputational *> \ingroup bdsqr
* *
* ===================================================================== * =====================================================================
SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
@ -255,8 +266,8 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL LOWER, ROTATE LOGICAL LOWER, ROTATE
INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
$ NM12, NM13, OLDLL, OLDM $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
$ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
$ SINR, SLL, SMAX, SMIN, SMINOA, $ SINR, SLL, SMAX, SMIN, SMINOA,
@ -389,20 +400,21 @@
40 CONTINUE 40 CONTINUE
50 CONTINUE 50 CONTINUE
SMINOA = SMINOA / SQRT( REAL( N ) ) SMINOA = SMINOA / SQRT( REAL( N ) )
THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) )
ELSE ELSE
* *
* Absolute accuracy desired * Absolute accuracy desired
* *
THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) )
END IF END IF
* *
* Prepare for main iteration loop for the singular values * Prepare for main iteration loop for the singular values
* (MAXIT is the maximum number of passes through the inner * (MAXIT is the maximum number of passes through the inner
* loop permitted before nonconvergence signalled.) * loop permitted before nonconvergence signalled.)
* *
MAXIT = MAXITR*N*N MAXITDIVN = MAXITR*N
ITER = 0 ITERDIVN = 0
ITER = -1
OLDLL = -1 OLDLL = -1
OLDM = -1 OLDM = -1
* *
@ -418,8 +430,12 @@
* *
IF( M.LE.1 ) IF( M.LE.1 )
$ GO TO 160 $ GO TO 160
IF( ITER.GT.MAXIT ) IF( ITER.GE.N ) THEN
$ GO TO 200 ITER = ITER - N
ITERDIVN = ITERDIVN + 1
IF( ITERDIVN.GE.MAXITDIVN )
$ GO TO 200
END IF
* *
* Find diagonal block of matrix to work on * Find diagonal block of matrix to work on
* *

View File

@ -204,6 +204,17 @@
*> algorithm through its inner loop. The algorithms stops *> algorithm through its inner loop. The algorithms stops
*> (and so fails to converge) if the number of passes *> (and so fails to converge) if the number of passes
*> through the inner loop exceeds MAXITR*N**2. *> through the inner loop exceeds MAXITR*N**2.
*>
*> \endverbatim
*
*> \par Note:
* ===========
*>
*> \verbatim
*> Bug report from Cezary Dendek.
*> On November 3rd 2023, the INTEGER variable MAXIT = MAXITR*N**2 is
*> removed since it can overflow pretty easily (for N larger or equal
*> than 18,919). We instead use MAXITDIVN = MAXITR*N.
*> \endverbatim *> \endverbatim
* *
* Authors: * Authors:
@ -214,7 +225,7 @@
*> \author Univ. of Colorado Denver *> \author Univ. of Colorado Denver
*> \author NAG Ltd. *> \author NAG Ltd.
* *
*> \ingroup complex16OTHERcomputational *> \ingroup bdsqr
* *
* ===================================================================== * =====================================================================
SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
@ -255,8 +266,8 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL LOWER, ROTATE LOGICAL LOWER, ROTATE
INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
$ NM12, NM13, OLDLL, OLDM $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
$ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
$ SINR, SLL, SMAX, SMIN, SMINOA, $ SINR, SLL, SMAX, SMIN, SMINOA,
@ -389,20 +400,21 @@
40 CONTINUE 40 CONTINUE
50 CONTINUE 50 CONTINUE
SMINOA = SMINOA / SQRT( DBLE( N ) ) SMINOA = SMINOA / SQRT( DBLE( N ) )
THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) )
ELSE ELSE
* *
* Absolute accuracy desired * Absolute accuracy desired
* *
THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) )
END IF END IF
* *
* Prepare for main iteration loop for the singular values * Prepare for main iteration loop for the singular values
* (MAXIT is the maximum number of passes through the inner * (MAXIT is the maximum number of passes through the inner
* loop permitted before nonconvergence signalled.) * loop permitted before nonconvergence signalled.)
* *
MAXIT = MAXITR*N*N MAXITDIVN = MAXITR*N
ITER = 0 ITERDIVN = 0
ITER = -1
OLDLL = -1 OLDLL = -1
OLDM = -1 OLDM = -1
* *
@ -418,8 +430,12 @@
* *
IF( M.LE.1 ) IF( M.LE.1 )
$ GO TO 160 $ GO TO 160
IF( ITER.GT.MAXIT ) IF( ITER.GE.N ) THEN
$ GO TO 200 ITER = ITER - N
ITERDIVN = ITERDIVN + 1
IF( ITERDIVN.GE.MAXITDIVN )
$ GO TO 200
END IF
* *
* Find diagonal block of matrix to work on * Find diagonal block of matrix to work on
* *