--- rpl/lapack/lapack/dbdsqr.f 2010/08/07 13:22:11 1.5
+++ rpl/lapack/lapack/dbdsqr.f 2023/08/07 08:38:47 1.19
@@ -1,10 +1,247 @@
+*> \brief \b DBDSQR
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DBDSQR + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \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
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
+* $ VT( LDVT, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DBDSQR computes the singular values and, optionally, the right and/or
+*> 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
+*> subroutine actually returns U*Q instead of Q, and, if right singular
+*> vectors are requested, this subroutine returns P**T*VT instead of
+*> P**T, for given real input matrices U and VT. When U and VT are the
+*> orthogonal matrices that reduce a general matrix A to bidiagonal
+*> form: A = U*B*VT, as computed by DGEBRD, then
+*>
+*> A = (U*Q) * S * (P**T*VT)
+*>
+*> is the SVD of A. Optionally, the subroutine may also compute Q**T*C
+*> for a given real input matrix C.
+*>
+*> See "Computing Small Singular Values of Bidiagonal Matrices With
+*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
+*> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
+*> no. 5, pp. 873-912, Sept 1990) and
+*> "Accurate singular values and differential qd algorithms," by
+*> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
+*> Department, University of California at Berkeley, July 1992
+*> for a detailed description of the algorithm.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': B is upper bidiagonal;
+*> = 'L': B is lower bidiagonal.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NCVT
+*> \verbatim
+*> NCVT is INTEGER
+*> The number of columns of the matrix VT. NCVT >= 0.
+*> \endverbatim
+*>
+*> \param[in] NRU
+*> \verbatim
+*> NRU is INTEGER
+*> The number of rows of the matrix U. NRU >= 0.
+*> \endverbatim
+*>
+*> \param[in] NCC
+*> \verbatim
+*> NCC is INTEGER
+*> The number of columns of the matrix C. NCC >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] D
+*> \verbatim
+*> D is DOUBLE PRECISION array, dimension (N)
+*> On entry, the n diagonal elements of the bidiagonal matrix B.
+*> On exit, if INFO=0, the singular values of B in decreasing
+*> order.
+*> \endverbatim
+*>
+*> \param[in,out] E
+*> \verbatim
+*> E is DOUBLE PRECISION array, dimension (N-1)
+*> On entry, the N-1 offdiagonal elements of the bidiagonal
+*> 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
+*> as input.
+*> \endverbatim
+*>
+*> \param[in,out] VT
+*> \verbatim
+*> VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
+*> On entry, an N-by-NCVT matrix VT.
+*> On exit, VT is overwritten by P**T * VT.
+*> Not referenced if NCVT = 0.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*> LDVT is INTEGER
+*> The leading dimension of the array VT.
+*> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
+*> \endverbatim
+*>
+*> \param[in,out] U
+*> \verbatim
+*> U is DOUBLE PRECISION array, dimension (LDU, N)
+*> On entry, an NRU-by-N matrix U.
+*> On exit, U is overwritten by U * Q.
+*> Not referenced if NRU = 0.
+*> \endverbatim
+*>
+*> \param[in] LDU
+*> \verbatim
+*> LDU is INTEGER
+*> The leading dimension of the array U. LDU >= max(1,NRU).
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*> C is DOUBLE PRECISION array, dimension (LDC, NCC)
+*> On entry, an N-by-NCC matrix C.
+*> On exit, C is overwritten by Q**T * C.
+*> Not referenced if NCC = 0.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*> LDC is INTEGER
+*> The leading dimension of the array C.
+*> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (4*(N-1))
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: If INFO = -i, the i-th argument had an illegal value
+*> > 0:
+*> if NCVT = NRU = NCC = 0,
+*> = 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
+*> (program created more than N unreduced blocks)
+*> else NCVT = NRU = NCC = 0,
+*> the algorithm did not converge; D and E contain the
+*> elements of a bidiagonal matrix which is orthogonally
+*> similar to the input matrix B; if INFO = i, i
+*> elements of E have not converged to zero.
+*> \endverbatim
+*
+*> \par Internal Parameters:
+* =========================
+*>
+*> \verbatim
+*> TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
+*> TOLMUL controls the convergence criterion of the QR loop.
+*> If it is positive, TOLMUL*EPS is the desired relative
+*> precision in the computed singular values.
+*> If it is negative, abs(TOLMUL*EPS*sigma_max) is the
+*> desired absolute accuracy in the computed singular
+*> values (corresponds to relative accuracy
+*> abs(TOLMUL*EPS) in the largest singular value.
+*> abs(TOLMUL) should be between 1 and 1/EPS, and preferably
+*> between 10 (for fast convergence) and .1/EPS
+*> (for there to be some accuracy in the results).
+*> Default is to lose at either one eighth or 2 of the
+*> available decimal digits in each computed singular value
+*> (whichever is smaller).
+*>
+*> MAXITR INTEGER, default = 6
+*> MAXITR controls the maximum number of passes of the
+*> 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.
+*
+*> \ingroup auxOTHERcomputational
+*
+* =====================================================================
SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
$ LDU, C, LDC, WORK, INFO )
*
-* -- LAPACK routine (version 3.2) --
+* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* January 2007
*
* .. Scalar Arguments ..
CHARACTER UPLO
@@ -15,139 +252,6 @@
$ VT( LDVT, * ), WORK( * )
* ..
*
-* Purpose
-* =======
-*
-* DBDSQR computes the singular values and, optionally, the right and/or
-* 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
-* subroutine actually returns U*Q instead of Q, and, if right singular
-* vectors are requested, this subroutine returns P**T*VT instead of
-* P**T, for given real input matrices U and VT. When U and VT are the
-* orthogonal matrices that reduce a general matrix A to bidiagonal
-* form: A = U*B*VT, as computed by DGEBRD, then
-*
-* A = (U*Q) * S * (P**T*VT)
-*
-* is the SVD of A. Optionally, the subroutine may also compute Q**T*C
-* for a given real input matrix C.
-*
-* See "Computing Small Singular Values of Bidiagonal Matrices With
-* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
-* LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
-* no. 5, pp. 873-912, Sept 1990) and
-* "Accurate singular values and differential qd algorithms," by
-* B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
-* Department, University of California at Berkeley, July 1992
-* for a detailed description of the algorithm.
-*
-* Arguments
-* =========
-*
-* UPLO (input) CHARACTER*1
-* = 'U': B is upper bidiagonal;
-* = 'L': B is lower bidiagonal.
-*
-* N (input) INTEGER
-* The order of the matrix B. N >= 0.
-*
-* NCVT (input) INTEGER
-* The number of columns of the matrix VT. NCVT >= 0.
-*
-* NRU (input) INTEGER
-* The number of rows of the matrix U. NRU >= 0.
-*
-* NCC (input) INTEGER
-* The number of columns of the matrix C. NCC >= 0.
-*
-* D (input/output) DOUBLE PRECISION array, dimension (N)
-* On entry, the n diagonal elements of the bidiagonal matrix B.
-* On exit, if INFO=0, the singular values of B in decreasing
-* order.
-*
-* E (input/output) DOUBLE PRECISION array, dimension (N-1)
-* On entry, the N-1 offdiagonal elements of the bidiagonal
-* 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
-* as input.
-*
-* VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
-* On entry, an N-by-NCVT matrix VT.
-* On exit, VT is overwritten by P**T * VT.
-* Not referenced if NCVT = 0.
-*
-* LDVT (input) INTEGER
-* The leading dimension of the array VT.
-* LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
-*
-* U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
-* On entry, an NRU-by-N matrix U.
-* On exit, U is overwritten by U * Q.
-* Not referenced if NRU = 0.
-*
-* LDU (input) INTEGER
-* The leading dimension of the array U. LDU >= max(1,NRU).
-*
-* C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
-* On entry, an N-by-NCC matrix C.
-* On exit, C is overwritten by Q**T * C.
-* Not referenced if NCC = 0.
-*
-* LDC (input) INTEGER
-* The leading dimension of the array C.
-* LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
-*
-* WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: If INFO = -i, the i-th argument had an illegal value
-* > 0:
-* if NCVT = NRU = NCC = 0,
-* = 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
-* (program created more than N unreduced blocks)
-* else NCVT = NRU = NCC = 0,
-* the algorithm did not converge; D and E contain the
-* elements of a bidiagonal matrix which is orthogonally
-* similar to the input matrix B; if INFO = i, i
-* elements of E have not converged to zero.
-*
-* Internal Parameters
-* ===================
-*
-* TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
-* TOLMUL controls the convergence criterion of the QR loop.
-* If it is positive, TOLMUL*EPS is the desired relative
-* precision in the computed singular values.
-* If it is negative, abs(TOLMUL*EPS*sigma_max) is the
-* desired absolute accuracy in the computed singular
-* values (corresponds to relative accuracy
-* abs(TOLMUL*EPS) in the largest singular value.
-* abs(TOLMUL) should be between 1 and 1/EPS, and preferably
-* between 10 (for fast convergence) and .1/EPS
-* (for there to be some accuracy in the results).
-* Default is to lose at either one eighth or 2 of the
-* available decimal digits in each computed singular value
-* (whichever is smaller).
-*
-* MAXITR INTEGER, default = 6
-* MAXITR controls the maximum number of passes of the
-* 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.
-*
* =====================================================================
*
* .. Parameters ..
@@ -170,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,
@@ -231,7 +335,11 @@
*
IF( .NOT.ROTATE ) THEN
CALL DLASQ1( N, D, E, WORK, INFO )
- RETURN
+*
+* If INFO equals 2, dqds didn't finish, try to finish
+*
+ IF( INFO .NE. 2 ) RETURN
+ INFO = 0
END IF
*
NM1 = N - 1
@@ -300,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
*
@@ -329,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
*