version 1.7, 2010/12/21 13:53:24
|
version 1.19, 2023/08/07 08:38:47
|
Line 1
|
Line 1
|
|
*> \brief \b DBDSQR |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DBDSQR + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f"> |
|
*> [TXT]</a> |
|
*> \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, |
SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, |
$ LDU, C, LDC, WORK, INFO ) |
$ LDU, C, LDC, WORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* January 2007 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER UPLO |
CHARACTER UPLO |
Line 15
|
Line 252
|
$ VT( LDVT, * ), WORK( * ) |
$ 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 .. |
* .. Parameters .. |
Line 170
|
Line 274
|
* .. |
* .. |
* .. 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, SMINL, SMINOA, |
$ SINR, SLL, SMAX, SMIN, SMINL, SMINOA, |
Line 231
|
Line 335
|
* |
* |
IF( .NOT.ROTATE ) THEN |
IF( .NOT.ROTATE ) THEN |
CALL DLASQ1( N, D, E, WORK, INFO ) |
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 |
END IF |
* |
* |
NM1 = N - 1 |
NM1 = N - 1 |
Line 300
|
Line 408
|
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 |
* |
* |
Line 329
|
Line 438
|
* |
* |
IF( M.LE.1 ) |
IF( M.LE.1 ) |
$ GO TO 160 |
$ 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 |
* Find diagonal block of matrix to work on |
* |
* |