LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dbdsdc.f
Go to the documentation of this file.
00001 *> \brief \b DBDSDC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DBDSDC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
00022 *                          WORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          COMPQ, UPLO
00026 *       INTEGER            INFO, LDU, LDVT, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IQ( * ), IWORK( * )
00030 *       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
00031 *      $                   VT( LDVT, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DBDSDC computes the singular value decomposition (SVD) of a real
00041 *> N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
00042 *> using a divide and conquer method, where S is a diagonal matrix
00043 *> with non-negative diagonal elements (the singular values of B), and
00044 *> U and VT are orthogonal matrices of left and right singular vectors,
00045 *> respectively. DBDSDC can be used to compute all singular values,
00046 *> and optionally, singular vectors or singular vectors in compact form.
00047 *>
00048 *> This code makes very mild assumptions about floating point
00049 *> arithmetic. It will work on machines with a guard digit in
00050 *> add/subtract, or on those binary machines without guard digits
00051 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
00052 *> It could conceivably fail on hexadecimal or decimal machines
00053 *> without guard digits, but we know of none.  See DLASD3 for details.
00054 *>
00055 *> The code currently calls DLASDQ if singular values only are desired.
00056 *> However, it can be slightly modified to compute singular values
00057 *> using the divide and conquer method.
00058 *> \endverbatim
00059 *
00060 *  Arguments:
00061 *  ==========
00062 *
00063 *> \param[in] UPLO
00064 *> \verbatim
00065 *>          UPLO is CHARACTER*1
00066 *>          = 'U':  B is upper bidiagonal.
00067 *>          = 'L':  B is lower bidiagonal.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] COMPQ
00071 *> \verbatim
00072 *>          COMPQ is CHARACTER*1
00073 *>          Specifies whether singular vectors are to be computed
00074 *>          as follows:
00075 *>          = 'N':  Compute singular values only;
00076 *>          = 'P':  Compute singular values and compute singular
00077 *>                  vectors in compact form;
00078 *>          = 'I':  Compute singular values and singular vectors.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] N
00082 *> \verbatim
00083 *>          N is INTEGER
00084 *>          The order of the matrix B.  N >= 0.
00085 *> \endverbatim
00086 *>
00087 *> \param[in,out] D
00088 *> \verbatim
00089 *>          D is DOUBLE PRECISION array, dimension (N)
00090 *>          On entry, the n diagonal elements of the bidiagonal matrix B.
00091 *>          On exit, if INFO=0, the singular values of B.
00092 *> \endverbatim
00093 *>
00094 *> \param[in,out] E
00095 *> \verbatim
00096 *>          E is DOUBLE PRECISION array, dimension (N-1)
00097 *>          On entry, the elements of E contain the offdiagonal
00098 *>          elements of the bidiagonal matrix whose SVD is desired.
00099 *>          On exit, E has been destroyed.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] U
00103 *> \verbatim
00104 *>          U is DOUBLE PRECISION array, dimension (LDU,N)
00105 *>          If  COMPQ = 'I', then:
00106 *>             On exit, if INFO = 0, U contains the left singular vectors
00107 *>             of the bidiagonal matrix.
00108 *>          For other values of COMPQ, U is not referenced.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDU
00112 *> \verbatim
00113 *>          LDU is INTEGER
00114 *>          The leading dimension of the array U.  LDU >= 1.
00115 *>          If singular vectors are desired, then LDU >= max( 1, N ).
00116 *> \endverbatim
00117 *>
00118 *> \param[out] VT
00119 *> \verbatim
00120 *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
00121 *>          If  COMPQ = 'I', then:
00122 *>             On exit, if INFO = 0, VT**T contains the right singular
00123 *>             vectors of the bidiagonal matrix.
00124 *>          For other values of COMPQ, VT is not referenced.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] LDVT
00128 *> \verbatim
00129 *>          LDVT is INTEGER
00130 *>          The leading dimension of the array VT.  LDVT >= 1.
00131 *>          If singular vectors are desired, then LDVT >= max( 1, N ).
00132 *> \endverbatim
00133 *>
00134 *> \param[out] Q
00135 *> \verbatim
00136 *>          Q is DOUBLE PRECISION array, dimension (LDQ)
00137 *>          If  COMPQ = 'P', then:
00138 *>             On exit, if INFO = 0, Q and IQ contain the left
00139 *>             and right singular vectors in a compact form,
00140 *>             requiring O(N log N) space instead of 2*N**2.
00141 *>             In particular, Q contains all the DOUBLE PRECISION data in
00142 *>             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
00143 *>             words of memory, where SMLSIZ is returned by ILAENV and
00144 *>             is equal to the maximum size of the subproblems at the
00145 *>             bottom of the computation tree (usually about 25).
00146 *>          For other values of COMPQ, Q is not referenced.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] IQ
00150 *> \verbatim
00151 *>          IQ is INTEGER array, dimension (LDIQ)
00152 *>          If  COMPQ = 'P', then:
00153 *>             On exit, if INFO = 0, Q and IQ contain the left
00154 *>             and right singular vectors in a compact form,
00155 *>             requiring O(N log N) space instead of 2*N**2.
00156 *>             In particular, IQ contains all INTEGER data in
00157 *>             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
00158 *>             words of memory, where SMLSIZ is returned by ILAENV and
00159 *>             is equal to the maximum size of the subproblems at the
00160 *>             bottom of the computation tree (usually about 25).
00161 *>          For other values of COMPQ, IQ is not referenced.
00162 *> \endverbatim
00163 *>
00164 *> \param[out] WORK
00165 *> \verbatim
00166 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00167 *>          If COMPQ = 'N' then LWORK >= (4 * N).
00168 *>          If COMPQ = 'P' then LWORK >= (6 * N).
00169 *>          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
00170 *> \endverbatim
00171 *>
00172 *> \param[out] IWORK
00173 *> \verbatim
00174 *>          IWORK is INTEGER array, dimension (8*N)
00175 *> \endverbatim
00176 *>
00177 *> \param[out] INFO
00178 *> \verbatim
00179 *>          INFO is INTEGER
00180 *>          = 0:  successful exit.
00181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00182 *>          > 0:  The algorithm failed to compute a singular value.
00183 *>                The update process of divide and conquer failed.
00184 *> \endverbatim
00185 *
00186 *  Authors:
00187 *  ========
00188 *
00189 *> \author Univ. of Tennessee 
00190 *> \author Univ. of California Berkeley 
00191 *> \author Univ. of Colorado Denver 
00192 *> \author NAG Ltd. 
00193 *
00194 *> \date November 2011
00195 *
00196 *> \ingroup auxOTHERcomputational
00197 *
00198 *> \par Contributors:
00199 *  ==================
00200 *>
00201 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00202 *>     California at Berkeley, USA
00203 *>
00204 *  =====================================================================
00205       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
00206      $                   WORK, IWORK, INFO )
00207 *
00208 *  -- LAPACK computational routine (version 3.4.0) --
00209 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00210 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00211 *     November 2011
00212 *
00213 *     .. Scalar Arguments ..
00214       CHARACTER          COMPQ, UPLO
00215       INTEGER            INFO, LDU, LDVT, N
00216 *     ..
00217 *     .. Array Arguments ..
00218       INTEGER            IQ( * ), IWORK( * )
00219       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
00220      $                   VT( LDVT, * ), WORK( * )
00221 *     ..
00222 *
00223 *  =====================================================================
00224 *  Changed dimension statement in comment describing E from (N) to
00225 *  (N-1).  Sven, 17 Feb 05.
00226 *  =====================================================================
00227 *
00228 *     .. Parameters ..
00229       DOUBLE PRECISION   ZERO, ONE, TWO
00230       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
00231 *     ..
00232 *     .. Local Scalars ..
00233       INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
00234      $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
00235      $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
00236      $                   SMLSZP, SQRE, START, WSTART, Z
00237       DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
00238 *     ..
00239 *     .. External Functions ..
00240       LOGICAL            LSAME
00241       INTEGER            ILAENV
00242       DOUBLE PRECISION   DLAMCH, DLANST
00243       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
00244 *     ..
00245 *     .. External Subroutines ..
00246       EXTERNAL           DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
00247      $                   DLASET, DLASR, DSWAP, XERBLA
00248 *     ..
00249 *     .. Intrinsic Functions ..
00250       INTRINSIC          ABS, DBLE, INT, LOG, SIGN
00251 *     ..
00252 *     .. Executable Statements ..
00253 *
00254 *     Test the input parameters.
00255 *
00256       INFO = 0
00257 *
00258       IUPLO = 0
00259       IF( LSAME( UPLO, 'U' ) )
00260      $   IUPLO = 1
00261       IF( LSAME( UPLO, 'L' ) )
00262      $   IUPLO = 2
00263       IF( LSAME( COMPQ, 'N' ) ) THEN
00264          ICOMPQ = 0
00265       ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
00266          ICOMPQ = 1
00267       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00268          ICOMPQ = 2
00269       ELSE
00270          ICOMPQ = -1
00271       END IF
00272       IF( IUPLO.EQ.0 ) THEN
00273          INFO = -1
00274       ELSE IF( ICOMPQ.LT.0 ) THEN
00275          INFO = -2
00276       ELSE IF( N.LT.0 ) THEN
00277          INFO = -3
00278       ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
00279      $         N ) ) ) THEN
00280          INFO = -7
00281       ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
00282      $         N ) ) ) THEN
00283          INFO = -9
00284       END IF
00285       IF( INFO.NE.0 ) THEN
00286          CALL XERBLA( 'DBDSDC', -INFO )
00287          RETURN
00288       END IF
00289 *
00290 *     Quick return if possible
00291 *
00292       IF( N.EQ.0 )
00293      $   RETURN
00294       SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
00295       IF( N.EQ.1 ) THEN
00296          IF( ICOMPQ.EQ.1 ) THEN
00297             Q( 1 ) = SIGN( ONE, D( 1 ) )
00298             Q( 1+SMLSIZ*N ) = ONE
00299          ELSE IF( ICOMPQ.EQ.2 ) THEN
00300             U( 1, 1 ) = SIGN( ONE, D( 1 ) )
00301             VT( 1, 1 ) = ONE
00302          END IF
00303          D( 1 ) = ABS( D( 1 ) )
00304          RETURN
00305       END IF
00306       NM1 = N - 1
00307 *
00308 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
00309 *     by applying Givens rotations on the left
00310 *
00311       WSTART = 1
00312       QSTART = 3
00313       IF( ICOMPQ.EQ.1 ) THEN
00314          CALL DCOPY( N, D, 1, Q( 1 ), 1 )
00315          CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
00316       END IF
00317       IF( IUPLO.EQ.2 ) THEN
00318          QSTART = 5
00319          WSTART = 2*N - 1
00320          DO 10 I = 1, N - 1
00321             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00322             D( I ) = R
00323             E( I ) = SN*D( I+1 )
00324             D( I+1 ) = CS*D( I+1 )
00325             IF( ICOMPQ.EQ.1 ) THEN
00326                Q( I+2*N ) = CS
00327                Q( I+3*N ) = SN
00328             ELSE IF( ICOMPQ.EQ.2 ) THEN
00329                WORK( I ) = CS
00330                WORK( NM1+I ) = -SN
00331             END IF
00332    10    CONTINUE
00333       END IF
00334 *
00335 *     If ICOMPQ = 0, use DLASDQ to compute the singular values.
00336 *
00337       IF( ICOMPQ.EQ.0 ) THEN
00338          CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
00339      $                LDU, WORK( WSTART ), INFO )
00340          GO TO 40
00341       END IF
00342 *
00343 *     If N is smaller than the minimum divide size SMLSIZ, then solve
00344 *     the problem with another solver.
00345 *
00346       IF( N.LE.SMLSIZ ) THEN
00347          IF( ICOMPQ.EQ.2 ) THEN
00348             CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
00349             CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
00350             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
00351      $                   LDU, WORK( WSTART ), INFO )
00352          ELSE IF( ICOMPQ.EQ.1 ) THEN
00353             IU = 1
00354             IVT = IU + N
00355             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
00356      $                   N )
00357             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
00358      $                   N )
00359             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
00360      $                   Q( IVT+( QSTART-1 )*N ), N,
00361      $                   Q( IU+( QSTART-1 )*N ), N,
00362      $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
00363      $                   INFO )
00364          END IF
00365          GO TO 40
00366       END IF
00367 *
00368       IF( ICOMPQ.EQ.2 ) THEN
00369          CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
00370          CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
00371       END IF
00372 *
00373 *     Scale.
00374 *
00375       ORGNRM = DLANST( 'M', N, D, E )
00376       IF( ORGNRM.EQ.ZERO )
00377      $   RETURN
00378       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
00379       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
00380 *
00381       EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
00382 *
00383       MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
00384       SMLSZP = SMLSIZ + 1
00385 *
00386       IF( ICOMPQ.EQ.1 ) THEN
00387          IU = 1
00388          IVT = 1 + SMLSIZ
00389          DIFL = IVT + SMLSZP
00390          DIFR = DIFL + MLVL
00391          Z = DIFR + MLVL*2
00392          IC = Z + MLVL
00393          IS = IC + 1
00394          POLES = IS + 1
00395          GIVNUM = POLES + 2*MLVL
00396 *
00397          K = 1
00398          GIVPTR = 2
00399          PERM = 3
00400          GIVCOL = PERM + MLVL
00401       END IF
00402 *
00403       DO 20 I = 1, N
00404          IF( ABS( D( I ) ).LT.EPS ) THEN
00405             D( I ) = SIGN( EPS, D( I ) )
00406          END IF
00407    20 CONTINUE
00408 *
00409       START = 1
00410       SQRE = 0
00411 *
00412       DO 30 I = 1, NM1
00413          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
00414 *
00415 *        Subproblem found. First determine its size and then
00416 *        apply divide and conquer on it.
00417 *
00418             IF( I.LT.NM1 ) THEN
00419 *
00420 *        A subproblem with E(I) small for I < NM1.
00421 *
00422                NSIZE = I - START + 1
00423             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
00424 *
00425 *        A subproblem with E(NM1) not too small but I = NM1.
00426 *
00427                NSIZE = N - START + 1
00428             ELSE
00429 *
00430 *        A subproblem with E(NM1) small. This implies an
00431 *        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
00432 *        first.
00433 *
00434                NSIZE = I - START + 1
00435                IF( ICOMPQ.EQ.2 ) THEN
00436                   U( N, N ) = SIGN( ONE, D( N ) )
00437                   VT( N, N ) = ONE
00438                ELSE IF( ICOMPQ.EQ.1 ) THEN
00439                   Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
00440                   Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
00441                END IF
00442                D( N ) = ABS( D( N ) )
00443             END IF
00444             IF( ICOMPQ.EQ.2 ) THEN
00445                CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
00446      $                      U( START, START ), LDU, VT( START, START ),
00447      $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
00448             ELSE
00449                CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
00450      $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
00451      $                      Q( START+( IVT+QSTART-2 )*N ),
00452      $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
00453      $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
00454      $                      Q( START+( Z+QSTART-2 )*N ),
00455      $                      Q( START+( POLES+QSTART-2 )*N ),
00456      $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
00457      $                      N, IQ( START+PERM*N ),
00458      $                      Q( START+( GIVNUM+QSTART-2 )*N ),
00459      $                      Q( START+( IC+QSTART-2 )*N ),
00460      $                      Q( START+( IS+QSTART-2 )*N ),
00461      $                      WORK( WSTART ), IWORK, INFO )
00462             END IF
00463             IF( INFO.NE.0 ) THEN
00464                RETURN
00465             END IF
00466             START = I + 1
00467          END IF
00468    30 CONTINUE
00469 *
00470 *     Unscale
00471 *
00472       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
00473    40 CONTINUE
00474 *
00475 *     Use Selection Sort to minimize swaps of singular vectors
00476 *
00477       DO 60 II = 2, N
00478          I = II - 1
00479          KK = I
00480          P = D( I )
00481          DO 50 J = II, N
00482             IF( D( J ).GT.P ) THEN
00483                KK = J
00484                P = D( J )
00485             END IF
00486    50    CONTINUE
00487          IF( KK.NE.I ) THEN
00488             D( KK ) = D( I )
00489             D( I ) = P
00490             IF( ICOMPQ.EQ.1 ) THEN
00491                IQ( I ) = KK
00492             ELSE IF( ICOMPQ.EQ.2 ) THEN
00493                CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
00494                CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
00495             END IF
00496          ELSE IF( ICOMPQ.EQ.1 ) THEN
00497             IQ( I ) = I
00498          END IF
00499    60 CONTINUE
00500 *
00501 *     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
00502 *
00503       IF( ICOMPQ.EQ.1 ) THEN
00504          IF( IUPLO.EQ.1 ) THEN
00505             IQ( N ) = 1
00506          ELSE
00507             IQ( N ) = 0
00508          END IF
00509       END IF
00510 *
00511 *     If B is lower bidiagonal, update U by those Givens rotations
00512 *     which rotated B to be upper bidiagonal
00513 *
00514       IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
00515      $   CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
00516 *
00517       RETURN
00518 *
00519 *     End of DBDSDC
00520 *
00521       END
 All Files Functions