SRC/zlals0.f

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