version 1.3, 2010/08/06 15:28:50
|
version 1.9, 2011/11/21 22:19:44
|
Line 1
|
Line 1
|
|
*> \brief \b ZBDSQR |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZBDSQR + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbdsqr.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbdsqr.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbdsqr.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, |
|
* LDU, C, LDC, RWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER UPLO |
|
* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION D( * ), E( * ), RWORK( * ) |
|
* COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZBDSQR 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**H |
|
*> |
|
*> 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**H*VT instead of |
|
*> P**H, for given complex input matrices U and VT. When U and VT are |
|
*> the unitary matrices that reduce a general matrix A to bidiagonal |
|
*> form: A = U*B*VT, as computed by ZGEBRD, then |
|
*> |
|
*> A = (U*Q) * S * (P**H*VT) |
|
*> |
|
*> is the SVD of A. Optionally, the subroutine may also compute Q**H*C |
|
*> for a given complex 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 COMPLEX*16 array, dimension (LDVT, NCVT) |
|
*> On entry, an N-by-NCVT matrix VT. |
|
*> On exit, VT is overwritten by P**H * 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDC, NCC) |
|
*> On entry, an N-by-NCC matrix C. |
|
*> On exit, C is overwritten by Q**H * 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] RWORK |
|
*> \verbatim |
|
*> RWORK is DOUBLE PRECISION array, dimension (2*N) |
|
*> if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: If INFO = -i, the i-th argument had an illegal value |
|
*> > 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 |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2011 |
|
* |
|
*> \ingroup complex16OTHERcomputational |
|
* |
|
* ===================================================================== |
SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, |
SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, |
$ LDU, C, LDC, RWORK, INFO ) |
$ LDU, C, LDC, RWORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine (version 3.4.0) -- |
* -- 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..-- |
* November 2006 |
* November 2011 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER UPLO |
CHARACTER UPLO |
Line 15
|
Line 237
|
COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * ) |
COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZBDSQR 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**H |
|
* |
|
* 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**H*VT instead of |
|
* P**H, for given complex input matrices U and VT. When U and VT are |
|
* the unitary matrices that reduce a general matrix A to bidiagonal |
|
* form: A = U*B*VT, as computed by ZGEBRD, then |
|
* |
|
* A = (U*Q) * S * (P**H*VT) |
|
* |
|
* is the SVD of A. Optionally, the subroutine may also compute Q**H*C |
|
* for a given complex 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) COMPLEX*16 array, dimension (LDVT, NCVT) |
|
* On entry, an N-by-NCVT matrix VT. |
|
* On exit, VT is overwritten by P**H * 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) COMPLEX*16 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) COMPLEX*16 array, dimension (LDC, NCC) |
|
* On entry, an N-by-NCC matrix C. |
|
* On exit, C is overwritten by Q**H * 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. |
|
* |
|
* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) |
|
* if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: If INFO = -i, the i-th argument had an illegal value |
|
* > 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 224
|
Line 320
|
* |
* |
IF( .NOT.ROTATE ) THEN |
IF( .NOT.ROTATE ) THEN |
CALL DLASQ1( N, D, E, RWORK, INFO ) |
CALL DLASQ1( N, D, E, RWORK, 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 |