1 SUBROUTINE SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, 2 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 3 $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO ) 4* 5* -- LAPACK routine (instrumented to count ops, version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* December 22, 1999 9* 10* .. Scalar Arguments .. 11 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL, 12 $ LDGNUM, NL, NR, NRHS, SQRE 13 REAL C, S 14* .. 15* .. Array Arguments .. 16 INTEGER GIVCOL( LDGCOL, * ), PERM( * ) 17 REAL B( LDB, * ), BX( LDBX, * ), DIFL( * ), 18 $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ), 19 $ POLES( LDGNUM, * ), WORK( * ), Z( * ) 20* .. 21* .. Common block to return operation count .. 22 COMMON / LATIME / OPS, ITCNT 23* .. 24* .. Scalars in Common .. 25 REAL ITCNT, OPS 26* .. 27* 28* Purpose 29* ======= 30* 31* SLALS0 applies back the multiplying factors of either the left or the 32* right singular vector matrix of a diagonal matrix appended by a row 33* to the right hand side matrix B in solving the least squares problem 34* using the divide-and-conquer SVD approach. 35* 36* For the left singular vector matrix, three types of orthogonal 37* matrices are involved: 38* 39* (1L) Givens rotations: the number of such rotations is GIVPTR; the 40* pairs of columns/rows they were applied to are stored in GIVCOL; 41* and the C- and S-values of these rotations are stored in GIVNUM. 42* 43* (2L) Permutation. The (NL+1)-st row of B is to be moved to the first 44* row, and for J=2:N, PERM(J)-th row of B is to be moved to the 45* J-th row. 46* 47* (3L) The left singular vector matrix of the remaining matrix. 48* 49* For the right singular vector matrix, four types of orthogonal 50* matrices are involved: 51* 52* (1R) The right singular vector matrix of the remaining matrix. 53* 54* (2R) If SQRE = 1, one extra Givens rotation to generate the right 55* null space. 56* 57* (3R) The inverse transformation of (2L). 58* 59* (4R) The inverse transformation of (1L). 60* 61* Arguments 62* ========= 63* 64* ICOMPQ (input) INTEGER 65* Specifies whether singular vectors are to be computed in 66* factored form: 67* = 0: Left singular vector matrix. 68* = 1: Right singular vector matrix. 69* 70* NL (input) INTEGER 71* The row dimension of the upper block. NL >= 1. 72* 73* NR (input) INTEGER 74* The row dimension of the lower block. NR >= 1. 75* 76* SQRE (input) INTEGER 77* = 0: the lower block is an NR-by-NR square matrix. 78* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 79* 80* The bidiagonal matrix has row dimension N = NL + NR + 1, 81* and column dimension M = N + SQRE. 82* 83* NRHS (input) INTEGER 84* The number of columns of B and BX. NRHS must be at least 1. 85* 86* B (input/output) REAL array, dimension ( LDB, NRHS ) 87* On input, B contains the right hand sides of the least 88* squares problem in rows 1 through M. On output, B contains 89* the solution X in rows 1 through N. 90* 91* LDB (input) INTEGER 92* The leading dimension of B. LDB must be at least 93* max(1,MAX( M, N ) ). 94* 95* BX (workspace) REAL array, dimension ( LDBX, NRHS ) 96* 97* LDBX (input) INTEGER 98* The leading dimension of BX. 99* 100* PERM (input) INTEGER array, dimension ( N ) 101* The permutations (from deflation and sorting) applied 102* to the two blocks. 103* 104* GIVPTR (input) INTEGER 105* The number of Givens rotations which took place in this 106* subproblem. 107* 108* GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 ) 109* Each pair of numbers indicates a pair of rows/columns 110* involved in a Givens rotation. 111* 112* LDGCOL (input) INTEGER 113* The leading dimension of GIVCOL, must be at least N. 114* 115* GIVNUM (input) REAL array, dimension ( LDGNUM, 2 ) 116* Each number indicates the C or S value used in the 117* corresponding Givens rotation. 118* 119* LDGNUM (input) INTEGER 120* The leading dimension of arrays DIFR, POLES and 121* GIVNUM, must be at least K. 122* 123* POLES (input) REAL array, dimension ( LDGNUM, 2 ) 124* On entry, POLES(1:K, 1) contains the new singular 125* values obtained from solving the secular equation, and 126* POLES(1:K, 2) is an array containing the poles in the secular 127* equation. 128* 129* DIFL (input) REAL array, dimension ( K ). 130* On entry, DIFL(I) is the distance between I-th updated 131* (undeflated) singular value and the I-th (undeflated) old 132* singular value. 133* 134* DIFR (input) REAL array, dimension ( LDGNUM, 2 ). 135* On entry, DIFR(I, 1) contains the distances between I-th 136* updated (undeflated) singular value and the I+1-th 137* (undeflated) old singular value. And DIFR(I, 2) is the 138* normalizing factor for the I-th right singular vector. 139* 140* Z (input) REAL array, dimension ( K ) 141* Contain the components of the deflation-adjusted updating row 142* vector. 143* 144* K (input) INTEGER 145* Contains the dimension of the non-deflated matrix, 146* This is the order of the related secular equation. 1 <= K <=N. 147* 148* C (input) REAL 149* C contains garbage if SQRE =0 and the C-value of a Givens 150* rotation related to the right null space if SQRE = 1. 151* 152* S (input) REAL 153* S contains garbage if SQRE =0 and the S-value of a Givens 154* rotation related to the right null space if SQRE = 1. 155* 156* WORK (workspace) REAL array, dimension ( K ) 157* 158* INFO (output) INTEGER 159* = 0: successful exit. 160* < 0: if INFO = -i, the i-th argument had an illegal value. 161* 162* ===================================================================== 163* 164* .. Parameters .. 165 REAL ONE, ZERO, NEGONE 166 PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 ) 167* .. 168* .. Local Scalars .. 169 INTEGER I, J, M, N, NLP1 170 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP 171* .. 172* .. External Subroutines .. 173 EXTERNAL SCOPY, SGEMV, SLACPY, SLASCL, SROT, SSCAL, 174 $ XERBLA 175* .. 176* .. External Functions .. 177 REAL SLAMC3, SNRM2, SOPBL2 178 EXTERNAL SLAMC3, SNRM2, SOPBL2 179* .. 180* .. Intrinsic Functions .. 181 INTRINSIC REAL 182* .. 183* .. Executable Statements .. 184* 185* Test the input parameters. 186* 187 INFO = 0 188* 189 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 190 INFO = -1 191 ELSE IF( NL.LT.1 ) THEN 192 INFO = -2 193 ELSE IF( NR.LT.1 ) THEN 194 INFO = -3 195 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 196 INFO = -4 197 END IF 198* 199 N = NL + NR + 1 200* 201 IF( NRHS.LT.1 ) THEN 202 INFO = -5 203 ELSE IF( LDB.LT.N ) THEN 204 INFO = -7 205 ELSE IF( LDBX.LT.N ) THEN 206 INFO = -9 207 ELSE IF( GIVPTR.LT.0 ) THEN 208 INFO = -11 209 ELSE IF( LDGCOL.LT.N ) THEN 210 INFO = -13 211 ELSE IF( LDGNUM.LT.N ) THEN 212 INFO = -15 213 ELSE IF( K.LT.1 ) THEN 214 INFO = -20 215 END IF 216 IF( INFO.NE.0 ) THEN 217 CALL XERBLA( 'SLALS0', -INFO ) 218 RETURN 219 END IF 220* 221 M = N + SQRE 222 NLP1 = NL + 1 223* 224 IF( ICOMPQ.EQ.0 ) THEN 225* 226* Apply back orthogonal transformations from the left. 227* 228* Step (1L): apply back the Givens rotations performed. 229* 230 OPS = OPS + REAL( 6*NRHS*GIVPTR ) 231 DO 10 I = 1, GIVPTR 232 CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 233 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 234 $ GIVNUM( I, 1 ) ) 235 10 CONTINUE 236* 237* Step (2L): permute rows of B. 238* 239 CALL SCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) 240 DO 20 I = 2, N 241 CALL SCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX ) 242 20 CONTINUE 243* 244* Step (3L): apply the inverse of the left singular vector 245* matrix to BX. 246* 247 IF( K.EQ.1 ) THEN 248 CALL SCOPY( NRHS, BX, LDBX, B, LDB ) 249 IF( Z( 1 ).LT.ZERO ) THEN 250 OPS = OPS + REAL( NRHS ) 251 CALL SSCAL( NRHS, NEGONE, B, LDB ) 252 END IF 253 ELSE 254 DO 50 J = 1, K 255 DIFLJ = DIFL( J ) 256 DJ = POLES( J, 1 ) 257 DSIGJ = -POLES( J, 2 ) 258 IF( J.LT.K ) THEN 259 DIFRJ = -DIFR( J, 1 ) 260 DSIGJP = -POLES( J+1, 2 ) 261 END IF 262 IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) ) 263 $ THEN 264 WORK( J ) = ZERO 265 ELSE 266 OPS = OPS + REAL( 4 ) 267 WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ / 268 $ ( POLES( J, 2 )+DJ ) 269 END IF 270 DO 30 I = 1, J - 1 271 IF( ( Z( I ).EQ.ZERO ) .OR. 272 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 273 WORK( I ) = ZERO 274 ELSE 275 OPS = OPS + REAL( 6 ) 276 WORK( I ) = POLES( I, 2 )*Z( I ) / 277 $ ( SLAMC3( POLES( I, 2 ), DSIGJ )- 278 $ DIFLJ ) / ( POLES( I, 2 )+DJ ) 279 END IF 280 30 CONTINUE 281 DO 40 I = J + 1, K 282 IF( ( Z( I ).EQ.ZERO ) .OR. 283 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 284 WORK( I ) = ZERO 285 ELSE 286 OPS = OPS + REAL( 6 ) 287 WORK( I ) = POLES( I, 2 )*Z( I ) / 288 $ ( SLAMC3( POLES( I, 2 ), DSIGJP )+ 289 $ DIFRJ ) / ( POLES( I, 2 )+DJ ) 290 END IF 291 40 CONTINUE 292 WORK( 1 ) = NEGONE 293 OPS = OPS + 2*K + NRHS + 294 $ SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) 295 TEMP = SNRM2( K, WORK, 1 ) 296 CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO, 297 $ B( J, 1 ), LDB ) 298 CALL SLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ), 299 $ LDB, INFO ) 300 50 CONTINUE 301 END IF 302* 303* Move the deflated rows of BX to B also. 304* 305 IF( K.LT.MAX( M, N ) ) 306 $ CALL SLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX, 307 $ B( K+1, 1 ), LDB ) 308 ELSE 309* 310* Apply back the right orthogonal transformations. 311* 312* Step (1R): apply back the new right singular vector matrix 313* to B. 314* 315 IF( K.EQ.1 ) THEN 316 CALL SCOPY( NRHS, B, LDB, BX, LDBX ) 317 ELSE 318 DO 80 J = 1, K 319 DSIGJ = POLES( J, 2 ) 320 IF( Z( J ).EQ.ZERO ) THEN 321 WORK( J ) = ZERO 322 ELSE 323 OPS = OPS + REAL( 4 ) 324 WORK( J ) = -Z( J ) / DIFL( J ) / 325 $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 ) 326 END IF 327 DO 60 I = 1, J - 1 328 IF( Z( J ).EQ.ZERO ) THEN 329 WORK( I ) = ZERO 330 ELSE 331 OPS = OPS + REAL( 6 ) 332 WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1, 333 $ 2 ) )-DIFR( I, 1 ) ) / 334 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 335 END IF 336 60 CONTINUE 337 DO 70 I = J + 1, K 338 IF( Z( J ).EQ.ZERO ) THEN 339 WORK( I ) = ZERO 340 ELSE 341 OPS = OPS + REAL( 6 ) 342 WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I, 343 $ 2 ) )-DIFL( I ) ) / 344 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 345 END IF 346 70 CONTINUE 347 OPS = OPS + SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) 348 CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO, 349 $ BX( J, 1 ), LDBX ) 350 80 CONTINUE 351 END IF 352* 353* Step (2R): if SQRE = 1, apply back the rotation that is 354* related to the right null space of the subproblem. 355* 356 IF( SQRE.EQ.1 ) THEN 357 OPS = OPS + REAL( 6*NRHS ) 358 CALL SCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) 359 CALL SROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) 360 END IF 361 IF( K.LT.MAX( M, N ) ) 362 $ CALL SLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, 363 $ BX( K+1, 1 ), LDBX ) 364* 365* Step (3R): permute rows of B. 366* 367 CALL SCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) 368 IF( SQRE.EQ.1 ) THEN 369 CALL SCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) 370 END IF 371 DO 90 I = 2, N 372 CALL SCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) 373 90 CONTINUE 374* 375* Step (4R): apply back the Givens rotations performed. 376* 377 OPS = OPS + REAL( 6*NRHS*GIVPTR ) 378 DO 100 I = GIVPTR, 1, -1 379 CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 380 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 381 $ -GIVNUM( I, 1 ) ) 382 100 CONTINUE 383 END IF 384* 385 RETURN 386* 387* End of SLALS0 388* 389 END 390