|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 |
SUBROUTINE SLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
$ LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, $ INFO ) * * -- LAPACK auxiliary routine (version 3.2.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * June 2010 * * .. Scalar Arguments .. INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR, $ SQRE * .. * .. Array Arguments .. INTEGER CTOT( * ), IDXC( * ) REAL D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ), $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ), $ Z( * ) * .. * * Purpose * ======= * * SLASD3 finds all the square roots of the roots of the secular * equation, as defined by the values in D and Z. It makes the * appropriate calls to SLASD4 and then updates the singular * vectors by matrix multiplication. * * This code makes very mild assumptions about floating point * arithmetic. It will work on machines with a guard digit in * add/subtract, or on those binary machines without guard digits * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. * It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * SLASD3 is called from SLASD1. * * Arguments * ========= * * NL (input) INTEGER * The row dimension of the upper block. NL >= 1. * * NR (input) INTEGER * The row dimension of the lower block. NR >= 1. * * SQRE (input) INTEGER * = 0: the lower block is an NR-by-NR square matrix. * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. * * The bidiagonal matrix has N = NL + NR + 1 rows and * M = N + SQRE >= N columns. * * K (input) INTEGER * The size of the secular equation, 1 =< K = < N. * * D (output) REAL array, dimension(K) * On exit the square roots of the roots of the secular equation, * in ascending order. * * Q (workspace) REAL array, * dimension at least (LDQ,K). * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= K. * * DSIGMA (input/output) REAL array, dimension(K) * The first K elements of this array contain the old roots * of the deflated updating problem. These are the poles * of the secular equation. * * U (output) REAL array, dimension (LDU, N) * The last N - K columns of this matrix contain the deflated * left singular vectors. * * LDU (input) INTEGER * The leading dimension of the array U. LDU >= N. * * U2 (input) REAL array, dimension (LDU2, N) * The first K columns of this matrix contain the non-deflated * left singular vectors for the split problem. * * LDU2 (input) INTEGER * The leading dimension of the array U2. LDU2 >= N. * * VT (output) REAL array, dimension (LDVT, M) * The last M - K columns of VT**T contain the deflated * right singular vectors. * * LDVT (input) INTEGER * The leading dimension of the array VT. LDVT >= N. * * VT2 (input/output) REAL array, dimension (LDVT2, N) * The first K columns of VT2**T contain the non-deflated * right singular vectors for the split problem. * * LDVT2 (input) INTEGER * The leading dimension of the array VT2. LDVT2 >= N. * * IDXC (input) INTEGER array, dimension (N) * The permutation used to arrange the columns of U (and rows of * VT) into three groups: the first group contains non-zero * entries only at and above (or before) NL +1; the second * contains non-zero entries only at and below (or after) NL+2; * and the third is dense. The first column of U and the row of * VT are treated separately, however. * * The rows of the singular vectors found by SLASD4 * must be likewise permuted before the matrix multiplies can * take place. * * CTOT (input) INTEGER array, dimension (4) * A count of the total number of the various types of columns * in U (or rows in VT), as described in IDXC. The fourth column * type is any column which has been deflated. * * Z (input/output) REAL array, dimension (K) * The first K elements of this array contain the components * of the deflation-adjusted updating row vector. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, a singular value did not converge * * Further Details * =============== * * Based on contributions by * Ming Gu and Huan Ren, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO, NEGONE PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, $ NEGONE = -1.0E+0 ) * .. * .. Local Scalars .. INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1 REAL RHO, TEMP * .. * .. External Functions .. REAL SLAMC3, SNRM2 EXTERNAL SLAMC3, SNRM2 * .. * .. External Subroutines .. EXTERNAL SCOPY, SGEMM, SLACPY, SLASCL, SLASD4, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, SIGN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( NL.LT.1 ) THEN INFO = -1 ELSE IF( NR.LT.1 ) THEN INFO = -2 ELSE IF( ( SQRE.NE.1 ) .AND. ( SQRE.NE.0 ) ) THEN INFO = -3 END IF * N = NL + NR + 1 M = N + SQRE NLP1 = NL + 1 NLP2 = NL + 2 * IF( ( K.LT.1 ) .OR. ( K.GT.N ) ) THEN INFO = -4 ELSE IF( LDQ.LT.K ) THEN INFO = -7 ELSE IF( LDU.LT.N ) THEN INFO = -10 ELSE IF( LDU2.LT.N ) THEN INFO = -12 ELSE IF( LDVT.LT.M ) THEN INFO = -14 ELSE IF( LDVT2.LT.M ) THEN INFO = -16 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SLASD3', -INFO ) RETURN END IF * * Quick return if possible * IF( K.EQ.1 ) THEN D( 1 ) = ABS( Z( 1 ) ) CALL SCOPY( M, VT2( 1, 1 ), LDVT2, VT( 1, 1 ), LDVT ) IF( Z( 1 ).GT.ZERO ) THEN CALL SCOPY( N, U2( 1, 1 ), 1, U( 1, 1 ), 1 ) ELSE DO 10 I = 1, N U( I, 1 ) = -U2( I, 1 ) 10 CONTINUE END IF RETURN END IF * * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can * be computed with high relative accuracy (barring over/underflow). * This is a problem on machines without a guard digit in * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), * which on any of these machines zeros out the bottommost * bit of DSIGMA(I) if it is 1; this makes the subsequent * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation * occurs. On binary machines with a guard digit (almost all * machines) it does not change DSIGMA(I) at all. On hexadecimal * and decimal machines with a guard digit, it slightly * changes the bottommost bits of DSIGMA(I). It does not account * for hexadecimal or decimal machines without guard digits * (we know of none). We use a subroutine call to compute * 2*DSIGMA(I) to prevent optimizing compilers from eliminating * this code. * DO 20 I = 1, K DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) 20 CONTINUE * * Keep a copy of Z. * CALL SCOPY( K, Z, 1, Q, 1 ) * * Normalize Z. * RHO = SNRM2( K, Z, 1 ) CALL SLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO ) RHO = RHO*RHO * * Find the new singular values. * DO 30 J = 1, K CALL SLASD4( K, J, DSIGMA, Z, U( 1, J ), RHO, D( J ), $ VT( 1, J ), INFO ) * * If the zero finder fails, the computation is terminated. * IF( INFO.NE.0 ) THEN RETURN END IF 30 CONTINUE * * Compute updated Z. * DO 60 I = 1, K Z( I ) = U( I, K )*VT( I, K ) DO 40 J = 1, I - 1 Z( I ) = Z( I )*( U( I, J )*VT( I, J ) / $ ( DSIGMA( I )-DSIGMA( J ) ) / $ ( DSIGMA( I )+DSIGMA( J ) ) ) 40 CONTINUE DO 50 J = I, K - 1 Z( I ) = Z( I )*( U( I, J )*VT( I, J ) / $ ( DSIGMA( I )-DSIGMA( J+1 ) ) / $ ( DSIGMA( I )+DSIGMA( J+1 ) ) ) 50 CONTINUE Z( I ) = SIGN( SQRT( ABS( Z( I ) ) ), Q( I, 1 ) ) 60 CONTINUE * * Compute left singular vectors of the modified diagonal matrix, * and store related information for the right singular vectors. * DO 90 I = 1, K VT( 1, I ) = Z( 1 ) / U( 1, I ) / VT( 1, I ) U( 1, I ) = NEGONE DO 70 J = 2, K VT( J, I ) = Z( J ) / U( J, I ) / VT( J, I ) U( J, I ) = DSIGMA( J )*VT( J, I ) 70 CONTINUE TEMP = SNRM2( K, U( 1, I ), 1 ) Q( 1, I ) = U( 1, I ) / TEMP DO 80 J = 2, K JC = IDXC( J ) Q( J, I ) = U( JC, I ) / TEMP 80 CONTINUE 90 CONTINUE * * Update the left singular vector matrix. * IF( K.EQ.2 ) THEN CALL SGEMM( 'N', 'N', N, K, K, ONE, U2, LDU2, Q, LDQ, ZERO, U, $ LDU ) GO TO 100 END IF IF( CTOT( 1 ).GT.0 ) THEN CALL SGEMM( 'N', 'N', NL, K, CTOT( 1 ), ONE, U2( 1, 2 ), LDU2, $ Q( 2, 1 ), LDQ, ZERO, U( 1, 1 ), LDU ) IF( CTOT( 3 ).GT.0 ) THEN KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) CALL SGEMM( 'N', 'N', NL, K, CTOT( 3 ), ONE, U2( 1, KTEMP ), $ LDU2, Q( KTEMP, 1 ), LDQ, ONE, U( 1, 1 ), LDU ) END IF ELSE IF( CTOT( 3 ).GT.0 ) THEN KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) CALL SGEMM( 'N', 'N', NL, K, CTOT( 3 ), ONE, U2( 1, KTEMP ), $ LDU2, Q( KTEMP, 1 ), LDQ, ZERO, U( 1, 1 ), LDU ) ELSE CALL SLACPY( 'F', NL, K, U2, LDU2, U, LDU ) END IF CALL SCOPY( K, Q( 1, 1 ), LDQ, U( NLP1, 1 ), LDU ) KTEMP = 2 + CTOT( 1 ) CTEMP = CTOT( 2 ) + CTOT( 3 ) CALL SGEMM( 'N', 'N', NR, K, CTEMP, ONE, U2( NLP2, KTEMP ), LDU2, $ Q( KTEMP, 1 ), LDQ, ZERO, U( NLP2, 1 ), LDU ) * * Generate the right singular vectors. * 100 CONTINUE DO 120 I = 1, K TEMP = SNRM2( K, VT( 1, I ), 1 ) Q( I, 1 ) = VT( 1, I ) / TEMP DO 110 J = 2, K JC = IDXC( J ) Q( I, J ) = VT( JC, I ) / TEMP 110 CONTINUE 120 CONTINUE * * Update the right singular vector matrix. * IF( K.EQ.2 ) THEN CALL SGEMM( 'N', 'N', K, M, K, ONE, Q, LDQ, VT2, LDVT2, ZERO, $ VT, LDVT ) RETURN END IF KTEMP = 1 + CTOT( 1 ) CALL SGEMM( 'N', 'N', K, NLP1, KTEMP, ONE, Q( 1, 1 ), LDQ, $ VT2( 1, 1 ), LDVT2, ZERO, VT( 1, 1 ), LDVT ) KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) IF( KTEMP.LE.LDVT2 ) $ CALL SGEMM( 'N', 'N', K, NLP1, CTOT( 3 ), ONE, Q( 1, KTEMP ), $ LDQ, VT2( KTEMP, 1 ), LDVT2, ONE, VT( 1, 1 ), $ LDVT ) * KTEMP = CTOT( 1 ) + 1 NRP1 = NR + SQRE IF( KTEMP.GT.1 ) THEN DO 130 I = 1, K Q( I, KTEMP ) = Q( I, 1 ) 130 CONTINUE DO 140 I = NLP2, M VT2( KTEMP, I ) = VT2( 1, I ) 140 CONTINUE END IF CTEMP = 1 + CTOT( 2 ) + CTOT( 3 ) CALL SGEMM( 'N', 'N', K, NRP1, CTEMP, ONE, Q( 1, KTEMP ), LDQ, $ VT2( KTEMP, NLP2 ), LDVT2, ZERO, VT( 1, NLP2 ), LDVT ) * RETURN * * End of SLASD3 * END |