--- rpl/lapack/lapack/zbdsqr.f 2010/12/21 13:53:42 1.7
+++ rpl/lapack/lapack/zbdsqr.f 2011/11/21 20:43:07 1.8
@@ -1,10 +1,232 @@
+*> \brief \b ZBDSQR
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZBDSQR + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \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,
$ 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, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
+* November 2011
*
* .. Scalar Arguments ..
CHARACTER UPLO
@@ -15,132 +237,6 @@
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 ..
@@ -224,7 +320,11 @@
*
IF( .NOT.ROTATE ) THEN
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
*
NM1 = N - 1