version 1.9, 2011/07/22 07:38:04
|
version 1.21, 2023/08/07 08:38:47
|
Line 1
|
Line 1
|
|
*> \brief \b DBDSDC |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DBDSDC + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, |
|
* WORK, IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER COMPQ, UPLO |
|
* INTEGER INFO, LDU, LDVT, N |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IQ( * ), IWORK( * ) |
|
* DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ), |
|
* $ VT( LDVT, * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DBDSDC computes the singular value decomposition (SVD) of a real |
|
*> N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, |
|
*> using a divide and conquer method, where S is a diagonal matrix |
|
*> with non-negative diagonal elements (the singular values of B), and |
|
*> U and VT are orthogonal matrices of left and right singular vectors, |
|
*> respectively. DBDSDC can be used to compute all singular values, |
|
*> and optionally, singular vectors or singular vectors in compact form. |
|
*> |
|
*> This code makes very mild assumptions about floating point |
|
*> arithmetic. It will work on machines with a guard digit in |
|
*> add/subtract, or on those binary machines without guard digits |
|
*> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. |
|
*> It could conceivably fail on hexadecimal or decimal machines |
|
*> without guard digits, but we know of none. See DLASD3 for details. |
|
*> |
|
*> The code currently calls DLASDQ if singular values only are desired. |
|
*> However, it can be slightly modified to compute singular values |
|
*> using the divide and conquer method. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] UPLO |
|
*> \verbatim |
|
*> UPLO is CHARACTER*1 |
|
*> = 'U': B is upper bidiagonal. |
|
*> = 'L': B is lower bidiagonal. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] COMPQ |
|
*> \verbatim |
|
*> COMPQ is CHARACTER*1 |
|
*> Specifies whether singular vectors are to be computed |
|
*> as follows: |
|
*> = 'N': Compute singular values only; |
|
*> = 'P': Compute singular values and compute singular |
|
*> vectors in compact form; |
|
*> = 'I': Compute singular values and singular vectors. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The order of the matrix B. N >= 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. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] E |
|
*> \verbatim |
|
*> E is DOUBLE PRECISION array, dimension (N-1) |
|
*> On entry, the elements of E contain the offdiagonal |
|
*> elements of the bidiagonal matrix whose SVD is desired. |
|
*> On exit, E has been destroyed. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] U |
|
*> \verbatim |
|
*> U is DOUBLE PRECISION array, dimension (LDU,N) |
|
*> If COMPQ = 'I', then: |
|
*> On exit, if INFO = 0, U contains the left singular vectors |
|
*> of the bidiagonal matrix. |
|
*> For other values of COMPQ, U is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDU |
|
*> \verbatim |
|
*> LDU is INTEGER |
|
*> The leading dimension of the array U. LDU >= 1. |
|
*> If singular vectors are desired, then LDU >= max( 1, N ). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] VT |
|
*> \verbatim |
|
*> VT is DOUBLE PRECISION array, dimension (LDVT,N) |
|
*> If COMPQ = 'I', then: |
|
*> On exit, if INFO = 0, VT**T contains the right singular |
|
*> vectors of the bidiagonal matrix. |
|
*> For other values of COMPQ, VT is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVT |
|
*> \verbatim |
|
*> LDVT is INTEGER |
|
*> The leading dimension of the array VT. LDVT >= 1. |
|
*> If singular vectors are desired, then LDVT >= max( 1, N ). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] Q |
|
*> \verbatim |
|
*> Q is DOUBLE PRECISION array, dimension (LDQ) |
|
*> If COMPQ = 'P', then: |
|
*> On exit, if INFO = 0, Q and IQ contain the left |
|
*> and right singular vectors in a compact form, |
|
*> requiring O(N log N) space instead of 2*N**2. |
|
*> In particular, Q contains all the DOUBLE PRECISION data in |
|
*> LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) |
|
*> words of memory, where SMLSIZ is returned by ILAENV and |
|
*> is equal to the maximum size of the subproblems at the |
|
*> bottom of the computation tree (usually about 25). |
|
*> For other values of COMPQ, Q is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IQ |
|
*> \verbatim |
|
*> IQ is INTEGER array, dimension (LDIQ) |
|
*> If COMPQ = 'P', then: |
|
*> On exit, if INFO = 0, Q and IQ contain the left |
|
*> and right singular vectors in a compact form, |
|
*> requiring O(N log N) space instead of 2*N**2. |
|
*> In particular, IQ contains all INTEGER data in |
|
*> LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) |
|
*> words of memory, where SMLSIZ is returned by ILAENV and |
|
*> is equal to the maximum size of the subproblems at the |
|
*> bottom of the computation tree (usually about 25). |
|
*> For other values of COMPQ, IQ is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
*> If COMPQ = 'N' then LWORK >= (4 * N). |
|
*> If COMPQ = 'P' then LWORK >= (6 * N). |
|
*> If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (8*N) |
|
*> \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 failed to compute a singular value. |
|
*> The update process of divide and conquer failed. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \ingroup auxOTHERcomputational |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Ming Gu and Huan Ren, Computer Science Division, University of |
|
*> California at Berkeley, USA |
|
*> |
|
* ===================================================================== |
SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, |
SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, |
$ WORK, IWORK, INFO ) |
$ WORK, IWORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.3.1) -- |
* -- 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..-- |
* -- April 2011 -- |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER COMPQ, UPLO |
CHARACTER COMPQ, UPLO |
Line 16
|
Line 217
|
$ VT( LDVT, * ), WORK( * ) |
$ VT( LDVT, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DBDSDC computes the singular value decomposition (SVD) of a real |
|
* N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, |
|
* using a divide and conquer method, where S is a diagonal matrix |
|
* with non-negative diagonal elements (the singular values of B), and |
|
* U and VT are orthogonal matrices of left and right singular vectors, |
|
* respectively. DBDSDC can be used to compute all singular values, |
|
* and optionally, singular vectors or singular vectors in compact form. |
|
* |
|
* This code makes very mild assumptions about floating point |
|
* arithmetic. It will work on machines with a guard digit in |
|
* add/subtract, or on those binary machines without guard digits |
|
* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. |
|
* It could conceivably fail on hexadecimal or decimal machines |
|
* without guard digits, but we know of none. See DLASD3 for details. |
|
* |
|
* The code currently calls DLASDQ if singular values only are desired. |
|
* However, it can be slightly modified to compute singular values |
|
* using the divide and conquer method. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* UPLO (input) CHARACTER*1 |
|
* = 'U': B is upper bidiagonal. |
|
* = 'L': B is lower bidiagonal. |
|
* |
|
* COMPQ (input) CHARACTER*1 |
|
* Specifies whether singular vectors are to be computed |
|
* as follows: |
|
* = 'N': Compute singular values only; |
|
* = 'P': Compute singular values and compute singular |
|
* vectors in compact form; |
|
* = 'I': Compute singular values and singular vectors. |
|
* |
|
* N (input) INTEGER |
|
* The order of the matrix B. N >= 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. |
|
* |
|
* E (input/output) DOUBLE PRECISION array, dimension (N-1) |
|
* On entry, the elements of E contain the offdiagonal |
|
* elements of the bidiagonal matrix whose SVD is desired. |
|
* On exit, E has been destroyed. |
|
* |
|
* U (output) DOUBLE PRECISION array, dimension (LDU,N) |
|
* If COMPQ = 'I', then: |
|
* On exit, if INFO = 0, U contains the left singular vectors |
|
* of the bidiagonal matrix. |
|
* For other values of COMPQ, U is not referenced. |
|
* |
|
* LDU (input) INTEGER |
|
* The leading dimension of the array U. LDU >= 1. |
|
* If singular vectors are desired, then LDU >= max( 1, N ). |
|
* |
|
* VT (output) DOUBLE PRECISION array, dimension (LDVT,N) |
|
* If COMPQ = 'I', then: |
|
* On exit, if INFO = 0, VT**T contains the right singular |
|
* vectors of the bidiagonal matrix. |
|
* For other values of COMPQ, VT is not referenced. |
|
* |
|
* LDVT (input) INTEGER |
|
* The leading dimension of the array VT. LDVT >= 1. |
|
* If singular vectors are desired, then LDVT >= max( 1, N ). |
|
* |
|
* Q (output) DOUBLE PRECISION array, dimension (LDQ) |
|
* If COMPQ = 'P', then: |
|
* On exit, if INFO = 0, Q and IQ contain the left |
|
* and right singular vectors in a compact form, |
|
* requiring O(N log N) space instead of 2*N**2. |
|
* In particular, Q contains all the DOUBLE PRECISION data in |
|
* LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) |
|
* words of memory, where SMLSIZ is returned by ILAENV and |
|
* is equal to the maximum size of the subproblems at the |
|
* bottom of the computation tree (usually about 25). |
|
* For other values of COMPQ, Q is not referenced. |
|
* |
|
* IQ (output) INTEGER array, dimension (LDIQ) |
|
* If COMPQ = 'P', then: |
|
* On exit, if INFO = 0, Q and IQ contain the left |
|
* and right singular vectors in a compact form, |
|
* requiring O(N log N) space instead of 2*N**2. |
|
* In particular, IQ contains all INTEGER data in |
|
* LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) |
|
* words of memory, where SMLSIZ is returned by ILAENV and |
|
* is equal to the maximum size of the subproblems at the |
|
* bottom of the computation tree (usually about 25). |
|
* For other values of COMPQ, IQ is not referenced. |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
* If COMPQ = 'N' then LWORK >= (4 * N). |
|
* If COMPQ = 'P' then LWORK >= (6 * N). |
|
* If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (8*N) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit. |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: The algorithm failed to compute a singular value. |
|
* The update process of divide and conquer failed. |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Based on contributions by |
|
* Ming Gu and Huan Ren, Computer Science Division, University of |
|
* California at Berkeley, USA |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* Changed dimension statement in comment describing E from (N) to |
* Changed dimension statement in comment describing E from (N) to |
* (N-1). Sven, 17 Feb 05. |
* (N-1). Sven, 17 Feb 05. |
Line 225
|
Line 313
|
END IF |
END IF |
IF( IUPLO.EQ.2 ) THEN |
IF( IUPLO.EQ.2 ) THEN |
QSTART = 5 |
QSTART = 5 |
WSTART = 2*N - 1 |
IF( ICOMPQ .EQ. 2 ) WSTART = 2*N - 1 |
DO 10 I = 1, N - 1 |
DO 10 I = 1, N - 1 |
CALL DLARTG( D( I ), E( I ), CS, SN, R ) |
CALL DLARTG( D( I ), E( I ), CS, SN, R ) |
D( I ) = R |
D( I ) = R |
Line 244
|
Line 332
|
* If ICOMPQ = 0, use DLASDQ to compute the singular values. |
* If ICOMPQ = 0, use DLASDQ to compute the singular values. |
* |
* |
IF( ICOMPQ.EQ.0 ) THEN |
IF( ICOMPQ.EQ.0 ) THEN |
|
* Ignore WSTART, instead using WORK( 1 ), since the two vectors |
|
* for CS and -SN above are added only if ICOMPQ == 2, |
|
* and adding them exceeds documented WORK size of 4*n. |
CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U, |
CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U, |
$ LDU, WORK( WSTART ), INFO ) |
$ LDU, WORK( 1 ), INFO ) |
GO TO 40 |
GO TO 40 |
END IF |
END IF |
* |
* |
Line 321
|
Line 412
|
DO 30 I = 1, NM1 |
DO 30 I = 1, NM1 |
IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN |
IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN |
* |
* |
* Subproblem found. First determine its size and then |
* Subproblem found. First determine its size and then |
* apply divide and conquer on it. |
* apply divide and conquer on it. |
* |
* |
IF( I.LT.NM1 ) THEN |
IF( I.LT.NM1 ) THEN |
* |
* |
* A subproblem with E(I) small for I < NM1. |
* A subproblem with E(I) small for I < NM1. |
* |
* |
NSIZE = I - START + 1 |
NSIZE = I - START + 1 |
ELSE IF( ABS( E( I ) ).GE.EPS ) THEN |
ELSE IF( ABS( E( I ) ).GE.EPS ) THEN |
* |
* |
* A subproblem with E(NM1) not too small but I = NM1. |
* A subproblem with E(NM1) not too small but I = NM1. |
* |
* |
NSIZE = N - START + 1 |
NSIZE = N - START + 1 |
ELSE |
ELSE |
* |
* |
* A subproblem with E(NM1) small. This implies an |
* A subproblem with E(NM1) small. This implies an |
* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem |
* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem |
* first. |
* first. |
* |
* |
NSIZE = I - START + 1 |
NSIZE = I - START + 1 |
IF( ICOMPQ.EQ.2 ) THEN |
IF( ICOMPQ.EQ.2 ) THEN |