|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, 00002 $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, 00003 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00004 $ C, S, INFO ) 00005 * 00006 * -- LAPACK auxiliary routine (version 3.2) -- 00007 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00008 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00009 * November 2006 00010 * 00011 * .. Scalar Arguments .. 00012 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00013 $ NR, SQRE 00014 DOUBLE PRECISION ALPHA, BETA, C, S 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ), 00018 $ IDXQ( * ), PERM( * ) 00019 DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ), 00020 $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ), 00021 $ ZW( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * DLASD7 merges the two sets of singular values together into a single 00028 * sorted set. Then it tries to deflate the size of the problem. There 00029 * are two ways in which deflation can occur: when two or more singular 00030 * values are close together or if there is a tiny entry in the Z 00031 * vector. For each such occurrence the order of the related 00032 * secular equation problem is reduced by one. 00033 * 00034 * DLASD7 is called from DLASD6. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * ICOMPQ (input) INTEGER 00040 * Specifies whether singular vectors are to be computed 00041 * in compact form, as follows: 00042 * = 0: Compute singular values only. 00043 * = 1: Compute singular vectors of upper 00044 * bidiagonal matrix in compact form. 00045 * 00046 * NL (input) INTEGER 00047 * The row dimension of the upper block. NL >= 1. 00048 * 00049 * NR (input) INTEGER 00050 * The row dimension of the lower block. NR >= 1. 00051 * 00052 * SQRE (input) INTEGER 00053 * = 0: the lower block is an NR-by-NR square matrix. 00054 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00055 * 00056 * The bidiagonal matrix has 00057 * N = NL + NR + 1 rows and 00058 * M = N + SQRE >= N columns. 00059 * 00060 * K (output) INTEGER 00061 * Contains the dimension of the non-deflated matrix, this is 00062 * the order of the related secular equation. 1 <= K <=N. 00063 * 00064 * D (input/output) DOUBLE PRECISION array, dimension ( N ) 00065 * On entry D contains the singular values of the two submatrices 00066 * to be combined. On exit D contains the trailing (N-K) updated 00067 * singular values (those which were deflated) sorted into 00068 * increasing order. 00069 * 00070 * Z (output) DOUBLE PRECISION array, dimension ( M ) 00071 * On exit Z contains the updating row vector in the secular 00072 * equation. 00073 * 00074 * ZW (workspace) DOUBLE PRECISION array, dimension ( M ) 00075 * Workspace for Z. 00076 * 00077 * VF (input/output) DOUBLE PRECISION array, dimension ( M ) 00078 * On entry, VF(1:NL+1) contains the first components of all 00079 * right singular vectors of the upper block; and VF(NL+2:M) 00080 * contains the first components of all right singular vectors 00081 * of the lower block. On exit, VF contains the first components 00082 * of all right singular vectors of the bidiagonal matrix. 00083 * 00084 * VFW (workspace) DOUBLE PRECISION array, dimension ( M ) 00085 * Workspace for VF. 00086 * 00087 * VL (input/output) DOUBLE PRECISION array, dimension ( M ) 00088 * On entry, VL(1:NL+1) contains the last components of all 00089 * right singular vectors of the upper block; and VL(NL+2:M) 00090 * contains the last components of all right singular vectors 00091 * of the lower block. On exit, VL contains the last components 00092 * of all right singular vectors of the bidiagonal matrix. 00093 * 00094 * VLW (workspace) DOUBLE PRECISION array, dimension ( M ) 00095 * Workspace for VL. 00096 * 00097 * ALPHA (input) DOUBLE PRECISION 00098 * Contains the diagonal element associated with the added row. 00099 * 00100 * BETA (input) DOUBLE PRECISION 00101 * Contains the off-diagonal element associated with the added 00102 * row. 00103 * 00104 * DSIGMA (output) DOUBLE PRECISION array, dimension ( N ) 00105 * Contains a copy of the diagonal elements (K-1 singular values 00106 * and one zero) in the secular equation. 00107 * 00108 * IDX (workspace) INTEGER array, dimension ( N ) 00109 * This will contain the permutation used to sort the contents of 00110 * D into ascending order. 00111 * 00112 * IDXP (workspace) INTEGER array, dimension ( N ) 00113 * This will contain the permutation used to place deflated 00114 * values of D at the end of the array. On output IDXP(2:K) 00115 * points to the nondeflated D-values and IDXP(K+1:N) 00116 * points to the deflated singular values. 00117 * 00118 * IDXQ (input) INTEGER array, dimension ( N ) 00119 * This contains the permutation which separately sorts the two 00120 * sub-problems in D into ascending order. Note that entries in 00121 * the first half of this permutation must first be moved one 00122 * position backward; and entries in the second half 00123 * must first have NL+1 added to their values. 00124 * 00125 * PERM (output) INTEGER array, dimension ( N ) 00126 * The permutations (from deflation and sorting) to be applied 00127 * to each singular block. Not referenced if ICOMPQ = 0. 00128 * 00129 * GIVPTR (output) INTEGER 00130 * The number of Givens rotations which took place in this 00131 * subproblem. Not referenced if ICOMPQ = 0. 00132 * 00133 * GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) 00134 * Each pair of numbers indicates a pair of columns to take place 00135 * in a Givens rotation. Not referenced if ICOMPQ = 0. 00136 * 00137 * LDGCOL (input) INTEGER 00138 * The leading dimension of GIVCOL, must be at least N. 00139 * 00140 * GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00141 * Each number indicates the C or S value to be used in the 00142 * corresponding Givens rotation. Not referenced if ICOMPQ = 0. 00143 * 00144 * LDGNUM (input) INTEGER 00145 * The leading dimension of GIVNUM, must be at least N. 00146 * 00147 * C (output) DOUBLE PRECISION 00148 * C contains garbage if SQRE =0 and the C-value of a Givens 00149 * rotation related to the right null space if SQRE = 1. 00150 * 00151 * S (output) DOUBLE PRECISION 00152 * S contains garbage if SQRE =0 and the S-value of a Givens 00153 * rotation related to the right null space if SQRE = 1. 00154 * 00155 * INFO (output) INTEGER 00156 * = 0: successful exit. 00157 * < 0: if INFO = -i, the i-th argument had an illegal value. 00158 * 00159 * Further Details 00160 * =============== 00161 * 00162 * Based on contributions by 00163 * Ming Gu and Huan Ren, Computer Science Division, University of 00164 * California at Berkeley, USA 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT 00170 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, 00171 $ EIGHT = 8.0D+0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 * 00175 INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N, 00176 $ NLP1, NLP2 00177 DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL DCOPY, DLAMRG, DROT, XERBLA 00181 * .. 00182 * .. External Functions .. 00183 DOUBLE PRECISION DLAMCH, DLAPY2 00184 EXTERNAL DLAMCH, DLAPY2 00185 * .. 00186 * .. Intrinsic Functions .. 00187 INTRINSIC ABS, MAX 00188 * .. 00189 * .. Executable Statements .. 00190 * 00191 * Test the input parameters. 00192 * 00193 INFO = 0 00194 N = NL + NR + 1 00195 M = N + SQRE 00196 * 00197 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00198 INFO = -1 00199 ELSE IF( NL.LT.1 ) THEN 00200 INFO = -2 00201 ELSE IF( NR.LT.1 ) THEN 00202 INFO = -3 00203 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00204 INFO = -4 00205 ELSE IF( LDGCOL.LT.N ) THEN 00206 INFO = -22 00207 ELSE IF( LDGNUM.LT.N ) THEN 00208 INFO = -24 00209 END IF 00210 IF( INFO.NE.0 ) THEN 00211 CALL XERBLA( 'DLASD7', -INFO ) 00212 RETURN 00213 END IF 00214 * 00215 NLP1 = NL + 1 00216 NLP2 = NL + 2 00217 IF( ICOMPQ.EQ.1 ) THEN 00218 GIVPTR = 0 00219 END IF 00220 * 00221 * Generate the first part of the vector Z and move the singular 00222 * values in the first part of D one position backward. 00223 * 00224 Z1 = ALPHA*VL( NLP1 ) 00225 VL( NLP1 ) = ZERO 00226 TAU = VF( NLP1 ) 00227 DO 10 I = NL, 1, -1 00228 Z( I+1 ) = ALPHA*VL( I ) 00229 VL( I ) = ZERO 00230 VF( I+1 ) = VF( I ) 00231 D( I+1 ) = D( I ) 00232 IDXQ( I+1 ) = IDXQ( I ) + 1 00233 10 CONTINUE 00234 VF( 1 ) = TAU 00235 * 00236 * Generate the second part of the vector Z. 00237 * 00238 DO 20 I = NLP2, M 00239 Z( I ) = BETA*VF( I ) 00240 VF( I ) = ZERO 00241 20 CONTINUE 00242 * 00243 * Sort the singular values into increasing order 00244 * 00245 DO 30 I = NLP2, N 00246 IDXQ( I ) = IDXQ( I ) + NLP1 00247 30 CONTINUE 00248 * 00249 * DSIGMA, IDXC, IDXC, and ZW are used as storage space. 00250 * 00251 DO 40 I = 2, N 00252 DSIGMA( I ) = D( IDXQ( I ) ) 00253 ZW( I ) = Z( IDXQ( I ) ) 00254 VFW( I ) = VF( IDXQ( I ) ) 00255 VLW( I ) = VL( IDXQ( I ) ) 00256 40 CONTINUE 00257 * 00258 CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) ) 00259 * 00260 DO 50 I = 2, N 00261 IDXI = 1 + IDX( I ) 00262 D( I ) = DSIGMA( IDXI ) 00263 Z( I ) = ZW( IDXI ) 00264 VF( I ) = VFW( IDXI ) 00265 VL( I ) = VLW( IDXI ) 00266 50 CONTINUE 00267 * 00268 * Calculate the allowable deflation tolerence 00269 * 00270 EPS = DLAMCH( 'Epsilon' ) 00271 TOL = MAX( ABS( ALPHA ), ABS( BETA ) ) 00272 TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL ) 00273 * 00274 * There are 2 kinds of deflation -- first a value in the z-vector 00275 * is small, second two (or more) singular values are very close 00276 * together (their difference is small). 00277 * 00278 * If the value in the z-vector is small, we simply permute the 00279 * array so that the corresponding singular value is moved to the 00280 * end. 00281 * 00282 * If two values in the D-vector are close, we perform a two-sided 00283 * rotation designed to make one of the corresponding z-vector 00284 * entries zero, and then permute the array so that the deflated 00285 * singular value is moved to the end. 00286 * 00287 * If there are multiple singular values then the problem deflates. 00288 * Here the number of equal singular values are found. As each equal 00289 * singular value is found, an elementary reflector is computed to 00290 * rotate the corresponding singular subspace so that the 00291 * corresponding components of Z are zero in this new basis. 00292 * 00293 K = 1 00294 K2 = N + 1 00295 DO 60 J = 2, N 00296 IF( ABS( Z( J ) ).LE.TOL ) THEN 00297 * 00298 * Deflate due to small z component. 00299 * 00300 K2 = K2 - 1 00301 IDXP( K2 ) = J 00302 IF( J.EQ.N ) 00303 $ GO TO 100 00304 ELSE 00305 JPREV = J 00306 GO TO 70 00307 END IF 00308 60 CONTINUE 00309 70 CONTINUE 00310 J = JPREV 00311 80 CONTINUE 00312 J = J + 1 00313 IF( J.GT.N ) 00314 $ GO TO 90 00315 IF( ABS( Z( J ) ).LE.TOL ) THEN 00316 * 00317 * Deflate due to small z component. 00318 * 00319 K2 = K2 - 1 00320 IDXP( K2 ) = J 00321 ELSE 00322 * 00323 * Check if singular values are close enough to allow deflation. 00324 * 00325 IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN 00326 * 00327 * Deflation is possible. 00328 * 00329 S = Z( JPREV ) 00330 C = Z( J ) 00331 * 00332 * Find sqrt(a**2+b**2) without overflow or 00333 * destructive underflow. 00334 * 00335 TAU = DLAPY2( C, S ) 00336 Z( J ) = TAU 00337 Z( JPREV ) = ZERO 00338 C = C / TAU 00339 S = -S / TAU 00340 * 00341 * Record the appropriate Givens rotation 00342 * 00343 IF( ICOMPQ.EQ.1 ) THEN 00344 GIVPTR = GIVPTR + 1 00345 IDXJP = IDXQ( IDX( JPREV )+1 ) 00346 IDXJ = IDXQ( IDX( J )+1 ) 00347 IF( IDXJP.LE.NLP1 ) THEN 00348 IDXJP = IDXJP - 1 00349 END IF 00350 IF( IDXJ.LE.NLP1 ) THEN 00351 IDXJ = IDXJ - 1 00352 END IF 00353 GIVCOL( GIVPTR, 2 ) = IDXJP 00354 GIVCOL( GIVPTR, 1 ) = IDXJ 00355 GIVNUM( GIVPTR, 2 ) = C 00356 GIVNUM( GIVPTR, 1 ) = S 00357 END IF 00358 CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S ) 00359 CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S ) 00360 K2 = K2 - 1 00361 IDXP( K2 ) = JPREV 00362 JPREV = J 00363 ELSE 00364 K = K + 1 00365 ZW( K ) = Z( JPREV ) 00366 DSIGMA( K ) = D( JPREV ) 00367 IDXP( K ) = JPREV 00368 JPREV = J 00369 END IF 00370 END IF 00371 GO TO 80 00372 90 CONTINUE 00373 * 00374 * Record the last singular value. 00375 * 00376 K = K + 1 00377 ZW( K ) = Z( JPREV ) 00378 DSIGMA( K ) = D( JPREV ) 00379 IDXP( K ) = JPREV 00380 * 00381 100 CONTINUE 00382 * 00383 * Sort the singular values into DSIGMA. The singular values which 00384 * were not deflated go into the first K slots of DSIGMA, except 00385 * that DSIGMA(1) is treated separately. 00386 * 00387 DO 110 J = 2, N 00388 JP = IDXP( J ) 00389 DSIGMA( J ) = D( JP ) 00390 VFW( J ) = VF( JP ) 00391 VLW( J ) = VL( JP ) 00392 110 CONTINUE 00393 IF( ICOMPQ.EQ.1 ) THEN 00394 DO 120 J = 2, N 00395 JP = IDXP( J ) 00396 PERM( J ) = IDXQ( IDX( JP )+1 ) 00397 IF( PERM( J ).LE.NLP1 ) THEN 00398 PERM( J ) = PERM( J ) - 1 00399 END IF 00400 120 CONTINUE 00401 END IF 00402 * 00403 * The deflated singular values go back into the last N - K slots of 00404 * D. 00405 * 00406 CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 ) 00407 * 00408 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and 00409 * VL(M). 00410 * 00411 DSIGMA( 1 ) = ZERO 00412 HLFTOL = TOL / TWO 00413 IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL ) 00414 $ DSIGMA( 2 ) = HLFTOL 00415 IF( M.GT.N ) THEN 00416 Z( 1 ) = DLAPY2( Z1, Z( M ) ) 00417 IF( Z( 1 ).LE.TOL ) THEN 00418 C = ONE 00419 S = ZERO 00420 Z( 1 ) = TOL 00421 ELSE 00422 C = Z1 / Z( 1 ) 00423 S = -Z( M ) / Z( 1 ) 00424 END IF 00425 CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S ) 00426 CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S ) 00427 ELSE 00428 IF( ABS( Z1 ).LE.TOL ) THEN 00429 Z( 1 ) = TOL 00430 ELSE 00431 Z( 1 ) = Z1 00432 END IF 00433 END IF 00434 * 00435 * Restore Z, VF, and VL. 00436 * 00437 CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 ) 00438 CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 ) 00439 CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 ) 00440 * 00441 RETURN 00442 * 00443 * End of DLASD7 00444 * 00445 END
1.7.3 
