SRC/slals0.f

   1 *> \brief \b SLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. 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 SLALS0 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slals0.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slals0.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slals0.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
  22 *                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
  23 *                          POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
  27 *      $                   LDGNUM, NL, NR, NRHS, SQRE
  28 *       REAL               C, S
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
  32 *       REAL               B( LDB, * ), BX( LDBX, * ), DIFL( * ),
  33 *      $                   DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
  34 *      $                   POLES( LDGNUM, * ), WORK( * ), Z( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> SLALS0 applies back the multiplying factors of either the left or the
  44 *> right singular vector matrix of a diagonal matrix appended by a row
  45 *> to the right hand side matrix B in solving the least squares problem
  46 *> using the divide-and-conquer SVD approach.
  47 *>
  48 *> For the left singular vector matrix, three types of orthogonal
  49 *> matrices are involved:
  50 *>
  51 *> (1L) Givens rotations: the number of such rotations is GIVPTR; the
  52 *>      pairs of columns/rows they were applied to are stored in GIVCOL;
  53 *>      and the C- and S-values of these rotations are stored in GIVNUM.
  54 *>
  55 *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
  56 *>      row, and for J=2:N, PERM(J)-th row of B is to be moved to the
  57 *>      J-th row.
  58 *>
  59 *> (3L) The left singular vector matrix of the remaining matrix.
  60 *>
  61 *> For the right singular vector matrix, four types of orthogonal
  62 *> matrices are involved:
  63 *>
  64 *> (1R) The right singular vector matrix of the remaining matrix.
  65 *>
  66 *> (2R) If SQRE = 1, one extra Givens rotation to generate the right
  67 *>      null space.
  68 *>
  69 *> (3R) The inverse transformation of (2L).
  70 *>
  71 *> (4R) The inverse transformation of (1L).
  72 *> \endverbatim
  73 *
  74 *  Arguments:
  75 *  ==========
  76 *
  77 *> \param[in] ICOMPQ
  78 *> \verbatim
  79 *>          ICOMPQ is INTEGER
  80 *>         Specifies whether singular vectors are to be computed in
  81 *>         factored form:
  82 *>         = 0: Left singular vector matrix.
  83 *>         = 1: Right singular vector matrix.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] NL
  87 *> \verbatim
  88 *>          NL is INTEGER
  89 *>         The row dimension of the upper block. NL >= 1.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] NR
  93 *> \verbatim
  94 *>          NR is INTEGER
  95 *>         The row dimension of the lower block. NR >= 1.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] SQRE
  99 *> \verbatim
 100 *>          SQRE is INTEGER
 101 *>         = 0: the lower block is an NR-by-NR square matrix.
 102 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 103 *>
 104 *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
 105 *>         and column dimension M = N + SQRE.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] NRHS
 109 *> \verbatim
 110 *>          NRHS is INTEGER
 111 *>         The number of columns of B and BX. NRHS must be at least 1.
 112 *> \endverbatim
 113 *>
 114 *> \param[in,out] B
 115 *> \verbatim
 116 *>          B is REAL array, dimension ( LDB, NRHS )
 117 *>         On input, B contains the right hand sides of the least
 118 *>         squares problem in rows 1 through M. On output, B contains
 119 *>         the solution X in rows 1 through N.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] LDB
 123 *> \verbatim
 124 *>          LDB is INTEGER
 125 *>         The leading dimension of B. LDB must be at least
 126 *>         max(1,MAX( M, N ) ).
 127 *> \endverbatim
 128 *>
 129 *> \param[out] BX
 130 *> \verbatim
 131 *>          BX is REAL array, dimension ( LDBX, NRHS )
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LDBX
 135 *> \verbatim
 136 *>          LDBX is INTEGER
 137 *>         The leading dimension of BX.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] PERM
 141 *> \verbatim
 142 *>          PERM is INTEGER array, dimension ( N )
 143 *>         The permutations (from deflation and sorting) applied
 144 *>         to the two blocks.
 145 *> \endverbatim
 146 *>
 147 *> \param[in] GIVPTR
 148 *> \verbatim
 149 *>          GIVPTR is INTEGER
 150 *>         The number of Givens rotations which took place in this
 151 *>         subproblem.
 152 *> \endverbatim
 153 *>
 154 *> \param[in] GIVCOL
 155 *> \verbatim
 156 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
 157 *>         Each pair of numbers indicates a pair of rows/columns
 158 *>         involved in a Givens rotation.
 159 *> \endverbatim
 160 *>
 161 *> \param[in] LDGCOL
 162 *> \verbatim
 163 *>          LDGCOL is INTEGER
 164 *>         The leading dimension of GIVCOL, must be at least N.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] GIVNUM
 168 *> \verbatim
 169 *>          GIVNUM is REAL array, dimension ( LDGNUM, 2 )
 170 *>         Each number indicates the C or S value used in the
 171 *>         corresponding Givens rotation.
 172 *> \endverbatim
 173 *>
 174 *> \param[in] LDGNUM
 175 *> \verbatim
 176 *>          LDGNUM is INTEGER
 177 *>         The leading dimension of arrays DIFR, POLES and
 178 *>         GIVNUM, must be at least K.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] POLES
 182 *> \verbatim
 183 *>          POLES is REAL array, dimension ( LDGNUM, 2 )
 184 *>         On entry, POLES(1:K, 1) contains the new singular
 185 *>         values obtained from solving the secular equation, and
 186 *>         POLES(1:K, 2) is an array containing the poles in the secular
 187 *>         equation.
 188 *> \endverbatim
 189 *>
 190 *> \param[in] DIFL
 191 *> \verbatim
 192 *>          DIFL is REAL array, dimension ( K ).
 193 *>         On entry, DIFL(I) is the distance between I-th updated
 194 *>         (undeflated) singular value and the I-th (undeflated) old
 195 *>         singular value.
 196 *> \endverbatim
 197 *>
 198 *> \param[in] DIFR
 199 *> \verbatim
 200 *>          DIFR is REAL array, dimension ( LDGNUM, 2 ).
 201 *>         On entry, DIFR(I, 1) contains the distances between I-th
 202 *>         updated (undeflated) singular value and the I+1-th
 203 *>         (undeflated) old singular value. And DIFR(I, 2) is the
 204 *>         normalizing factor for the I-th right singular vector.
 205 *> \endverbatim
 206 *>
 207 *> \param[in] Z
 208 *> \verbatim
 209 *>          Z is REAL array, dimension ( K )
 210 *>         Contain the components of the deflation-adjusted updating row
 211 *>         vector.
 212 *> \endverbatim
 213 *>
 214 *> \param[in] K
 215 *> \verbatim
 216 *>          K is INTEGER
 217 *>         Contains the dimension of the non-deflated matrix,
 218 *>         This is the order of the related secular equation. 1 <= K <=N.
 219 *> \endverbatim
 220 *>
 221 *> \param[in] C
 222 *> \verbatim
 223 *>          C is REAL
 224 *>         C contains garbage if SQRE =0 and the C-value of a Givens
 225 *>         rotation related to the right null space if SQRE = 1.
 226 *> \endverbatim
 227 *>
 228 *> \param[in] S
 229 *> \verbatim
 230 *>          S is REAL
 231 *>         S contains garbage if SQRE =0 and the S-value of a Givens
 232 *>         rotation related to the right null space if SQRE = 1.
 233 *> \endverbatim
 234 *>
 235 *> \param[out] WORK
 236 *> \verbatim
 237 *>          WORK is REAL array, dimension ( K )
 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 November 2015
 256 *
 257 *> \ingroup realOTHERcomputational
 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 SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
 268      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
 269      $                   POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
 270 *
 271 *  -- LAPACK computational routine (version 3.6.0) --
 272 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 273 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 274 *     November 2015
 275 *
 276 *     .. Scalar Arguments ..
 277       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
 278      $                   LDGNUM, NL, NR, NRHS, SQRE
 279       REAL               C, S
 280 *     ..
 281 *     .. Array Arguments ..
 282       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
 283       REAL               B( LDB, * ), BX( LDBX, * ), DIFL( * ),
 284      $                   DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
 285      $                   POLES( LDGNUM, * ), WORK( * ), Z( * )
 286 *     ..
 287 *
 288 *  =====================================================================
 289 *
 290 *     .. Parameters ..
 291       REAL               ONE, ZERO, NEGONE
 292       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
 293 *     ..
 294 *     .. Local Scalars ..
 295       INTEGER            I, J, M, N, NLP1
 296       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
 297 *     ..
 298 *     .. External Subroutines ..
 299       EXTERNAL           SCOPY, SGEMV, SLACPY, SLASCL, SROT, SSCAL,
 300      $                   XERBLA
 301 *     ..
 302 *     .. External Functions ..
 303       REAL               SLAMC3, SNRM2
 304       EXTERNAL           SLAMC3, SNRM2
 305 *     ..
 306 *     .. Intrinsic Functions ..
 307       INTRINSIC          MAX
 308 *     ..
 309 *     .. Executable Statements ..
 310 *
 311 *     Test the input parameters.
 312 *
 313       INFO = 0
 314       N = NL + NR + 1
 315 *
 316       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 317          INFO = -1
 318       ELSE IF( NL.LT.1 ) THEN
 319          INFO = -2
 320       ELSE IF( NR.LT.1 ) THEN
 321          INFO = -3
 322       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
 323          INFO = -4
 324       ELSE IF( NRHS.LT.1 ) THEN
 325          INFO = -5
 326       ELSE IF( LDB.LT.N ) THEN
 327          INFO = -7
 328       ELSE IF( LDBX.LT.N ) THEN
 329          INFO = -9
 330       ELSE IF( GIVPTR.LT.0 ) THEN
 331          INFO = -11
 332       ELSE IF( LDGCOL.LT.N ) THEN
 333          INFO = -13
 334       ELSE IF( LDGNUM.LT.N ) THEN
 335          INFO = -15
 336       ELSE IF( K.LT.1 ) THEN
 337          INFO = -20
 338       END IF
 339       IF( INFO.NE.0 ) THEN
 340          CALL XERBLA( 'SLALS0', -INFO )
 341          RETURN
 342       END IF
 343 *
 344       M = N + SQRE
 345       NLP1 = NL + 1
 346 *
 347       IF( ICOMPQ.EQ.0 ) THEN
 348 *
 349 *        Apply back orthogonal transformations from the left.
 350 *
 351 *        Step (1L): apply back the Givens rotations performed.
 352 *
 353          DO 10 I = 1, GIVPTR
 354             CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
 355      $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
 356      $                 GIVNUM( I, 1 ) )
 357    10    CONTINUE
 358 *
 359 *        Step (2L): permute rows of B.
 360 *
 361          CALL SCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
 362          DO 20 I = 2, N
 363             CALL SCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
 364    20    CONTINUE
 365 *
 366 *        Step (3L): apply the inverse of the left singular vector
 367 *        matrix to BX.
 368 *
 369          IF( K.EQ.1 ) THEN
 370             CALL SCOPY( NRHS, BX, LDBX, B, LDB )
 371             IF( Z( 1 ).LT.ZERO ) THEN
 372                CALL SSCAL( NRHS, NEGONE, B, LDB )
 373             END IF
 374          ELSE
 375             DO 50 J = 1, K
 376                DIFLJ = DIFL( J )
 377                DJ = POLES( J, 1 )
 378                DSIGJ = -POLES( J, 2 )
 379                IF( J.LT.K ) THEN
 380                   DIFRJ = -DIFR( J, 1 )
 381                   DSIGJP = -POLES( J+1, 2 )
 382                END IF
 383                IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
 384      $              THEN
 385                   WORK( J ) = ZERO
 386                ELSE
 387                   WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
 388      $                        ( POLES( J, 2 )+DJ )
 389                END IF
 390                DO 30 I = 1, J - 1
 391                   IF( ( Z( I ).EQ.ZERO ) .OR.
 392      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
 393                      WORK( I ) = ZERO
 394                   ELSE
 395                      WORK( I ) = POLES( I, 2 )*Z( I ) /
 396      $                           ( SLAMC3( POLES( I, 2 ), DSIGJ )-
 397      $                           DIFLJ ) / ( POLES( I, 2 )+DJ )
 398                   END IF
 399    30          CONTINUE
 400                DO 40 I = J + 1, K
 401                   IF( ( Z( I ).EQ.ZERO ) .OR.
 402      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
 403                      WORK( I ) = ZERO
 404                   ELSE
 405                      WORK( I ) = POLES( I, 2 )*Z( I ) /
 406      $                           ( SLAMC3( POLES( I, 2 ), DSIGJP )+
 407      $                           DIFRJ ) / ( POLES( I, 2 )+DJ )
 408                   END IF
 409    40          CONTINUE
 410                WORK( 1 ) = NEGONE
 411                TEMP = SNRM2( K, WORK, 1 )
 412                CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
 413      $                     B( J, 1 ), LDB )
 414                CALL SLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
 415      $                      LDB, INFO )
 416    50       CONTINUE
 417          END IF
 418 *
 419 *        Move the deflated rows of BX to B also.
 420 *
 421          IF( K.LT.MAX( M, N ) )
 422      $      CALL SLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
 423      $                   B( K+1, 1 ), LDB )
 424       ELSE
 425 *
 426 *        Apply back the right orthogonal transformations.
 427 *
 428 *        Step (1R): apply back the new right singular vector matrix
 429 *        to B.
 430 *
 431          IF( K.EQ.1 ) THEN
 432             CALL SCOPY( NRHS, B, LDB, BX, LDBX )
 433          ELSE
 434             DO 80 J = 1, K
 435                DSIGJ = POLES( J, 2 )
 436                IF( Z( J ).EQ.ZERO ) THEN
 437                   WORK( J ) = ZERO
 438                ELSE
 439                   WORK( J ) = -Z( J ) / DIFL( J ) /
 440      $                        ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
 441                END IF
 442                DO 60 I = 1, J - 1
 443                   IF( Z( J ).EQ.ZERO ) THEN
 444                      WORK( I ) = ZERO
 445                   ELSE
 446                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1,
 447      $                           2 ) )-DIFR( I, 1 ) ) /
 448      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
 449                   END IF
 450    60          CONTINUE
 451                DO 70 I = J + 1, K
 452                   IF( Z( J ).EQ.ZERO ) THEN
 453                      WORK( I ) = ZERO
 454                   ELSE
 455                      WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I,
 456      $                           2 ) )-DIFL( I ) ) /
 457      $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
 458                   END IF
 459    70          CONTINUE
 460                CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
 461      $                     BX( J, 1 ), LDBX )
 462    80       CONTINUE
 463          END IF
 464 *
 465 *        Step (2R): if SQRE = 1, apply back the rotation that is
 466 *        related to the right null space of the subproblem.
 467 *
 468          IF( SQRE.EQ.1 ) THEN
 469             CALL SCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
 470             CALL SROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
 471          END IF
 472          IF( K.LT.MAX( M, N ) )
 473      $      CALL SLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
 474      $                   LDBX )
 475 *
 476 *        Step (3R): permute rows of B.
 477 *
 478          CALL SCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
 479          IF( SQRE.EQ.1 ) THEN
 480             CALL SCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
 481          END IF
 482          DO 90 I = 2, N
 483             CALL SCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
 484    90    CONTINUE
 485 *
 486 *        Step (4R): apply back the Givens rotations performed.
 487 *
 488          DO 100 I = GIVPTR, 1, -1
 489             CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
 490      $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
 491      $                 -GIVNUM( I, 1 ) )
 492   100    CONTINUE
 493       END IF
 494 *
 495       RETURN
 496 *
 497 *     End of SLALS0
 498 *
 499       END