SRC/dlalsa.f

   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.
 231 *>         The dimension must be at least N.
 232 *> \endverbatim
 233 *>
 234 *> \param[out] IWORK
 235 *> \verbatim
 236 *>          IWORK is INTEGER array.
 237 *>         The dimension must be at least 3 * N
 238 *> \endverbatim
 239 *>
 240 *> \param[out] INFO
 241 *> \verbatim
 242 *>          INFO is INTEGER
 243 *>          = 0:  successful exit.
 244 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 245 *> \endverbatim
 246 *
 247 *  Authors:
 248 *  ========
 249 *
 250 *> \author Univ. of Tennessee
 251 *> \author Univ. of California Berkeley
 252 *> \author Univ. of Colorado Denver
 253 *> \author NAG Ltd.
 254 *
 255 *> \date September 2012
 256 *
 257 *> \ingroup doubleOTHERcomputational
 258 *
 259 *> \par Contributors:
 260 *  ==================
 261 *>
 262 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
 263 *>       California at Berkeley, USA \n
 264 *>     Osni Marques, LBNL/NERSC, USA \n
 265 *
 266 *  =====================================================================
 267       SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
 268      $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
 269      $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
 270      $                   IWORK, INFO )
 271 *
 272 *  -- LAPACK computational routine (version 3.4.2) --
 273 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 274 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 275 *     September 2012
 276 *
 277 *     .. Scalar Arguments ..
 278       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
 279      $                   SMLSIZ
 280 *     ..
 281 *     .. Array Arguments ..
 282       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
 283      $                   K( * ), PERM( LDGCOL, * )
 284       DOUBLE PRECISION   B( LDB, * ), BX( LDBX, * ), C( * ),
 285      $                   DIFL( LDU, * ), DIFR( LDU, * ),
 286      $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
 287      $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
 288      $                   Z( LDU, * )
 289 *     ..
 290 *
 291 *  =====================================================================
 292 *
 293 *     .. Parameters ..
 294       DOUBLE PRECISION   ZERO, ONE
 295       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 296 *     ..
 297 *     .. Local Scalars ..
 298       INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
 299      $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
 300      $                   NR, NRF, NRP1, SQRE
 301 *     ..
 302 *     .. External Subroutines ..
 303       EXTERNAL           DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
 304 *     ..
 305 *     .. Executable Statements ..
 306 *
 307 *     Test the input parameters.
 308 *
 309       INFO = 0
 310 *
 311       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 312          INFO = -1
 313       ELSE IF( SMLSIZ.LT.3 ) THEN
 314          INFO = -2
 315       ELSE IF( N.LT.SMLSIZ ) THEN
 316          INFO = -3
 317       ELSE IF( NRHS.LT.1 ) THEN
 318          INFO = -4
 319       ELSE IF( LDB.LT.N ) THEN
 320          INFO = -6
 321       ELSE IF( LDBX.LT.N ) THEN
 322          INFO = -8
 323       ELSE IF( LDU.LT.N ) THEN
 324          INFO = -10
 325       ELSE IF( LDGCOL.LT.N ) THEN
 326          INFO = -19
 327       END IF
 328       IF( INFO.NE.0 ) THEN
 329          CALL XERBLA( 'DLALSA', -INFO )
 330          RETURN
 331       END IF
 332 *
 333 *     Book-keeping and  setting up the computation tree.
 334 *
 335       INODE = 1
 336       NDIML = INODE + N
 337       NDIMR = NDIML + N
 338 *
 339       CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
 340      $             IWORK( NDIMR ), SMLSIZ )
 341 *
 342 *     The following code applies back the left singular vector factors.
 343 *     For applying back the right singular vector factors, go to 50.
 344 *
 345       IF( ICOMPQ.EQ.1 ) THEN
 346          GO TO 50
 347       END IF
 348 *
 349 *     The nodes on the bottom level of the tree were solved
 350 *     by DLASDQ. The corresponding left and right singular vector
 351 *     matrices are in explicit form. First apply back the left
 352 *     singular vector matrices.
 353 *
 354       NDB1 = ( ND+1 ) / 2
 355       DO 10 I = NDB1, ND
 356 *
 357 *        IC : center row of each node
 358 *        NL : number of rows of left  subproblem
 359 *        NR : number of rows of right subproblem
 360 *        NLF: starting row of the left   subproblem
 361 *        NRF: starting row of the right  subproblem
 362 *
 363          I1 = I - 1
 364          IC = IWORK( INODE+I1 )
 365          NL = IWORK( NDIML+I1 )
 366          NR = IWORK( NDIMR+I1 )
 367          NLF = IC - NL
 368          NRF = IC + 1
 369          CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
 370      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
 371          CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
 372      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
 373    10 CONTINUE
 374 *
 375 *     Next copy the rows of B that correspond to unchanged rows
 376 *     in the bidiagonal matrix to BX.
 377 *
 378       DO 20 I = 1, ND
 379          IC = IWORK( INODE+I-1 )
 380          CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
 381    20 CONTINUE
 382 *
 383 *     Finally go through the left singular vector matrices of all
 384 *     the other subproblems bottom-up on the tree.
 385 *
 386       J = 2**NLVL
 387       SQRE = 0
 388 *
 389       DO 40 LVL = NLVL, 1, -1
 390          LVL2 = 2*LVL - 1
 391 *
 392 *        find the first node LF and last node LL on
 393 *        the current level LVL
 394 *
 395          IF( LVL.EQ.1 ) THEN
 396             LF = 1
 397             LL = 1
 398          ELSE
 399             LF = 2**( LVL-1 )
 400             LL = 2*LF - 1
 401          END IF
 402          DO 30 I = LF, LL
 403             IM1 = I - 1
 404             IC = IWORK( INODE+IM1 )
 405             NL = IWORK( NDIML+IM1 )
 406             NR = IWORK( NDIMR+IM1 )
 407             NLF = IC - NL
 408             NRF = IC + 1
 409             J = J - 1
 410             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
 411      $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
 412      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
 413      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
 414      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
 415      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
 416      $                   INFO )
 417    30    CONTINUE
 418    40 CONTINUE
 419       GO TO 90
 420 *
 421 *     ICOMPQ = 1: applying back the right singular vector factors.
 422 *
 423    50 CONTINUE
 424 *
 425 *     First now go through the right singular vector matrices of all
 426 *     the tree nodes top-down.
 427 *
 428       J = 0
 429       DO 70 LVL = 1, NLVL
 430          LVL2 = 2*LVL - 1
 431 *
 432 *        Find the first node LF and last node LL on
 433 *        the current level LVL.
 434 *
 435          IF( LVL.EQ.1 ) THEN
 436             LF = 1
 437             LL = 1
 438          ELSE
 439             LF = 2**( LVL-1 )
 440             LL = 2*LF - 1
 441          END IF
 442          DO 60 I = LL, LF, -1
 443             IM1 = I - 1
 444             IC = IWORK( INODE+IM1 )
 445             NL = IWORK( NDIML+IM1 )
 446             NR = IWORK( NDIMR+IM1 )
 447             NLF = IC - NL
 448             NRF = IC + 1
 449             IF( I.EQ.LL ) THEN
 450                SQRE = 0
 451             ELSE
 452                SQRE = 1
 453             END IF
 454             J = J + 1
 455             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
 456      $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
 457      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
 458      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
 459      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
 460      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
 461      $                   INFO )
 462    60    CONTINUE
 463    70 CONTINUE
 464 *
 465 *     The nodes on the bottom level of the tree were solved
 466 *     by DLASDQ. The corresponding right singular vector
 467 *     matrices are in explicit form. Apply them back.
 468 *
 469       NDB1 = ( ND+1 ) / 2
 470       DO 80 I = NDB1, ND
 471          I1 = I - 1
 472          IC = IWORK( INODE+I1 )
 473          NL = IWORK( NDIML+I1 )
 474          NR = IWORK( NDIMR+I1 )
 475          NLP1 = NL + 1
 476          IF( I.EQ.ND ) THEN
 477             NRP1 = NR
 478          ELSE
 479             NRP1 = NR + 1
 480          END IF
 481          NLF = IC - NL
 482          NRF = IC + 1
 483          CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
 484      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
 485          CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
 486      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
 487    80 CONTINUE
 488 *
 489    90 CONTINUE
 490 *
 491       RETURN
 492 *
 493 *     End of DLALSA
 494 *
 495       END