1 SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 2 $ RANK, WORK, RWORK, IWORK, INFO ) 3* 4* -- LAPACK routine (instrumented to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* October 31, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER UPLO 11 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ 12 REAL RCOND 13* .. 14* .. Array Arguments .. 15 INTEGER IWORK( * ) 16 REAL D( * ), E( * ), RWORK( * ) 17 COMPLEX B( LDB, * ), WORK( * ) 18* .. 19* .. Common block to return operation count .. 20 COMMON / LATIME / OPS, ITCNT 21* .. 22* .. Scalars in Common .. 23 REAL ITCNT, OPS 24* .. 25* 26* Purpose 27* ======= 28* 29* CLALSD uses the singular value decomposition of A to solve the least 30* squares problem of finding X to minimize the Euclidean norm of each 31* column of A*X-B, where A is N-by-N upper bidiagonal, and X and B 32* are N-by-NRHS. The solution X overwrites B. 33* 34* The singular values of A smaller than RCOND times the largest 35* singular value are treated as zero in solving the least squares 36* problem; in this case a minimum norm solution is returned. 37* The actual singular values are returned in D in ascending order. 38* 39* This code makes very mild assumptions about floating point 40* arithmetic. It will work on machines with a guard digit in 41* add/subtract, or on those binary machines without guard digits 42* which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. 43* It could conceivably fail on hexadecimal or decimal machines 44* without guard digits, but we know of none. 45* 46* Arguments 47* ========= 48* 49* UPLO (input) CHARACTER*1 50* = 'U': D and E define an upper bidiagonal matrix. 51* = 'L': D and E define a lower bidiagonal matrix. 52* 53* SMLSIZ (input) INTEGER 54* The maximum size of the subproblems at the bottom of the 55* computation tree. 56* 57* N (input) INTEGER 58* The dimension of the bidiagonal matrix. N >= 0. 59* 60* NRHS (input) INTEGER 61* The number of columns of B. NRHS must be at least 1. 62* 63* D (input/output) REAL array, dimension (N) 64* On entry D contains the main diagonal of the bidiagonal 65* matrix. On exit, if INFO = 0, D contains its singular values. 66* 67* E (input) REAL array, dimension (N-1) 68* Contains the super-diagonal entries of the bidiagonal matrix. 69* On exit, E has been destroyed. 70* 71* B (input/output) REAL array, dimension (LDB,NRHS) 72* On input, B contains the right hand sides of the least 73* squares problem. On output, B contains the solution X. 74* 75* LDB (input) INTEGER 76* The leading dimension of B in the calling subprogram. 77* LDB must be at least max(1,N). 78* 79* RCOND (input) REAL 80* The singular values of A less than or equal to RCOND times 81* the largest singular value are treated as zero in solving 82* the least squares problem. If RCOND is negative, 83* machine precision is used instead. 84* For example, if diag(S)*X=B were the least squares problem, 85* where diag(S) is a diagonal matrix of singular values, the 86* solution would be X(i) = B(i) / S(i) if S(i) is greater than 87* RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to 88* RCOND*max(S). 89* 90* RANK (output) INTEGER 91* The number of singular values of A greater than RCOND times 92* the largest singular value. 93* 94* WORK (workspace) COMPLEX array, dimension at least 95* (N * NRHS). 96* 97* RWORK (workspace) REAL array, dimension at least 98* (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2), 99* where 100* NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 101* 102* IWORK (workspace) INTEGER array, dimension at least 103* (3*N*NLVL + 11*N). 104* 105* INFO (output) INTEGER 106* = 0: successful exit. 107* < 0: if INFO = -i, the i-th argument had an illegal value. 108* > 0: The algorithm failed to compute an singular value while 109* working on the submatrix lying in rows and columns 110* INFO/(N+1) through MOD(INFO,N+1). 111* 112* ===================================================================== 113* 114* .. Parameters .. 115 REAL ZERO, ONE, TWO 116 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 117 COMPLEX CZERO 118 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ) ) 119* .. 120* .. Local Scalars .. 121 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 122 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB, 123 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG, 124 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB, 125 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1, 126 $ U, VT, Z 127 REAL CS, EPS, ORGNRM, R, SN, TOL 128* .. 129* .. External Subroutines .. 130 EXTERNAL CSROT, CCOPY, CLACPY, 131 $ CLALSA, CLASCL, CLASET, SGEMM, 132 $ SLARTG, SLASCL, SLASDA, SLASDQ, 133 $ SLASET, SLASRT, XERBLA 134* .. 135* .. External Functions .. 136 INTEGER ISAMAX 137 REAL SLAMCH, SLANST, SOPBL3 138 EXTERNAL ISAMAX, SLAMCH, SLANST, SOPBL3 139* .. 140* .. Intrinsic Functions .. 141 INTRINSIC CMPLX, REAL, AIMAG, ABS, INT, LOG, SIGN 142* .. 143* .. Executable Statements .. 144* 145* Test the input parameters. 146* 147 INFO = 0 148* 149 IF( N.LT.0 ) THEN 150 INFO = -3 151 ELSE IF( NRHS.LT.1 ) THEN 152 INFO = -4 153 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 154 INFO = -8 155 END IF 156 IF( INFO.NE.0 ) THEN 157 CALL XERBLA( 'CLALSD', -INFO ) 158 RETURN 159 END IF 160* 161 EPS = SLAMCH( 'Epsilon' ) 162* 163* Set up the tolerance. 164* 165 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 166 RCOND = EPS 167 END IF 168* 169 RANK = 0 170* 171* Quick return if possible. 172* 173 IF( N.EQ.0 ) THEN 174 RETURN 175 ELSE IF( N.EQ.1 ) THEN 176 IF( D( 1 ).EQ.ZERO ) THEN 177 CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB ) 178 ELSE 179 RANK = 1 180 OPS = OPS + REAL( 2*NRHS ) 181 CALL CLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 182 D( 1 ) = ABS( D( 1 ) ) 183 END IF 184 RETURN 185 END IF 186* 187* Rotate the matrix if it is lower bidiagonal. 188* 189 IF( UPLO.EQ.'L' ) THEN 190 OPS = OPS + REAL( 6*( N-1 ) ) 191 DO 10 I = 1, N - 1 192 CALL SLARTG( D( I ), E( I ), CS, SN, R ) 193 D( I ) = R 194 E( I ) = SN*D( I+1 ) 195 D( I+1 ) = CS*D( I+1 ) 196 IF( NRHS.EQ.1 ) THEN 197 OPS = OPS + REAL( 12 ) 198 CALL CSROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 199 ELSE 200 RWORK( I*2-1 ) = CS 201 RWORK( I*2 ) = SN 202 END IF 203 10 CONTINUE 204 IF( NRHS.GT.1 ) THEN 205 OPS = OPS + REAL( 12*( N-1 )*NRHS ) 206 DO 30 I = 1, NRHS 207 DO 20 J = 1, N - 1 208 CS = RWORK( J*2-1 ) 209 SN = RWORK( J*2 ) 210 CALL CSROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 211 20 CONTINUE 212 30 CONTINUE 213 END IF 214 END IF 215* 216* Scale. 217* 218 NM1 = N - 1 219 ORGNRM = SLANST( 'M', N, D, E ) 220 IF( ORGNRM.EQ.ZERO ) THEN 221 CALL CLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB ) 222 RETURN 223 END IF 224* 225 OPS = OPS + REAL( N + NM1 ) 226 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 227 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 228* 229* If N is smaller than the minimum divide size SMLSIZ, then solve 230* the problem with another solver. 231* 232 IF( N.LE.SMLSIZ ) THEN 233 IRWU = 1 234 IRWVT = IRWU + N*N 235 IRWWRK = IRWVT + N*N 236 IRWRB = IRWWRK 237 IRWIB = IRWRB + N*NRHS 238 IRWB = IRWIB + N*NRHS 239 CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N ) 240 CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N ) 241 CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N, 242 $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1, 243 $ RWORK( IRWWRK ), INFO ) 244 IF( INFO.NE.0 ) THEN 245 RETURN 246 END IF 247* 248* In the real version, B is passed to SLASDQ and multiplied 249* internally by Q'. Here B is complex and that product is 250* computed below in two steps (real and imaginary parts). 251* 252 J = IRWB - 1 253 DO 50 JCOL = 1, NRHS 254 DO 40 JROW = 1, N 255 J = J + 1 256 RWORK( J ) = REAL( B( JROW, JCOL ) ) 257 40 CONTINUE 258 50 CONTINUE 259 OPS = OPS + SOPBL3( 'SGEMM ', N, NRHS, N ) 260 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 261 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 262 J = IRWB - 1 263 DO 70 JCOL = 1, NRHS 264 DO 60 JROW = 1, N 265 J = J + 1 266 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 267 60 CONTINUE 268 70 CONTINUE 269 OPS = OPS + SOPBL3( 'SGEMM ', N, NRHS, N ) 270 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 271 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 272 JREAL = IRWRB - 1 273 JIMAG = IRWIB - 1 274 DO 90 JCOL = 1, NRHS 275 DO 80 JROW = 1, N 276 JREAL = JREAL + 1 277 JIMAG = JIMAG + 1 278 B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 279 $ RWORK( JIMAG ) ) 280 80 CONTINUE 281 90 CONTINUE 282* 283 OPS = OPS + REAL( 1 ) 284 TOL = RCOND*ABS( D( ISAMAX( N, D, 1 ) ) ) 285 DO 100 I = 1, N 286 IF( D( I ).LE.TOL ) THEN 287 CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 288 ELSE 289 OPS = OPS + REAL( 6*NRHS ) 290 CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 291 $ LDB, INFO ) 292 RANK = RANK + 1 293 END IF 294 100 CONTINUE 295* 296* Since B is complex, the following call to SGEMM is performed 297* in two steps (real and imaginary parts). That is for V * B 298* (in the real version of the code V' is stored in WORK). 299* 300* CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 301* $ WORK( NWORK ), N ) 302* 303 J = IRWB - 1 304 DO 120 JCOL = 1, NRHS 305 DO 110 JROW = 1, N 306 J = J + 1 307 RWORK( J ) = REAL( B( JROW, JCOL ) ) 308 110 CONTINUE 309 120 CONTINUE 310 OPS = OPS + SOPBL3( 'SGEMM ', N, NRHS, N ) 311 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 312 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 313 J = IRWB - 1 314 DO 140 JCOL = 1, NRHS 315 DO 130 JROW = 1, N 316 J = J + 1 317 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 318 130 CONTINUE 319 140 CONTINUE 320 OPS = OPS + SOPBL3( 'SGEMM ', N, NRHS, N ) 321 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 322 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 323 JREAL = IRWRB - 1 324 JIMAG = IRWIB - 1 325 DO 160 JCOL = 1, NRHS 326 DO 150 JROW = 1, N 327 JREAL = JREAL + 1 328 JIMAG = JIMAG + 1 329 B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 330 $ RWORK( JIMAG ) ) 331 150 CONTINUE 332 160 CONTINUE 333* 334* Unscale. 335* 336 OPS = OPS + REAL( N + 6*N*NRHS ) 337 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 338 CALL SLASRT( 'D', N, D, INFO ) 339 CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 340* 341 RETURN 342 END IF 343* 344* Book-keeping and setting up some constants. 345* 346 NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 347* 348 SMLSZP = SMLSIZ + 1 349* 350 U = 1 351 VT = 1 + SMLSIZ*N 352 DIFL = VT + SMLSZP*N 353 DIFR = DIFL + NLVL*N 354 Z = DIFR + NLVL*N*2 355 C = Z + NLVL*N 356 S = C + N 357 POLES = S + N 358 GIVNUM = POLES + 2*NLVL*N 359 NRWORK = GIVNUM + 2*NLVL*N 360 BX = 1 361* 362 IRWRB = NRWORK 363 IRWIB = IRWRB + SMLSIZ*NRHS 364 IRWB = IRWIB + SMLSIZ*NRHS 365* 366 SIZEI = 1 + N 367 K = SIZEI + N 368 GIVPTR = K + N 369 PERM = GIVPTR + N 370 GIVCOL = PERM + NLVL*N 371 IWK = GIVCOL + NLVL*N*2 372* 373 ST = 1 374 SQRE = 0 375 ICMPQ1 = 1 376 ICMPQ2 = 0 377 NSUB = 0 378* 379 DO 170 I = 1, N 380 IF( ABS( D( I ) ).LT.EPS ) THEN 381 D( I ) = SIGN( EPS, D( I ) ) 382 END IF 383 170 CONTINUE 384* 385 DO 240 I = 1, NM1 386 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 387 NSUB = NSUB + 1 388 IWORK( NSUB ) = ST 389* 390* Subproblem found. First determine its size and then 391* apply divide and conquer on it. 392* 393 IF( I.LT.NM1 ) THEN 394* 395* A subproblem with E(I) small for I < NM1. 396* 397 NSIZE = I - ST + 1 398 IWORK( SIZEI+NSUB-1 ) = NSIZE 399 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 400* 401* A subproblem with E(NM1) not too small but I = NM1. 402* 403 NSIZE = N - ST + 1 404 IWORK( SIZEI+NSUB-1 ) = NSIZE 405 ELSE 406* 407* A subproblem with E(NM1) small. This implies an 408* 1-by-1 subproblem at D(N), which is not solved 409* explicitly. 410* 411 NSIZE = I - ST + 1 412 IWORK( SIZEI+NSUB-1 ) = NSIZE 413 NSUB = NSUB + 1 414 IWORK( NSUB ) = N 415 IWORK( SIZEI+NSUB-1 ) = 1 416 CALL CCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 417 END IF 418 ST1 = ST - 1 419 IF( NSIZE.EQ.1 ) THEN 420* 421* This is a 1-by-1 subproblem and is not solved 422* explicitly. 423* 424 CALL CCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 425 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 426* 427* This is a small subproblem and is solved by SLASDQ. 428* 429 CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 430 $ RWORK( VT+ST1 ), N ) 431 CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 432 $ RWORK( U+ST1 ), N ) 433 CALL SLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ), 434 $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ), 435 $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ), 436 $ INFO ) 437 IF( INFO.NE.0 ) THEN 438 RETURN 439 END IF 440* 441* In the real version, B is passed to SLASDQ and multiplied 442* internally by Q'. Here B is complex and that product is 443* computed below in two steps (real and imaginary parts). 444* 445 J = IRWB - 1 446 DO 190 JCOL = 1, NRHS 447 DO 180 JROW = ST, ST + NSIZE - 1 448 J = J + 1 449 RWORK( J ) = REAL( B( JROW, JCOL ) ) 450 180 CONTINUE 451 190 CONTINUE 452 OPS = OPS + SOPBL3( 'SGEMM ', NSIZE, NRHS, NSIZE ) 453 CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 454 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 455 $ ZERO, RWORK( IRWRB ), NSIZE ) 456 J = IRWB - 1 457 DO 210 JCOL = 1, NRHS 458 DO 200 JROW = ST, ST + NSIZE - 1 459 J = J + 1 460 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 461 200 CONTINUE 462 210 CONTINUE 463 OPS = OPS + SOPBL3( 'SGEMM ', NSIZE, NRHS, NSIZE ) 464 CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 465 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 466 $ ZERO, RWORK( IRWIB ), NSIZE ) 467 JREAL = IRWRB - 1 468 JIMAG = IRWIB - 1 469 DO 230 JCOL = 1, NRHS 470 DO 220 JROW = ST, ST + NSIZE - 1 471 JREAL = JREAL + 1 472 JIMAG = JIMAG + 1 473 B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 474 $ RWORK( JIMAG ) ) 475 220 CONTINUE 476 230 CONTINUE 477* 478 CALL CLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 479 $ WORK( BX+ST1 ), N ) 480 ELSE 481* 482* A large problem. Solve it using divide and conquer. 483* 484 CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 485 $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ), 486 $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ), 487 $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ), 488 $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 489 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 490 $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ), 491 $ RWORK( S+ST1 ), RWORK( NRWORK ), 492 $ IWORK( IWK ), INFO ) 493 IF( INFO.NE.0 ) THEN 494 RETURN 495 END IF 496 BXST = BX + ST1 497 CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 498 $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N, 499 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 500 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 501 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 502 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 503 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 504 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 505 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 506 IF( INFO.NE.0 ) THEN 507 RETURN 508 END IF 509 END IF 510 ST = I + 1 511 END IF 512 240 CONTINUE 513* 514* Apply the singular values and treat the tiny ones as zero. 515* 516 OPS = OPS + REAL( 1 ) 517 TOL = RCOND*ABS( D( ISAMAX( N, D, 1 ) ) ) 518* 519 DO 250 I = 1, N 520* 521* Some of the elements in D can be negative because 1-by-1 522* subproblems were not solved explicitly. 523* 524 IF( ABS( D( I ) ).LE.TOL ) THEN 525 CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N ) 526 ELSE 527 RANK = RANK + 1 528 OPS = OPS + REAL( 6*NRHS ) 529 CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 530 $ WORK( BX+I-1 ), N, INFO ) 531 END IF 532 D( I ) = ABS( D( I ) ) 533 250 CONTINUE 534* 535* Now apply back the right singular vectors. 536* 537 ICMPQ2 = 1 538 DO 320 I = 1, NSUB 539 ST = IWORK( I ) 540 ST1 = ST - 1 541 NSIZE = IWORK( SIZEI+I-1 ) 542 BXST = BX + ST1 543 IF( NSIZE.EQ.1 ) THEN 544 CALL CCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 545 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 546* 547* Since B and BX are complex, the following call to SGEMM 548* is performed in two steps (real and imaginary parts). 549* 550* CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 551* $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO, 552* $ B( ST, 1 ), LDB ) 553* 554 J = BXST - N - 1 555 JREAL = IRWB - 1 556 DO 270 JCOL = 1, NRHS 557 J = J + N 558 DO 260 JROW = 1, NSIZE 559 JREAL = JREAL + 1 560 RWORK( JREAL ) = REAL( WORK( J+JROW ) ) 561 260 CONTINUE 562 270 CONTINUE 563 OPS = OPS + SOPBL3( 'SGEMM ', NSIZE, NRHS, NSIZE ) 564 CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 565 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 566 $ RWORK( IRWRB ), NSIZE ) 567 J = BXST - N - 1 568 JIMAG = IRWB - 1 569 DO 290 JCOL = 1, NRHS 570 J = J + N 571 DO 280 JROW = 1, NSIZE 572 JIMAG = JIMAG + 1 573 RWORK( JIMAG ) = AIMAG( WORK( J+JROW ) ) 574 280 CONTINUE 575 290 CONTINUE 576 OPS = OPS + SOPBL3( 'SGEMM ', NSIZE, NRHS, NSIZE ) 577 CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 578 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 579 $ RWORK( IRWIB ), NSIZE ) 580 JREAL = IRWRB - 1 581 JIMAG = IRWIB - 1 582 DO 310 JCOL = 1, NRHS 583 DO 300 JROW = ST, ST + NSIZE - 1 584 JREAL = JREAL + 1 585 JIMAG = JIMAG + 1 586 B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 587 $ RWORK( JIMAG ) ) 588 300 CONTINUE 589 310 CONTINUE 590 ELSE 591 CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 592 $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N, 593 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 594 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 595 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 596 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 597 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 598 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 599 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 600 IF( INFO.NE.0 ) THEN 601 RETURN 602 END IF 603 END IF 604 320 CONTINUE 605* 606* Unscale and sort the singular values. 607* 608 OPS = OPS + REAL( N + 6*N*NRHS ) 609 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 610 CALL SLASRT( 'D', N, D, INFO ) 611 CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 612* 613 RETURN 614* 615* End of CLALSD 616* 617 END 618