![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASDQ 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASDQ + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasdq.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasdq.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasdq.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, 00022 * U, LDU, C, LDC, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ), 00030 * $ VT( LDVT, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SLASDQ computes the singular value decomposition (SVD) of a real 00040 *> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal 00041 *> E, accumulating the transformations if desired. Letting B denote 00042 *> the input bidiagonal matrix, the algorithm computes orthogonal 00043 *> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose 00044 *> of P). The singular values S are overwritten on D. 00045 *> 00046 *> The input matrix U is changed to U * Q if desired. 00047 *> The input matrix VT is changed to P**T * VT if desired. 00048 *> The input matrix C is changed to Q**T * C if desired. 00049 *> 00050 *> See "Computing Small Singular Values of Bidiagonal Matrices With 00051 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 00052 *> LAPACK Working Note #3, for a detailed description of the algorithm. 00053 *> \endverbatim 00054 * 00055 * Arguments: 00056 * ========== 00057 * 00058 *> \param[in] UPLO 00059 *> \verbatim 00060 *> UPLO is CHARACTER*1 00061 *> On entry, UPLO specifies whether the input bidiagonal matrix 00062 *> is upper or lower bidiagonal, and wether it is square are 00063 *> not. 00064 *> UPLO = 'U' or 'u' B is upper bidiagonal. 00065 *> UPLO = 'L' or 'l' B is lower bidiagonal. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] SQRE 00069 *> \verbatim 00070 *> SQRE is INTEGER 00071 *> = 0: then the input matrix is N-by-N. 00072 *> = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and 00073 *> (N+1)-by-N if UPLU = 'L'. 00074 *> 00075 *> The bidiagonal matrix has 00076 *> N = NL + NR + 1 rows and 00077 *> M = N + SQRE >= N columns. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] N 00081 *> \verbatim 00082 *> N is INTEGER 00083 *> On entry, N specifies the number of rows and columns 00084 *> in the matrix. N must be at least 0. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] NCVT 00088 *> \verbatim 00089 *> NCVT is INTEGER 00090 *> On entry, NCVT specifies the number of columns of 00091 *> the matrix VT. NCVT must be at least 0. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] NRU 00095 *> \verbatim 00096 *> NRU is INTEGER 00097 *> On entry, NRU specifies the number of rows of 00098 *> the matrix U. NRU must be at least 0. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] NCC 00102 *> \verbatim 00103 *> NCC is INTEGER 00104 *> On entry, NCC specifies the number of columns of 00105 *> the matrix C. NCC must be at least 0. 00106 *> \endverbatim 00107 *> 00108 *> \param[in,out] D 00109 *> \verbatim 00110 *> D is REAL array, dimension (N) 00111 *> On entry, D contains the diagonal entries of the 00112 *> bidiagonal matrix whose SVD is desired. On normal exit, 00113 *> D contains the singular values in ascending order. 00114 *> \endverbatim 00115 *> 00116 *> \param[in,out] E 00117 *> \verbatim 00118 *> E is REAL array. 00119 *> dimension is (N-1) if SQRE = 0 and N if SQRE = 1. 00120 *> On entry, the entries of E contain the offdiagonal entries 00121 *> of the bidiagonal matrix whose SVD is desired. On normal 00122 *> exit, E will contain 0. If the algorithm does not converge, 00123 *> D and E will contain the diagonal and superdiagonal entries 00124 *> of a bidiagonal matrix orthogonally equivalent to the one 00125 *> given as input. 00126 *> \endverbatim 00127 *> 00128 *> \param[in,out] VT 00129 *> \verbatim 00130 *> VT is REAL array, dimension (LDVT, NCVT) 00131 *> On entry, contains a matrix which on exit has been 00132 *> premultiplied by P**T, dimension N-by-NCVT if SQRE = 0 00133 *> and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). 00134 *> \endverbatim 00135 *> 00136 *> \param[in] LDVT 00137 *> \verbatim 00138 *> LDVT is INTEGER 00139 *> On entry, LDVT specifies the leading dimension of VT as 00140 *> declared in the calling (sub) program. LDVT must be at 00141 *> least 1. If NCVT is nonzero LDVT must also be at least N. 00142 *> \endverbatim 00143 *> 00144 *> \param[in,out] U 00145 *> \verbatim 00146 *> U is REAL array, dimension (LDU, N) 00147 *> On entry, contains a matrix which on exit has been 00148 *> postmultiplied by Q, dimension NRU-by-N if SQRE = 0 00149 *> and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). 00150 *> \endverbatim 00151 *> 00152 *> \param[in] LDU 00153 *> \verbatim 00154 *> LDU is INTEGER 00155 *> On entry, LDU specifies the leading dimension of U as 00156 *> declared in the calling (sub) program. LDU must be at 00157 *> least max( 1, NRU ) . 00158 *> \endverbatim 00159 *> 00160 *> \param[in,out] C 00161 *> \verbatim 00162 *> C is REAL array, dimension (LDC, NCC) 00163 *> On entry, contains an N-by-NCC matrix which on exit 00164 *> has been premultiplied by Q**T dimension N-by-NCC if SQRE = 0 00165 *> and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). 00166 *> \endverbatim 00167 *> 00168 *> \param[in] LDC 00169 *> \verbatim 00170 *> LDC is INTEGER 00171 *> On entry, LDC specifies the leading dimension of C as 00172 *> declared in the calling (sub) program. LDC must be at 00173 *> least 1. If NCC is nonzero, LDC must also be at least N. 00174 *> \endverbatim 00175 *> 00176 *> \param[out] WORK 00177 *> \verbatim 00178 *> WORK is REAL array, dimension (4*N) 00179 *> Workspace. Only referenced if one of NCVT, NRU, or NCC is 00180 *> nonzero, and if N is at least 2. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] INFO 00184 *> \verbatim 00185 *> INFO is INTEGER 00186 *> On exit, a value of 0 indicates a successful exit. 00187 *> If INFO < 0, argument number -INFO is illegal. 00188 *> If INFO > 0, the algorithm did not converge, and INFO 00189 *> specifies how many superdiagonals did not converge. 00190 *> \endverbatim 00191 * 00192 * Authors: 00193 * ======== 00194 * 00195 *> \author Univ. of Tennessee 00196 *> \author Univ. of California Berkeley 00197 *> \author Univ. of Colorado Denver 00198 *> \author NAG Ltd. 00199 * 00200 *> \date November 2011 00201 * 00202 *> \ingroup auxOTHERauxiliary 00203 * 00204 *> \par Contributors: 00205 * ================== 00206 *> 00207 *> Ming Gu and Huan Ren, Computer Science Division, University of 00208 *> California at Berkeley, USA 00209 *> 00210 * ===================================================================== 00211 SUBROUTINE SLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, 00212 $ U, LDU, C, LDC, WORK, INFO ) 00213 * 00214 * -- LAPACK auxiliary routine (version 3.4.0) -- 00215 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00217 * November 2011 00218 * 00219 * .. Scalar Arguments .. 00220 CHARACTER UPLO 00221 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE 00222 * .. 00223 * .. Array Arguments .. 00224 REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ), 00225 $ VT( LDVT, * ), WORK( * ) 00226 * .. 00227 * 00228 * ===================================================================== 00229 * 00230 * .. Parameters .. 00231 REAL ZERO 00232 PARAMETER ( ZERO = 0.0E+0 ) 00233 * .. 00234 * .. Local Scalars .. 00235 LOGICAL ROTATE 00236 INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 00237 REAL CS, R, SMIN, SN 00238 * .. 00239 * .. External Subroutines .. 00240 EXTERNAL SBDSQR, SLARTG, SLASR, SSWAP, XERBLA 00241 * .. 00242 * .. External Functions .. 00243 LOGICAL LSAME 00244 EXTERNAL LSAME 00245 * .. 00246 * .. Intrinsic Functions .. 00247 INTRINSIC MAX 00248 * .. 00249 * .. Executable Statements .. 00250 * 00251 * Test the input parameters. 00252 * 00253 INFO = 0 00254 IUPLO = 0 00255 IF( LSAME( UPLO, 'U' ) ) 00256 $ IUPLO = 1 00257 IF( LSAME( UPLO, 'L' ) ) 00258 $ IUPLO = 2 00259 IF( IUPLO.EQ.0 ) THEN 00260 INFO = -1 00261 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00262 INFO = -2 00263 ELSE IF( N.LT.0 ) THEN 00264 INFO = -3 00265 ELSE IF( NCVT.LT.0 ) THEN 00266 INFO = -4 00267 ELSE IF( NRU.LT.0 ) THEN 00268 INFO = -5 00269 ELSE IF( NCC.LT.0 ) THEN 00270 INFO = -6 00271 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 00272 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 00273 INFO = -10 00274 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 00275 INFO = -12 00276 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 00277 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 00278 INFO = -14 00279 END IF 00280 IF( INFO.NE.0 ) THEN 00281 CALL XERBLA( 'SLASDQ', -INFO ) 00282 RETURN 00283 END IF 00284 IF( N.EQ.0 ) 00285 $ RETURN 00286 * 00287 * ROTATE is true if any singular vectors desired, false otherwise 00288 * 00289 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 00290 NP1 = N + 1 00291 SQRE1 = SQRE 00292 * 00293 * If matrix non-square upper bidiagonal, rotate to be lower 00294 * bidiagonal. The rotations are on the right. 00295 * 00296 IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN 00297 DO 10 I = 1, N - 1 00298 CALL SLARTG( D( I ), E( I ), CS, SN, R ) 00299 D( I ) = R 00300 E( I ) = SN*D( I+1 ) 00301 D( I+1 ) = CS*D( I+1 ) 00302 IF( ROTATE ) THEN 00303 WORK( I ) = CS 00304 WORK( N+I ) = SN 00305 END IF 00306 10 CONTINUE 00307 CALL SLARTG( D( N ), E( N ), CS, SN, R ) 00308 D( N ) = R 00309 E( N ) = ZERO 00310 IF( ROTATE ) THEN 00311 WORK( N ) = CS 00312 WORK( N+N ) = SN 00313 END IF 00314 IUPLO = 2 00315 SQRE1 = 0 00316 * 00317 * Update singular vectors if desired. 00318 * 00319 IF( NCVT.GT.0 ) 00320 $ CALL SLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), 00321 $ WORK( NP1 ), VT, LDVT ) 00322 END IF 00323 * 00324 * If matrix lower bidiagonal, rotate to be upper bidiagonal 00325 * by applying Givens rotations on the left. 00326 * 00327 IF( IUPLO.EQ.2 ) THEN 00328 DO 20 I = 1, N - 1 00329 CALL SLARTG( D( I ), E( I ), CS, SN, R ) 00330 D( I ) = R 00331 E( I ) = SN*D( I+1 ) 00332 D( I+1 ) = CS*D( I+1 ) 00333 IF( ROTATE ) THEN 00334 WORK( I ) = CS 00335 WORK( N+I ) = SN 00336 END IF 00337 20 CONTINUE 00338 * 00339 * If matrix (N+1)-by-N lower bidiagonal, one additional 00340 * rotation is needed. 00341 * 00342 IF( SQRE1.EQ.1 ) THEN 00343 CALL SLARTG( D( N ), E( N ), CS, SN, R ) 00344 D( N ) = R 00345 IF( ROTATE ) THEN 00346 WORK( N ) = CS 00347 WORK( N+N ) = SN 00348 END IF 00349 END IF 00350 * 00351 * Update singular vectors if desired. 00352 * 00353 IF( NRU.GT.0 ) THEN 00354 IF( SQRE1.EQ.0 ) THEN 00355 CALL SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), 00356 $ WORK( NP1 ), U, LDU ) 00357 ELSE 00358 CALL SLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), 00359 $ WORK( NP1 ), U, LDU ) 00360 END IF 00361 END IF 00362 IF( NCC.GT.0 ) THEN 00363 IF( SQRE1.EQ.0 ) THEN 00364 CALL SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), 00365 $ WORK( NP1 ), C, LDC ) 00366 ELSE 00367 CALL SLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), 00368 $ WORK( NP1 ), C, LDC ) 00369 END IF 00370 END IF 00371 END IF 00372 * 00373 * Call SBDSQR to compute the SVD of the reduced real 00374 * N-by-N upper bidiagonal matrix. 00375 * 00376 CALL SBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, 00377 $ LDC, WORK, INFO ) 00378 * 00379 * Sort the singular values into ascending order (insertion sort on 00380 * singular values, but only one transposition per singular vector) 00381 * 00382 DO 40 I = 1, N 00383 * 00384 * Scan for smallest D(I). 00385 * 00386 ISUB = I 00387 SMIN = D( I ) 00388 DO 30 J = I + 1, N 00389 IF( D( J ).LT.SMIN ) THEN 00390 ISUB = J 00391 SMIN = D( J ) 00392 END IF 00393 30 CONTINUE 00394 IF( ISUB.NE.I ) THEN 00395 * 00396 * Swap singular values and vectors. 00397 * 00398 D( ISUB ) = D( I ) 00399 D( I ) = SMIN 00400 IF( NCVT.GT.0 ) 00401 $ CALL SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) 00402 IF( NRU.GT.0 ) 00403 $ CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) 00404 IF( NCC.GT.0 ) 00405 $ CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) 00406 END IF 00407 40 CONTINUE 00408 * 00409 RETURN 00410 * 00411 * End of SLASDQ 00412 * 00413 END