1*> \brief \b DLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLALSA + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsa.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsa.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsa.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 22* LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 23* GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, 24* IWORK, INFO ) 25* 26* .. Scalar Arguments .. 27* INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 28* $ SMLSIZ 29* .. 30* .. Array Arguments .. 31* INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 32* $ K( * ), PERM( LDGCOL, * ) 33* DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ), 34* $ DIFL( LDU, * ), DIFR( LDU, * ), 35* $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ), 36* $ U( LDU, * ), VT( LDU, * ), WORK( * ), 37* $ Z( LDU, * ) 38* .. 39* 40* 41*> \par Purpose: 42* ============= 43*> 44*> \verbatim 45*> 46*> DLALSA is an itermediate step in solving the least squares problem 47*> by computing the SVD of the coefficient matrix in compact form (The 48*> singular vectors are computed as products of simple orthorgonal 49*> matrices.). 50*> 51*> If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector 52*> matrix of an upper bidiagonal matrix to the right hand side; and if 53*> ICOMPQ = 1, DLALSA applies the right singular vector matrix to the 54*> right hand side. The singular vector matrices were generated in 55*> compact form by DLALSA. 56*> \endverbatim 57* 58* Arguments: 59* ========== 60* 61*> \param[in] ICOMPQ 62*> \verbatim 63*> ICOMPQ is INTEGER 64*> Specifies whether the left or the right singular vector 65*> matrix is involved. 66*> = 0: Left singular vector matrix 67*> = 1: Right singular vector matrix 68*> \endverbatim 69*> 70*> \param[in] SMLSIZ 71*> \verbatim 72*> SMLSIZ is INTEGER 73*> The maximum size of the subproblems at the bottom of the 74*> computation tree. 75*> \endverbatim 76*> 77*> \param[in] N 78*> \verbatim 79*> N is INTEGER 80*> The row and column dimensions of the upper bidiagonal matrix. 81*> \endverbatim 82*> 83*> \param[in] NRHS 84*> \verbatim 85*> NRHS is INTEGER 86*> The number of columns of B and BX. NRHS must be at least 1. 87*> \endverbatim 88*> 89*> \param[in,out] B 90*> \verbatim 91*> B is DOUBLE PRECISION array, dimension ( LDB, NRHS ) 92*> On input, B contains the right hand sides of the least 93*> squares problem in rows 1 through M. 94*> On output, B contains the solution X in rows 1 through N. 95*> \endverbatim 96*> 97*> \param[in] LDB 98*> \verbatim 99*> LDB is INTEGER 100*> The leading dimension of B in the calling subprogram. 101*> LDB must be at least max(1,MAX( M, N ) ). 102*> \endverbatim 103*> 104*> \param[out] BX 105*> \verbatim 106*> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS ) 107*> On exit, the result of applying the left or right singular 108*> vector matrix to B. 109*> \endverbatim 110*> 111*> \param[in] LDBX 112*> \verbatim 113*> LDBX is INTEGER 114*> The leading dimension of BX. 115*> \endverbatim 116*> 117*> \param[in] U 118*> \verbatim 119*> U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ). 120*> On entry, U contains the left singular vector matrices of all 121*> subproblems at the bottom level. 122*> \endverbatim 123*> 124*> \param[in] LDU 125*> \verbatim 126*> LDU is INTEGER, LDU = > N. 127*> The leading dimension of arrays U, VT, DIFL, DIFR, 128*> POLES, GIVNUM, and Z. 129*> \endverbatim 130*> 131*> \param[in] VT 132*> \verbatim 133*> VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ). 134*> On entry, VT**T contains the right singular vector matrices of 135*> all subproblems at the bottom level. 136*> \endverbatim 137*> 138*> \param[in] K 139*> \verbatim 140*> K is INTEGER array, dimension ( N ). 141*> \endverbatim 142*> 143*> \param[in] DIFL 144*> \verbatim 145*> DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ). 146*> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. 147*> \endverbatim 148*> 149*> \param[in] DIFR 150*> \verbatim 151*> DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 152*> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record 153*> distances between singular values on the I-th level and 154*> singular values on the (I -1)-th level, and DIFR(*, 2 * I) 155*> record the normalizing factors of the right singular vectors 156*> matrices of subproblems on I-th level. 157*> \endverbatim 158*> 159*> \param[in] Z 160*> \verbatim 161*> Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ). 162*> On entry, Z(1, I) contains the components of the deflation- 163*> adjusted updating row vector for subproblems on the I-th 164*> level. 165*> \endverbatim 166*> 167*> \param[in] POLES 168*> \verbatim 169*> POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 170*> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old 171*> singular values involved in the secular equations on the I-th 172*> level. 173*> \endverbatim 174*> 175*> \param[in] GIVPTR 176*> \verbatim 177*> GIVPTR is INTEGER array, dimension ( N ). 178*> On entry, GIVPTR( I ) records the number of Givens 179*> rotations performed on the I-th problem on the computation 180*> tree. 181*> \endverbatim 182*> 183*> \param[in] GIVCOL 184*> \verbatim 185*> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ). 186*> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the 187*> locations of Givens rotations performed on the I-th level on 188*> the computation tree. 189*> \endverbatim 190*> 191*> \param[in] LDGCOL 192*> \verbatim 193*> LDGCOL is INTEGER, LDGCOL = > N. 194*> The leading dimension of arrays GIVCOL and PERM. 195*> \endverbatim 196*> 197*> \param[in] PERM 198*> \verbatim 199*> PERM is INTEGER array, dimension ( LDGCOL, NLVL ). 200*> On entry, PERM(*, I) records permutations done on the I-th 201*> level of the computation tree. 202*> \endverbatim 203*> 204*> \param[in] GIVNUM 205*> \verbatim 206*> GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 207*> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- 208*> values of Givens rotations performed on the I-th level on the 209*> computation tree. 210*> \endverbatim 211*> 212*> \param[in] C 213*> \verbatim 214*> C is DOUBLE PRECISION array, dimension ( N ). 215*> On entry, if the I-th subproblem is not square, 216*> C( I ) contains the C-value of a Givens rotation related to 217*> the right null space of the I-th subproblem. 218*> \endverbatim 219*> 220*> \param[in] S 221*> \verbatim 222*> S is DOUBLE PRECISION array, dimension ( N ). 223*> On entry, if the I-th subproblem is not square, 224*> S( I ) contains the S-value of a Givens rotation related to 225*> the right null space of the I-th subproblem. 226*> \endverbatim 227*> 228*> \param[out] WORK 229*> \verbatim 230*> WORK is DOUBLE PRECISION array, dimension (N) 231*> \endverbatim 232*> 233*> \param[out] IWORK 234*> \verbatim 235*> IWORK is INTEGER array, dimension (3*N) 236*> \endverbatim 237*> 238*> \param[out] INFO 239*> \verbatim 240*> INFO is INTEGER 241*> = 0: successful exit. 242*> < 0: if INFO = -i, the i-th argument had an illegal value. 243*> \endverbatim 244* 245* Authors: 246* ======== 247* 248*> \author Univ. of Tennessee 249*> \author Univ. of California Berkeley 250*> \author Univ. of Colorado Denver 251*> \author NAG Ltd. 252* 253*> \ingroup doubleOTHERcomputational 254* 255*> \par Contributors: 256* ================== 257*> 258*> Ming Gu and Ren-Cang Li, Computer Science Division, University of 259*> California at Berkeley, USA \n 260*> Osni Marques, LBNL/NERSC, USA \n 261* 262* ===================================================================== 263 SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 264 $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 265 $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, 266 $ IWORK, INFO ) 267* 268* -- LAPACK computational routine -- 269* -- LAPACK is a software package provided by Univ. of Tennessee, -- 270* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 271* 272* .. Scalar Arguments .. 273 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 274 $ SMLSIZ 275* .. 276* .. Array Arguments .. 277 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 278 $ K( * ), PERM( LDGCOL, * ) 279 DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ), 280 $ DIFL( LDU, * ), DIFR( LDU, * ), 281 $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ), 282 $ U( LDU, * ), VT( LDU, * ), WORK( * ), 283 $ Z( LDU, * ) 284* .. 285* 286* ===================================================================== 287* 288* .. Parameters .. 289 DOUBLE PRECISION ZERO, ONE 290 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 291* .. 292* .. Local Scalars .. 293 INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2, 294 $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL, 295 $ NR, NRF, NRP1, SQRE 296* .. 297* .. External Subroutines .. 298 EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA 299* .. 300* .. Executable Statements .. 301* 302* Test the input parameters. 303* 304 INFO = 0 305* 306 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 307 INFO = -1 308 ELSE IF( SMLSIZ.LT.3 ) THEN 309 INFO = -2 310 ELSE IF( N.LT.SMLSIZ ) THEN 311 INFO = -3 312 ELSE IF( NRHS.LT.1 ) THEN 313 INFO = -4 314 ELSE IF( LDB.LT.N ) THEN 315 INFO = -6 316 ELSE IF( LDBX.LT.N ) THEN 317 INFO = -8 318 ELSE IF( LDU.LT.N ) THEN 319 INFO = -10 320 ELSE IF( LDGCOL.LT.N ) THEN 321 INFO = -19 322 END IF 323 IF( INFO.NE.0 ) THEN 324 CALL XERBLA( 'DLALSA', -INFO ) 325 RETURN 326 END IF 327* 328* Book-keeping and setting up the computation tree. 329* 330 INODE = 1 331 NDIML = INODE + N 332 NDIMR = NDIML + N 333* 334 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 335 $ IWORK( NDIMR ), SMLSIZ ) 336* 337* The following code applies back the left singular vector factors. 338* For applying back the right singular vector factors, go to 50. 339* 340 IF( ICOMPQ.EQ.1 ) THEN 341 GO TO 50 342 END IF 343* 344* The nodes on the bottom level of the tree were solved 345* by DLASDQ. The corresponding left and right singular vector 346* matrices are in explicit form. First apply back the left 347* singular vector matrices. 348* 349 NDB1 = ( ND+1 ) / 2 350 DO 10 I = NDB1, ND 351* 352* IC : center row of each node 353* NL : number of rows of left subproblem 354* NR : number of rows of right subproblem 355* NLF: starting row of the left subproblem 356* NRF: starting row of the right subproblem 357* 358 I1 = I - 1 359 IC = IWORK( INODE+I1 ) 360 NL = IWORK( NDIML+I1 ) 361 NR = IWORK( NDIMR+I1 ) 362 NLF = IC - NL 363 NRF = IC + 1 364 CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 365 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 366 CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 367 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 368 10 CONTINUE 369* 370* Next copy the rows of B that correspond to unchanged rows 371* in the bidiagonal matrix to BX. 372* 373 DO 20 I = 1, ND 374 IC = IWORK( INODE+I-1 ) 375 CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX ) 376 20 CONTINUE 377* 378* Finally go through the left singular vector matrices of all 379* the other subproblems bottom-up on the tree. 380* 381 J = 2**NLVL 382 SQRE = 0 383* 384 DO 40 LVL = NLVL, 1, -1 385 LVL2 = 2*LVL - 1 386* 387* find the first node LF and last node LL on 388* the current level LVL 389* 390 IF( LVL.EQ.1 ) THEN 391 LF = 1 392 LL = 1 393 ELSE 394 LF = 2**( LVL-1 ) 395 LL = 2*LF - 1 396 END IF 397 DO 30 I = LF, LL 398 IM1 = I - 1 399 IC = IWORK( INODE+IM1 ) 400 NL = IWORK( NDIML+IM1 ) 401 NR = IWORK( NDIMR+IM1 ) 402 NLF = IC - NL 403 NRF = IC + 1 404 J = J - 1 405 CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX, 406 $ B( NLF, 1 ), LDB, PERM( NLF, LVL ), 407 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 408 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 409 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 410 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 411 $ INFO ) 412 30 CONTINUE 413 40 CONTINUE 414 GO TO 90 415* 416* ICOMPQ = 1: applying back the right singular vector factors. 417* 418 50 CONTINUE 419* 420* First now go through the right singular vector matrices of all 421* the tree nodes top-down. 422* 423 J = 0 424 DO 70 LVL = 1, NLVL 425 LVL2 = 2*LVL - 1 426* 427* Find the first node LF and last node LL on 428* the current level LVL. 429* 430 IF( LVL.EQ.1 ) THEN 431 LF = 1 432 LL = 1 433 ELSE 434 LF = 2**( LVL-1 ) 435 LL = 2*LF - 1 436 END IF 437 DO 60 I = LL, LF, -1 438 IM1 = I - 1 439 IC = IWORK( INODE+IM1 ) 440 NL = IWORK( NDIML+IM1 ) 441 NR = IWORK( NDIMR+IM1 ) 442 NLF = IC - NL 443 NRF = IC + 1 444 IF( I.EQ.LL ) THEN 445 SQRE = 0 446 ELSE 447 SQRE = 1 448 END IF 449 J = J + 1 450 CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB, 451 $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ), 452 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 453 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 454 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 455 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 456 $ INFO ) 457 60 CONTINUE 458 70 CONTINUE 459* 460* The nodes on the bottom level of the tree were solved 461* by DLASDQ. The corresponding right singular vector 462* matrices are in explicit form. Apply them back. 463* 464 NDB1 = ( ND+1 ) / 2 465 DO 80 I = NDB1, ND 466 I1 = I - 1 467 IC = IWORK( INODE+I1 ) 468 NL = IWORK( NDIML+I1 ) 469 NR = IWORK( NDIMR+I1 ) 470 NLP1 = NL + 1 471 IF( I.EQ.ND ) THEN 472 NRP1 = NR 473 ELSE 474 NRP1 = NR + 1 475 END IF 476 NLF = IC - NL 477 NRF = IC + 1 478 CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 479 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 480 CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 481 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 482 80 CONTINUE 483* 484 90 CONTINUE 485* 486 RETURN 487* 488* End of DLALSA 489* 490 END 491