--- rpl/lapack/lapack/dbdsqr.f 2012/12/14 14:22:27 1.11
+++ rpl/lapack/lapack/dbdsqr.f 2023/08/07 08:38:47 1.19
@@ -2,25 +2,25 @@
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DBDSQR + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download DBDSQR + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
* LDU, C, LDC, WORK, INFO )
-*
+*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
@@ -29,7 +29,7 @@
* DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
* $ VT( LDVT, * ), WORK( * )
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -40,9 +40,9 @@
*> left singular vectors from the singular value decomposition (SVD) of
*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
*> zero-shift QR algorithm. The SVD of B has the form
-*>
+*>
*> B = Q * S * P**T
-*>
+*>
*> where S is the diagonal matrix of singular values, Q is an orthogonal
*> matrix of left singular vectors, and P is an orthogonal matrix of
*> right singular vectors. If left singular vectors are requested, this
@@ -113,7 +113,7 @@
*> \verbatim
*> E is DOUBLE PRECISION array, dimension (N-1)
*> On entry, the N-1 offdiagonal elements of the bidiagonal
-*> matrix B.
+*> matrix B.
*> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
*> will contain the diagonal and superdiagonal elements of a
*> bidiagonal matrix orthogonally equivalent to the one given
@@ -166,7 +166,7 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is DOUBLE PRECISION array, dimension (4*N)
+*> WORK is DOUBLE PRECISION array, dimension (4*(N-1))
*> \endverbatim
*>
*> \param[out] INFO
@@ -179,7 +179,7 @@
*> = 1, a split was marked by a positive value in E
*> = 2, current block of Z not diagonalized after 30*N
*> iterations (in inner while loop)
-*> = 3, termination criterion of outer while loop not met
+*> = 3, termination criterion of outer while loop not met
*> (program created more than N unreduced blocks)
*> else NCVT = NRU = NCC = 0,
*> the algorithm did not converge; D and E contain the
@@ -212,17 +212,26 @@
*> algorithm through its inner loop. The algorithms stops
*> (and so fails to converge) if the number of passes
*> 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:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
-*
-*> \date November 2011
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
*> \ingroup auxOTHERcomputational
*
@@ -230,10 +239,9 @@
SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
$ LDU, C, LDC, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2011
*
* .. Scalar Arguments ..
CHARACTER UPLO
@@ -266,8 +274,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,
@@ -329,7 +337,7 @@
CALL DLASQ1( N, D, E, WORK, INFO )
*
* If INFO equals 2, dqds didn't finish, try to finish
-*
+*
IF( INFO .NE. 2 ) RETURN
INFO = 0
END IF
@@ -400,20 +408,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 +438,13 @@
*
IF( M.LE.1 )
$ GO TO 160
- IF( ITER.GT.MAXIT )
- $ GO TO 200
+*
+ 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
*