LAPACK: SRC/dbdsdc.f Source File

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