![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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