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