From 7908efafc8c2c8f58bb52eabde879d117519ead3 Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Fri, 24 Mar 2017 22:05:22 +0100 Subject: [PATCH] Fix integer overflow in LAPACK DBDSQR, SBDSQR (#1135) * Fix integer overflow in DBDSQR As noted in lapack issue 135, an integer overflow in the calculation of the iteration limit could lead to an immediate return without any iterations having been performed if the input matrix is sufficiently big. * Fix integer overflow in SBDSQR As noted in lapack issue 135, an integer overflow in the calculation of the iteration limit could lead to an immediate return without any iterations having been performed if the input matrix is sufficiently big. * Fix integer overflow in threshold calculation Related to lapack issue 135, the threshold calculation can overflow as well as the multiplication is evaluated from left to right. Without explicit parentheses, the calculation would overflow for N >= 18919 * Fix integer overflow in threshold calculation Related to lapack issue 135, the threshold calculation can overflow as well as the multiplication is evaluated from left to right. Without explicit parentheses, the calculation would overflow for N >= 18919 --- lapack-netlib/SRC/dbdsqr.f | 32 ++++++++++++++++++++++++-------- lapack-netlib/SRC/sbdsqr.f | 30 +++++++++++++++++++++++------- 2 files changed, 47 insertions(+), 15 deletions(-) diff --git a/lapack-netlib/SRC/dbdsqr.f b/lapack-netlib/SRC/dbdsqr.f index b9894d897..c4cfbb3f7 100644 --- a/lapack-netlib/SRC/dbdsqr.f +++ b/lapack-netlib/SRC/dbdsqr.f @@ -214,6 +214,16 @@ *> through the inner loop exceeds MAXITR*N**2. *> \endverbatim * +*> \par Note: +* =========== +*> +*> \verbatim +*> Bug report from Cezary Dendek. +*> On March 23rd 2017, 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 +* * Authors: * ======== * @@ -266,8 +276,8 @@ * .. * .. Local Scalars .. LOGICAL LOWER, ROTATE - INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, - $ NM12, NM13, OLDLL, OLDM + INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M, + $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA, @@ -300,7 +310,7 @@ ELSE IF( NRU.LT.0 ) THEN INFO = -4 ELSE IF( NCC.LT.0 ) THEN - INFO = -5 + INFO = -5 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN INFO = -9 @@ -400,20 +410,21 @@ 40 CONTINUE 50 CONTINUE SMINOA = SMINOA / SQRT( DBLE( N ) ) - THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) + THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) ) ELSE * * Absolute accuracy desired * - THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) + THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) ) END IF * * Prepare for main iteration loop for the singular values * (MAXIT is the maximum number of passes through the inner * loop permitted before nonconvergence signalled.) * - MAXIT = MAXITR*N*N - ITER = 0 + MAXITDIVN = MAXITR*N + ITERDIVN = 0 + ITER = -1 OLDLL = -1 OLDM = -1 * @@ -429,8 +440,13 @@ * IF( M.LE.1 ) $ GO TO 160 - IF( ITER.GT.MAXIT ) +* + IF( ITER.GE.N ) THEN + ITER = ITER - N + ITERDIVN = ITERDIVN + 1 + IF (ITERDIVN.GE.MAXITDIVN ) $ GO TO 200 + END IF * * Find diagonal block of matrix to work on * diff --git a/lapack-netlib/SRC/sbdsqr.f b/lapack-netlib/SRC/sbdsqr.f index 7dc3dd7c9..e80ac4ea9 100644 --- a/lapack-netlib/SRC/sbdsqr.f +++ b/lapack-netlib/SRC/sbdsqr.f @@ -214,6 +214,16 @@ *> through the inner loop exceeds MAXITR*N**2. *> \endverbatim * +*> \par Note: +* =========== +*> +*> \verbatim +*> Bug report from Cezary Dendek. +*> On March 23rd 2017, 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 +* * Authors: * ======== * @@ -266,8 +276,8 @@ * .. * .. Local Scalars .. LOGICAL LOWER, ROTATE - INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, - $ NM12, NM13, OLDLL, OLDM + INTEGER I, IDIR, ISUB, ITER, ITERDIVN J, LL, LLL, M, + $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA, @@ -400,20 +410,21 @@ 40 CONTINUE 50 CONTINUE SMINOA = SMINOA / SQRT( REAL( N ) ) - THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) + THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) ) ELSE * * Absolute accuracy desired * - THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) + THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) ) END IF * * Prepare for main iteration loop for the singular values * (MAXIT is the maximum number of passes through the inner * loop permitted before nonconvergence signalled.) * - MAXIT = MAXITR*N*N - ITER = 0 + MAXITDIVN = MAXITR*N + ITERDIVN = 0 + ITER = -1 OLDLL = -1 OLDM = -1 * @@ -429,8 +440,13 @@ * IF( M.LE.1 ) $ GO TO 160 - IF( ITER.GT.MAXIT ) +* + IF( ITER.GE.N ) THEN + ITER = ITER - N + ITERDIVN = ITERDIVN + 1 + IF( ITERDIVN.GE.MAXITDIVN ) $ GO TO 200 + END IF * * Find diagonal block of matrix to work on *