1*> \brief \b SLALSD 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 SLALSD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slalsd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slalsd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slalsd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 22* RANK, WORK, 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 B( LDB, * ), D( * ), E( * ), WORK( * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> SLALSD uses the singular value decomposition of A to solve the least 41*> squares problem of finding X to minimize the Euclidean norm of each 42*> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B 43*> are N-by-NRHS. The solution X overwrites B. 44*> 45*> The singular values of A smaller than RCOND times the largest 46*> singular value are treated as zero in solving the least squares 47*> problem; in this case a minimum norm solution is returned. 48*> The actual singular values are returned in D in ascending order. 49*> 50*> This code makes very mild assumptions about floating point 51*> arithmetic. It will work on machines with a guard digit in 52*> add/subtract, or on those binary machines without guard digits 53*> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. 54*> It could conceivably fail on hexadecimal or decimal machines 55*> without guard digits, but we know of none. 56*> \endverbatim 57* 58* Arguments: 59* ========== 60* 61*> \param[in] UPLO 62*> \verbatim 63*> UPLO is CHARACTER*1 64*> = 'U': D and E define an upper bidiagonal matrix. 65*> = 'L': D and E define a lower bidiagonal matrix. 66*> \endverbatim 67*> 68*> \param[in] SMLSIZ 69*> \verbatim 70*> SMLSIZ is INTEGER 71*> The maximum size of the subproblems at the bottom of the 72*> computation tree. 73*> \endverbatim 74*> 75*> \param[in] N 76*> \verbatim 77*> N is INTEGER 78*> The dimension of the bidiagonal matrix. N >= 0. 79*> \endverbatim 80*> 81*> \param[in] NRHS 82*> \verbatim 83*> NRHS is INTEGER 84*> The number of columns of B. NRHS must be at least 1. 85*> \endverbatim 86*> 87*> \param[in,out] D 88*> \verbatim 89*> D is REAL array, dimension (N) 90*> On entry D contains the main diagonal of the bidiagonal 91*> matrix. On exit, if INFO = 0, D contains its singular values. 92*> \endverbatim 93*> 94*> \param[in,out] E 95*> \verbatim 96*> E is REAL array, dimension (N-1) 97*> Contains the super-diagonal entries of the bidiagonal matrix. 98*> On exit, E has been destroyed. 99*> \endverbatim 100*> 101*> \param[in,out] B 102*> \verbatim 103*> B is REAL array, dimension (LDB,NRHS) 104*> On input, B contains the right hand sides of the least 105*> squares problem. On output, B contains the solution X. 106*> \endverbatim 107*> 108*> \param[in] LDB 109*> \verbatim 110*> LDB is INTEGER 111*> The leading dimension of B in the calling subprogram. 112*> LDB must be at least max(1,N). 113*> \endverbatim 114*> 115*> \param[in] RCOND 116*> \verbatim 117*> RCOND is REAL 118*> The singular values of A less than or equal to RCOND times 119*> the largest singular value are treated as zero in solving 120*> the least squares problem. If RCOND is negative, 121*> machine precision is used instead. 122*> For example, if diag(S)*X=B were the least squares problem, 123*> where diag(S) is a diagonal matrix of singular values, the 124*> solution would be X(i) = B(i) / S(i) if S(i) is greater than 125*> RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to 126*> RCOND*max(S). 127*> \endverbatim 128*> 129*> \param[out] RANK 130*> \verbatim 131*> RANK is INTEGER 132*> The number of singular values of A greater than RCOND times 133*> the largest singular value. 134*> \endverbatim 135*> 136*> \param[out] WORK 137*> \verbatim 138*> WORK is REAL array, dimension at least 139*> (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2), 140*> where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1). 141*> \endverbatim 142*> 143*> \param[out] IWORK 144*> \verbatim 145*> IWORK is INTEGER array, dimension at least 146*> (3*N*NLVL + 11*N) 147*> \endverbatim 148*> 149*> \param[out] INFO 150*> \verbatim 151*> INFO is INTEGER 152*> = 0: successful exit. 153*> < 0: if INFO = -i, the i-th argument had an illegal value. 154*> > 0: The algorithm failed to compute a singular value while 155*> working on the submatrix lying in rows and columns 156*> INFO/(N+1) through MOD(INFO,N+1). 157*> \endverbatim 158* 159* Authors: 160* ======== 161* 162*> \author Univ. of Tennessee 163*> \author Univ. of California Berkeley 164*> \author Univ. of Colorado Denver 165*> \author NAG Ltd. 166* 167*> \date December 2016 168* 169*> \ingroup realOTHERcomputational 170* 171*> \par Contributors: 172* ================== 173*> 174*> Ming Gu and Ren-Cang Li, Computer Science Division, University of 175*> California at Berkeley, USA \n 176*> Osni Marques, LBNL/NERSC, USA \n 177* 178* ===================================================================== 179 SUBROUTINE SLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 180 $ RANK, WORK, IWORK, INFO ) 181* 182* -- LAPACK computational routine (version 3.7.0) -- 183* -- LAPACK is a software package provided by Univ. of Tennessee, -- 184* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 185* December 2016 186* 187* .. Scalar Arguments .. 188 CHARACTER UPLO 189 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ 190 REAL RCOND 191* .. 192* .. Array Arguments .. 193 INTEGER IWORK( * ) 194 REAL B( LDB, * ), D( * ), E( * ), WORK( * ) 195* .. 196* 197* ===================================================================== 198* 199* .. Parameters .. 200 REAL ZERO, ONE, TWO 201 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 202* .. 203* .. Local Scalars .. 204 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 205 $ GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL, 206 $ NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI, 207 $ SMLSZP, SQRE, ST, ST1, U, VT, Z 208 REAL CS, EPS, ORGNRM, R, RCND, SN, TOL 209* .. 210* .. External Functions .. 211 INTEGER ISAMAX 212 REAL SLAMCH, SLANST 213 EXTERNAL ISAMAX, SLAMCH, SLANST 214* .. 215* .. External Subroutines .. 216 EXTERNAL SCOPY, SGEMM, SLACPY, SLALSA, SLARTG, SLASCL, 217 $ SLASDA, SLASDQ, SLASET, SLASRT, SROT, XERBLA 218* .. 219* .. Intrinsic Functions .. 220 INTRINSIC ABS, INT, LOG, REAL, SIGN 221* .. 222* .. Executable Statements .. 223* 224* Test the input parameters. 225* 226 INFO = 0 227* 228 IF( N.LT.0 ) THEN 229 INFO = -3 230 ELSE IF( NRHS.LT.1 ) THEN 231 INFO = -4 232 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 233 INFO = -8 234 END IF 235 IF( INFO.NE.0 ) THEN 236 CALL XERBLA( 'SLALSD', -INFO ) 237 RETURN 238 END IF 239* 240 EPS = SLAMCH( 'Epsilon' ) 241* 242* Set up the tolerance. 243* 244 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 245 RCND = EPS 246 ELSE 247 RCND = RCOND 248 END IF 249* 250 RANK = 0 251* 252* Quick return if possible. 253* 254 IF( N.EQ.0 ) THEN 255 RETURN 256 ELSE IF( N.EQ.1 ) THEN 257 IF( D( 1 ).EQ.ZERO ) THEN 258 CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB ) 259 ELSE 260 RANK = 1 261 CALL SLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 262 D( 1 ) = ABS( D( 1 ) ) 263 END IF 264 RETURN 265 END IF 266* 267* Rotate the matrix if it is lower bidiagonal. 268* 269 IF( UPLO.EQ.'L' ) THEN 270 DO 10 I = 1, N - 1 271 CALL SLARTG( D( I ), E( I ), CS, SN, R ) 272 D( I ) = R 273 E( I ) = SN*D( I+1 ) 274 D( I+1 ) = CS*D( I+1 ) 275 IF( NRHS.EQ.1 ) THEN 276 CALL SROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 277 ELSE 278 WORK( I*2-1 ) = CS 279 WORK( I*2 ) = SN 280 END IF 281 10 CONTINUE 282 IF( NRHS.GT.1 ) THEN 283 DO 30 I = 1, NRHS 284 DO 20 J = 1, N - 1 285 CS = WORK( J*2-1 ) 286 SN = WORK( J*2 ) 287 CALL SROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 288 20 CONTINUE 289 30 CONTINUE 290 END IF 291 END IF 292* 293* Scale. 294* 295 NM1 = N - 1 296 ORGNRM = SLANST( 'M', N, D, E ) 297 IF( ORGNRM.EQ.ZERO ) THEN 298 CALL SLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB ) 299 RETURN 300 END IF 301* 302 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 303 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 304* 305* If N is smaller than the minimum divide size SMLSIZ, then solve 306* the problem with another solver. 307* 308 IF( N.LE.SMLSIZ ) THEN 309 NWORK = 1 + N*N 310 CALL SLASET( 'A', N, N, ZERO, ONE, WORK, N ) 311 CALL SLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B, 312 $ LDB, WORK( NWORK ), INFO ) 313 IF( INFO.NE.0 ) THEN 314 RETURN 315 END IF 316 TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) ) 317 DO 40 I = 1, N 318 IF( D( I ).LE.TOL ) THEN 319 CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 320 ELSE 321 CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 322 $ LDB, INFO ) 323 RANK = RANK + 1 324 END IF 325 40 CONTINUE 326 CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 327 $ WORK( NWORK ), N ) 328 CALL SLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB ) 329* 330* Unscale. 331* 332 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 333 CALL SLASRT( 'D', N, D, INFO ) 334 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 335* 336 RETURN 337 END IF 338* 339* Book-keeping and setting up some constants. 340* 341 NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 342* 343 SMLSZP = SMLSIZ + 1 344* 345 U = 1 346 VT = 1 + SMLSIZ*N 347 DIFL = VT + SMLSZP*N 348 DIFR = DIFL + NLVL*N 349 Z = DIFR + NLVL*N*2 350 C = Z + NLVL*N 351 S = C + N 352 POLES = S + N 353 GIVNUM = POLES + 2*NLVL*N 354 BX = GIVNUM + 2*NLVL*N 355 NWORK = BX + N*NRHS 356* 357 SIZEI = 1 + N 358 K = SIZEI + N 359 GIVPTR = K + N 360 PERM = GIVPTR + N 361 GIVCOL = PERM + NLVL*N 362 IWK = GIVCOL + NLVL*N*2 363* 364 ST = 1 365 SQRE = 0 366 ICMPQ1 = 1 367 ICMPQ2 = 0 368 NSUB = 0 369* 370 DO 50 I = 1, N 371 IF( ABS( D( I ) ).LT.EPS ) THEN 372 D( I ) = SIGN( EPS, D( I ) ) 373 END IF 374 50 CONTINUE 375* 376 DO 60 I = 1, NM1 377 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 378 NSUB = NSUB + 1 379 IWORK( NSUB ) = ST 380* 381* Subproblem found. First determine its size and then 382* apply divide and conquer on it. 383* 384 IF( I.LT.NM1 ) THEN 385* 386* A subproblem with E(I) small for I < NM1. 387* 388 NSIZE = I - ST + 1 389 IWORK( SIZEI+NSUB-1 ) = NSIZE 390 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 391* 392* A subproblem with E(NM1) not too small but I = NM1. 393* 394 NSIZE = N - ST + 1 395 IWORK( SIZEI+NSUB-1 ) = NSIZE 396 ELSE 397* 398* A subproblem with E(NM1) small. This implies an 399* 1-by-1 subproblem at D(N), which is not solved 400* explicitly. 401* 402 NSIZE = I - ST + 1 403 IWORK( SIZEI+NSUB-1 ) = NSIZE 404 NSUB = NSUB + 1 405 IWORK( NSUB ) = N 406 IWORK( SIZEI+NSUB-1 ) = 1 407 CALL SCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 408 END IF 409 ST1 = ST - 1 410 IF( NSIZE.EQ.1 ) THEN 411* 412* This is a 1-by-1 subproblem and is not solved 413* explicitly. 414* 415 CALL SCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 416 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 417* 418* This is a small subproblem and is solved by SLASDQ. 419* 420 CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 421 $ WORK( VT+ST1 ), N ) 422 CALL SLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ), 423 $ E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ), 424 $ N, B( ST, 1 ), LDB, WORK( NWORK ), INFO ) 425 IF( INFO.NE.0 ) THEN 426 RETURN 427 END IF 428 CALL SLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 429 $ WORK( BX+ST1 ), N ) 430 ELSE 431* 432* A large problem. Solve it using divide and conquer. 433* 434 CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 435 $ E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ), 436 $ IWORK( K+ST1 ), WORK( DIFL+ST1 ), 437 $ WORK( DIFR+ST1 ), WORK( Z+ST1 ), 438 $ WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 439 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 440 $ WORK( GIVNUM+ST1 ), WORK( C+ST1 ), 441 $ WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ), 442 $ INFO ) 443 IF( INFO.NE.0 ) THEN 444 RETURN 445 END IF 446 BXST = BX + ST1 447 CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 448 $ LDB, WORK( BXST ), N, WORK( U+ST1 ), N, 449 $ WORK( VT+ST1 ), IWORK( K+ST1 ), 450 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ), 451 $ WORK( Z+ST1 ), WORK( POLES+ST1 ), 452 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 453 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ), 454 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ), 455 $ IWORK( IWK ), INFO ) 456 IF( INFO.NE.0 ) THEN 457 RETURN 458 END IF 459 END IF 460 ST = I + 1 461 END IF 462 60 CONTINUE 463* 464* Apply the singular values and treat the tiny ones as zero. 465* 466 TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) ) 467* 468 DO 70 I = 1, N 469* 470* Some of the elements in D can be negative because 1-by-1 471* subproblems were not solved explicitly. 472* 473 IF( ABS( D( I ) ).LE.TOL ) THEN 474 CALL SLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N ) 475 ELSE 476 RANK = RANK + 1 477 CALL SLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 478 $ WORK( BX+I-1 ), N, INFO ) 479 END IF 480 D( I ) = ABS( D( I ) ) 481 70 CONTINUE 482* 483* Now apply back the right singular vectors. 484* 485 ICMPQ2 = 1 486 DO 80 I = 1, NSUB 487 ST = IWORK( I ) 488 ST1 = ST - 1 489 NSIZE = IWORK( SIZEI+I-1 ) 490 BXST = BX + ST1 491 IF( NSIZE.EQ.1 ) THEN 492 CALL SCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 493 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 494 CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 495 $ WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO, 496 $ B( ST, 1 ), LDB ) 497 ELSE 498 CALL SLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 499 $ B( ST, 1 ), LDB, WORK( U+ST1 ), N, 500 $ WORK( VT+ST1 ), IWORK( K+ST1 ), 501 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ), 502 $ WORK( Z+ST1 ), WORK( POLES+ST1 ), 503 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 504 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ), 505 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ), 506 $ IWORK( IWK ), INFO ) 507 IF( INFO.NE.0 ) THEN 508 RETURN 509 END IF 510 END IF 511 80 CONTINUE 512* 513* Unscale and sort the singular values. 514* 515 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 516 CALL SLASRT( 'D', N, D, INFO ) 517 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 518* 519 RETURN 520* 521* End of SLALSD 522* 523 END 524