![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASD7 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASD7 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd7.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd7.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd7.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, 00022 * VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, 00023 * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00024 * C, S, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00028 * $ NR, SQRE 00029 * REAL ALPHA, BETA, C, S 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ), 00033 * $ IDXQ( * ), PERM( * ) 00034 * REAL D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ), 00035 * $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ), 00036 * $ ZW( * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> SLASD7 merges the two sets of singular values together into a single 00046 *> sorted set. Then it tries to deflate the size of the problem. There 00047 *> are two ways in which deflation can occur: when two or more singular 00048 *> values are close together or if there is a tiny entry in the Z 00049 *> vector. For each such occurrence the order of the related 00050 *> secular equation problem is reduced by one. 00051 *> 00052 *> SLASD7 is called from SLASD6. 00053 *> \endverbatim 00054 * 00055 * Arguments: 00056 * ========== 00057 * 00058 *> \param[in] ICOMPQ 00059 *> \verbatim 00060 *> ICOMPQ is INTEGER 00061 *> Specifies whether singular vectors are to be computed 00062 *> in compact form, as follows: 00063 *> = 0: Compute singular values only. 00064 *> = 1: Compute singular vectors of upper 00065 *> bidiagonal matrix in compact form. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] NL 00069 *> \verbatim 00070 *> NL is INTEGER 00071 *> The row dimension of the upper block. NL >= 1. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] NR 00075 *> \verbatim 00076 *> NR is INTEGER 00077 *> The row dimension of the lower block. NR >= 1. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] SQRE 00081 *> \verbatim 00082 *> SQRE is INTEGER 00083 *> = 0: the lower block is an NR-by-NR square matrix. 00084 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00085 *> 00086 *> The bidiagonal matrix has 00087 *> N = NL + NR + 1 rows and 00088 *> M = N + SQRE >= N columns. 00089 *> \endverbatim 00090 *> 00091 *> \param[out] K 00092 *> \verbatim 00093 *> K is INTEGER 00094 *> Contains the dimension of the non-deflated matrix, this is 00095 *> the order of the related secular equation. 1 <= K <=N. 00096 *> \endverbatim 00097 *> 00098 *> \param[in,out] D 00099 *> \verbatim 00100 *> D is REAL array, dimension ( N ) 00101 *> On entry D contains the singular values of the two submatrices 00102 *> to be combined. On exit D contains the trailing (N-K) updated 00103 *> singular values (those which were deflated) sorted into 00104 *> increasing order. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Z 00108 *> \verbatim 00109 *> Z is REAL array, dimension ( M ) 00110 *> On exit Z contains the updating row vector in the secular 00111 *> equation. 00112 *> \endverbatim 00113 *> 00114 *> \param[out] ZW 00115 *> \verbatim 00116 *> ZW is REAL array, dimension ( M ) 00117 *> Workspace for Z. 00118 *> \endverbatim 00119 *> 00120 *> \param[in,out] VF 00121 *> \verbatim 00122 *> VF is REAL array, dimension ( M ) 00123 *> On entry, VF(1:NL+1) contains the first components of all 00124 *> right singular vectors of the upper block; and VF(NL+2:M) 00125 *> contains the first components of all right singular vectors 00126 *> of the lower block. On exit, VF contains the first components 00127 *> of all right singular vectors of the bidiagonal matrix. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] VFW 00131 *> \verbatim 00132 *> VFW is REAL array, dimension ( M ) 00133 *> Workspace for VF. 00134 *> \endverbatim 00135 *> 00136 *> \param[in,out] VL 00137 *> \verbatim 00138 *> VL is REAL array, dimension ( M ) 00139 *> On entry, VL(1:NL+1) contains the last components of all 00140 *> right singular vectors of the upper block; and VL(NL+2:M) 00141 *> contains the last components of all right singular vectors 00142 *> of the lower block. On exit, VL contains the last components 00143 *> of all right singular vectors of the bidiagonal matrix. 00144 *> \endverbatim 00145 *> 00146 *> \param[out] VLW 00147 *> \verbatim 00148 *> VLW is REAL array, dimension ( M ) 00149 *> Workspace for VL. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] ALPHA 00153 *> \verbatim 00154 *> ALPHA is REAL 00155 *> Contains the diagonal element associated with the added row. 00156 *> \endverbatim 00157 *> 00158 *> \param[in] BETA 00159 *> \verbatim 00160 *> BETA is REAL 00161 *> Contains the off-diagonal element associated with the added 00162 *> row. 00163 *> \endverbatim 00164 *> 00165 *> \param[out] DSIGMA 00166 *> \verbatim 00167 *> DSIGMA is REAL array, dimension ( N ) 00168 *> Contains a copy of the diagonal elements (K-1 singular values 00169 *> and one zero) in the secular equation. 00170 *> \endverbatim 00171 *> 00172 *> \param[out] IDX 00173 *> \verbatim 00174 *> IDX is INTEGER array, dimension ( N ) 00175 *> This will contain the permutation used to sort the contents of 00176 *> D into ascending order. 00177 *> \endverbatim 00178 *> 00179 *> \param[out] IDXP 00180 *> \verbatim 00181 *> IDXP is INTEGER array, dimension ( N ) 00182 *> This will contain the permutation used to place deflated 00183 *> values of D at the end of the array. On output IDXP(2:K) 00184 *> points to the nondeflated D-values and IDXP(K+1:N) 00185 *> points to the deflated singular values. 00186 *> \endverbatim 00187 *> 00188 *> \param[in] IDXQ 00189 *> \verbatim 00190 *> IDXQ is INTEGER array, dimension ( N ) 00191 *> This contains the permutation which separately sorts the two 00192 *> sub-problems in D into ascending order. Note that entries in 00193 *> the first half of this permutation must first be moved one 00194 *> position backward; and entries in the second half 00195 *> must first have NL+1 added to their values. 00196 *> \endverbatim 00197 *> 00198 *> \param[out] PERM 00199 *> \verbatim 00200 *> PERM is INTEGER array, dimension ( N ) 00201 *> The permutations (from deflation and sorting) to be applied 00202 *> to each singular block. Not referenced if ICOMPQ = 0. 00203 *> \endverbatim 00204 *> 00205 *> \param[out] GIVPTR 00206 *> \verbatim 00207 *> GIVPTR is INTEGER 00208 *> The number of Givens rotations which took place in this 00209 *> subproblem. Not referenced if ICOMPQ = 0. 00210 *> \endverbatim 00211 *> 00212 *> \param[out] GIVCOL 00213 *> \verbatim 00214 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 ) 00215 *> Each pair of numbers indicates a pair of columns to take place 00216 *> in a Givens rotation. Not referenced if ICOMPQ = 0. 00217 *> \endverbatim 00218 *> 00219 *> \param[in] LDGCOL 00220 *> \verbatim 00221 *> LDGCOL is INTEGER 00222 *> The leading dimension of GIVCOL, must be at least N. 00223 *> \endverbatim 00224 *> 00225 *> \param[out] GIVNUM 00226 *> \verbatim 00227 *> GIVNUM is REAL array, dimension ( LDGNUM, 2 ) 00228 *> Each number indicates the C or S value to be used in the 00229 *> corresponding Givens rotation. Not referenced if ICOMPQ = 0. 00230 *> \endverbatim 00231 *> 00232 *> \param[in] LDGNUM 00233 *> \verbatim 00234 *> LDGNUM is INTEGER 00235 *> The leading dimension of GIVNUM, must be at least N. 00236 *> \endverbatim 00237 *> 00238 *> \param[out] C 00239 *> \verbatim 00240 *> C is REAL 00241 *> C contains garbage if SQRE =0 and the C-value of a Givens 00242 *> rotation related to the right null space if SQRE = 1. 00243 *> \endverbatim 00244 *> 00245 *> \param[out] S 00246 *> \verbatim 00247 *> S is REAL 00248 *> S contains garbage if SQRE =0 and the S-value of a Givens 00249 *> rotation related to the right null space if SQRE = 1. 00250 *> \endverbatim 00251 *> 00252 *> \param[out] INFO 00253 *> \verbatim 00254 *> INFO is INTEGER 00255 *> = 0: successful exit. 00256 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00257 *> \endverbatim 00258 * 00259 * Authors: 00260 * ======== 00261 * 00262 *> \author Univ. of Tennessee 00263 *> \author Univ. of California Berkeley 00264 *> \author Univ. of Colorado Denver 00265 *> \author NAG Ltd. 00266 * 00267 *> \date November 2011 00268 * 00269 *> \ingroup auxOTHERauxiliary 00270 * 00271 *> \par Contributors: 00272 * ================== 00273 *> 00274 *> Ming Gu and Huan Ren, Computer Science Division, University of 00275 *> California at Berkeley, USA 00276 *> 00277 * ===================================================================== 00278 SUBROUTINE SLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, 00279 $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, 00280 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00281 $ C, S, INFO ) 00282 * 00283 * -- LAPACK auxiliary routine (version 3.4.0) -- 00284 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00285 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00286 * November 2011 00287 * 00288 * .. Scalar Arguments .. 00289 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00290 $ NR, SQRE 00291 REAL ALPHA, BETA, C, S 00292 * .. 00293 * .. Array Arguments .. 00294 INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ), 00295 $ IDXQ( * ), PERM( * ) 00296 REAL D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ), 00297 $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ), 00298 $ ZW( * ) 00299 * .. 00300 * 00301 * ===================================================================== 00302 * 00303 * .. Parameters .. 00304 REAL ZERO, ONE, TWO, EIGHT 00305 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0, 00306 $ EIGHT = 8.0E+0 ) 00307 * .. 00308 * .. Local Scalars .. 00309 * 00310 INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N, 00311 $ NLP1, NLP2 00312 REAL EPS, HLFTOL, TAU, TOL, Z1 00313 * .. 00314 * .. External Subroutines .. 00315 EXTERNAL SCOPY, SLAMRG, SROT, XERBLA 00316 * .. 00317 * .. External Functions .. 00318 REAL SLAMCH, SLAPY2 00319 EXTERNAL SLAMCH, SLAPY2 00320 * .. 00321 * .. Intrinsic Functions .. 00322 INTRINSIC ABS, MAX 00323 * .. 00324 * .. Executable Statements .. 00325 * 00326 * Test the input parameters. 00327 * 00328 INFO = 0 00329 N = NL + NR + 1 00330 M = N + SQRE 00331 * 00332 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00333 INFO = -1 00334 ELSE IF( NL.LT.1 ) THEN 00335 INFO = -2 00336 ELSE IF( NR.LT.1 ) THEN 00337 INFO = -3 00338 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00339 INFO = -4 00340 ELSE IF( LDGCOL.LT.N ) THEN 00341 INFO = -22 00342 ELSE IF( LDGNUM.LT.N ) THEN 00343 INFO = -24 00344 END IF 00345 IF( INFO.NE.0 ) THEN 00346 CALL XERBLA( 'SLASD7', -INFO ) 00347 RETURN 00348 END IF 00349 * 00350 NLP1 = NL + 1 00351 NLP2 = NL + 2 00352 IF( ICOMPQ.EQ.1 ) THEN 00353 GIVPTR = 0 00354 END IF 00355 * 00356 * Generate the first part of the vector Z and move the singular 00357 * values in the first part of D one position backward. 00358 * 00359 Z1 = ALPHA*VL( NLP1 ) 00360 VL( NLP1 ) = ZERO 00361 TAU = VF( NLP1 ) 00362 DO 10 I = NL, 1, -1 00363 Z( I+1 ) = ALPHA*VL( I ) 00364 VL( I ) = ZERO 00365 VF( I+1 ) = VF( I ) 00366 D( I+1 ) = D( I ) 00367 IDXQ( I+1 ) = IDXQ( I ) + 1 00368 10 CONTINUE 00369 VF( 1 ) = TAU 00370 * 00371 * Generate the second part of the vector Z. 00372 * 00373 DO 20 I = NLP2, M 00374 Z( I ) = BETA*VF( I ) 00375 VF( I ) = ZERO 00376 20 CONTINUE 00377 * 00378 * Sort the singular values into increasing order 00379 * 00380 DO 30 I = NLP2, N 00381 IDXQ( I ) = IDXQ( I ) + NLP1 00382 30 CONTINUE 00383 * 00384 * DSIGMA, IDXC, IDXC, and ZW are used as storage space. 00385 * 00386 DO 40 I = 2, N 00387 DSIGMA( I ) = D( IDXQ( I ) ) 00388 ZW( I ) = Z( IDXQ( I ) ) 00389 VFW( I ) = VF( IDXQ( I ) ) 00390 VLW( I ) = VL( IDXQ( I ) ) 00391 40 CONTINUE 00392 * 00393 CALL SLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) ) 00394 * 00395 DO 50 I = 2, N 00396 IDXI = 1 + IDX( I ) 00397 D( I ) = DSIGMA( IDXI ) 00398 Z( I ) = ZW( IDXI ) 00399 VF( I ) = VFW( IDXI ) 00400 VL( I ) = VLW( IDXI ) 00401 50 CONTINUE 00402 * 00403 * Calculate the allowable deflation tolerence 00404 * 00405 EPS = SLAMCH( 'Epsilon' ) 00406 TOL = MAX( ABS( ALPHA ), ABS( BETA ) ) 00407 TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL ) 00408 * 00409 * There are 2 kinds of deflation -- first a value in the z-vector 00410 * is small, second two (or more) singular values are very close 00411 * together (their difference is small). 00412 * 00413 * If the value in the z-vector is small, we simply permute the 00414 * array so that the corresponding singular value is moved to the 00415 * end. 00416 * 00417 * If two values in the D-vector are close, we perform a two-sided 00418 * rotation designed to make one of the corresponding z-vector 00419 * entries zero, and then permute the array so that the deflated 00420 * singular value is moved to the end. 00421 * 00422 * If there are multiple singular values then the problem deflates. 00423 * Here the number of equal singular values are found. As each equal 00424 * singular value is found, an elementary reflector is computed to 00425 * rotate the corresponding singular subspace so that the 00426 * corresponding components of Z are zero in this new basis. 00427 * 00428 K = 1 00429 K2 = N + 1 00430 DO 60 J = 2, N 00431 IF( ABS( Z( J ) ).LE.TOL ) THEN 00432 * 00433 * Deflate due to small z component. 00434 * 00435 K2 = K2 - 1 00436 IDXP( K2 ) = J 00437 IF( J.EQ.N ) 00438 $ GO TO 100 00439 ELSE 00440 JPREV = J 00441 GO TO 70 00442 END IF 00443 60 CONTINUE 00444 70 CONTINUE 00445 J = JPREV 00446 80 CONTINUE 00447 J = J + 1 00448 IF( J.GT.N ) 00449 $ GO TO 90 00450 IF( ABS( Z( J ) ).LE.TOL ) THEN 00451 * 00452 * Deflate due to small z component. 00453 * 00454 K2 = K2 - 1 00455 IDXP( K2 ) = J 00456 ELSE 00457 * 00458 * Check if singular values are close enough to allow deflation. 00459 * 00460 IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN 00461 * 00462 * Deflation is possible. 00463 * 00464 S = Z( JPREV ) 00465 C = Z( J ) 00466 * 00467 * Find sqrt(a**2+b**2) without overflow or 00468 * destructive underflow. 00469 * 00470 TAU = SLAPY2( C, S ) 00471 Z( J ) = TAU 00472 Z( JPREV ) = ZERO 00473 C = C / TAU 00474 S = -S / TAU 00475 * 00476 * Record the appropriate Givens rotation 00477 * 00478 IF( ICOMPQ.EQ.1 ) THEN 00479 GIVPTR = GIVPTR + 1 00480 IDXJP = IDXQ( IDX( JPREV )+1 ) 00481 IDXJ = IDXQ( IDX( J )+1 ) 00482 IF( IDXJP.LE.NLP1 ) THEN 00483 IDXJP = IDXJP - 1 00484 END IF 00485 IF( IDXJ.LE.NLP1 ) THEN 00486 IDXJ = IDXJ - 1 00487 END IF 00488 GIVCOL( GIVPTR, 2 ) = IDXJP 00489 GIVCOL( GIVPTR, 1 ) = IDXJ 00490 GIVNUM( GIVPTR, 2 ) = C 00491 GIVNUM( GIVPTR, 1 ) = S 00492 END IF 00493 CALL SROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S ) 00494 CALL SROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S ) 00495 K2 = K2 - 1 00496 IDXP( K2 ) = JPREV 00497 JPREV = J 00498 ELSE 00499 K = K + 1 00500 ZW( K ) = Z( JPREV ) 00501 DSIGMA( K ) = D( JPREV ) 00502 IDXP( K ) = JPREV 00503 JPREV = J 00504 END IF 00505 END IF 00506 GO TO 80 00507 90 CONTINUE 00508 * 00509 * Record the last singular value. 00510 * 00511 K = K + 1 00512 ZW( K ) = Z( JPREV ) 00513 DSIGMA( K ) = D( JPREV ) 00514 IDXP( K ) = JPREV 00515 * 00516 100 CONTINUE 00517 * 00518 * Sort the singular values into DSIGMA. The singular values which 00519 * were not deflated go into the first K slots of DSIGMA, except 00520 * that DSIGMA(1) is treated separately. 00521 * 00522 DO 110 J = 2, N 00523 JP = IDXP( J ) 00524 DSIGMA( J ) = D( JP ) 00525 VFW( J ) = VF( JP ) 00526 VLW( J ) = VL( JP ) 00527 110 CONTINUE 00528 IF( ICOMPQ.EQ.1 ) THEN 00529 DO 120 J = 2, N 00530 JP = IDXP( J ) 00531 PERM( J ) = IDXQ( IDX( JP )+1 ) 00532 IF( PERM( J ).LE.NLP1 ) THEN 00533 PERM( J ) = PERM( J ) - 1 00534 END IF 00535 120 CONTINUE 00536 END IF 00537 * 00538 * The deflated singular values go back into the last N - K slots of 00539 * D. 00540 * 00541 CALL SCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 ) 00542 * 00543 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and 00544 * VL(M). 00545 * 00546 DSIGMA( 1 ) = ZERO 00547 HLFTOL = TOL / TWO 00548 IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL ) 00549 $ DSIGMA( 2 ) = HLFTOL 00550 IF( M.GT.N ) THEN 00551 Z( 1 ) = SLAPY2( Z1, Z( M ) ) 00552 IF( Z( 1 ).LE.TOL ) THEN 00553 C = ONE 00554 S = ZERO 00555 Z( 1 ) = TOL 00556 ELSE 00557 C = Z1 / Z( 1 ) 00558 S = -Z( M ) / Z( 1 ) 00559 END IF 00560 CALL SROT( 1, VF( M ), 1, VF( 1 ), 1, C, S ) 00561 CALL SROT( 1, VL( M ), 1, VL( 1 ), 1, C, S ) 00562 ELSE 00563 IF( ABS( Z1 ).LE.TOL ) THEN 00564 Z( 1 ) = TOL 00565 ELSE 00566 Z( 1 ) = Z1 00567 END IF 00568 END IF 00569 * 00570 * Restore Z, VF, and VL. 00571 * 00572 CALL SCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 ) 00573 CALL SCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 ) 00574 CALL SCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 ) 00575 * 00576 RETURN 00577 * 00578 * End of SLASD7 00579 * 00580 END