SRC/clalsa.f

   1 *> \brief \b CLALSA 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 CLALSA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsa.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsa.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsa.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLALSA( 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, RWORK,
  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 *       REAL               C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
  34 *      $                   GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
  35 *      $                   S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
  36 *       COMPLEX            B( LDB, * ), BX( LDBX, * )
  37 *       ..
  38 *
  39 *
  40 *> \par Purpose:
  41 *  =============
  42 *>
  43 *> \verbatim
  44 *>
  45 *> CLALSA is an itermediate step in solving the least squares problem
  46 *> by computing the SVD of the coefficient matrix in compact form (The
  47 *> singular vectors are computed as products of simple orthorgonal
  48 *> matrices.).
  49 *>
  50 *> If ICOMPQ = 0, CLALSA applies the inverse of the left singular vector
  51 *> matrix of an upper bidiagonal matrix to the right hand side; and if
  52 *> ICOMPQ = 1, CLALSA applies the right singular vector matrix to the
  53 *> right hand side. The singular vector matrices were generated in
  54 *> compact form by CLALSA.
  55 *> \endverbatim
  56 *
  57 *  Arguments:
  58 *  ==========
  59 *
  60 *> \param[in] ICOMPQ
  61 *> \verbatim
  62 *>          ICOMPQ is INTEGER
  63 *>         Specifies whether the left or the right singular vector
  64 *>         matrix is involved.
  65 *>         = 0: Left singular vector matrix
  66 *>         = 1: Right singular vector matrix
  67 *> \endverbatim
  68 *>
  69 *> \param[in] SMLSIZ
  70 *> \verbatim
  71 *>          SMLSIZ is INTEGER
  72 *>         The maximum size of the subproblems at the bottom of the
  73 *>         computation tree.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] N
  77 *> \verbatim
  78 *>          N is INTEGER
  79 *>         The row and column dimensions of the upper bidiagonal matrix.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] NRHS
  83 *> \verbatim
  84 *>          NRHS is INTEGER
  85 *>         The number of columns of B and BX. NRHS must be at least 1.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] B
  89 *> \verbatim
  90 *>          B is COMPLEX array, dimension ( LDB, NRHS )
  91 *>         On input, B contains the right hand sides of the least
  92 *>         squares problem in rows 1 through M.
  93 *>         On output, B contains the solution X in rows 1 through N.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] LDB
  97 *> \verbatim
  98 *>          LDB is INTEGER
  99 *>         The leading dimension of B in the calling subprogram.
 100 *>         LDB must be at least max(1,MAX( M, N ) ).
 101 *> \endverbatim
 102 *>
 103 *> \param[out] BX
 104 *> \verbatim
 105 *>          BX is COMPLEX array, dimension ( LDBX, NRHS )
 106 *>         On exit, the result of applying the left or right singular
 107 *>         vector matrix to B.
 108 *> \endverbatim
 109 *>
 110 *> \param[in] LDBX
 111 *> \verbatim
 112 *>          LDBX is INTEGER
 113 *>         The leading dimension of BX.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] U
 117 *> \verbatim
 118 *>          U is REAL array, dimension ( LDU, SMLSIZ ).
 119 *>         On entry, U contains the left singular vector matrices of all
 120 *>         subproblems at the bottom level.
 121 *> \endverbatim
 122 *>
 123 *> \param[in] LDU
 124 *> \verbatim
 125 *>          LDU is INTEGER, LDU = > N.
 126 *>         The leading dimension of arrays U, VT, DIFL, DIFR,
 127 *>         POLES, GIVNUM, and Z.
 128 *> \endverbatim
 129 *>
 130 *> \param[in] VT
 131 *> \verbatim
 132 *>          VT is REAL array, dimension ( LDU, SMLSIZ+1 ).
 133 *>         On entry, VT**H contains the right singular vector matrices of
 134 *>         all subproblems at the bottom level.
 135 *> \endverbatim
 136 *>
 137 *> \param[in] K
 138 *> \verbatim
 139 *>          K is INTEGER array, dimension ( N ).
 140 *> \endverbatim
 141 *>
 142 *> \param[in] DIFL
 143 *> \verbatim
 144 *>          DIFL is REAL array, dimension ( LDU, NLVL ).
 145 *>         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
 146 *> \endverbatim
 147 *>
 148 *> \param[in] DIFR
 149 *> \verbatim
 150 *>          DIFR is REAL array, dimension ( LDU, 2 * NLVL ).
 151 *>         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
 152 *>         distances between singular values on the I-th level and
 153 *>         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
 154 *>         record the normalizing factors of the right singular vectors
 155 *>         matrices of subproblems on I-th level.
 156 *> \endverbatim
 157 *>
 158 *> \param[in] Z
 159 *> \verbatim
 160 *>          Z is REAL array, dimension ( LDU, NLVL ).
 161 *>         On entry, Z(1, I) contains the components of the deflation-
 162 *>         adjusted updating row vector for subproblems on the I-th
 163 *>         level.
 164 *> \endverbatim
 165 *>
 166 *> \param[in] POLES
 167 *> \verbatim
 168 *>          POLES is REAL array, dimension ( LDU, 2 * NLVL ).
 169 *>         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
 170 *>         singular values involved in the secular equations on the I-th
 171 *>         level.
 172 *> \endverbatim
 173 *>
 174 *> \param[in] GIVPTR
 175 *> \verbatim
 176 *>          GIVPTR is INTEGER array, dimension ( N ).
 177 *>         On entry, GIVPTR( I ) records the number of Givens
 178 *>         rotations performed on the I-th problem on the computation
 179 *>         tree.
 180 *> \endverbatim
 181 *>
 182 *> \param[in] GIVCOL
 183 *> \verbatim
 184 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
 185 *>         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
 186 *>         locations of Givens rotations performed on the I-th level on
 187 *>         the computation tree.
 188 *> \endverbatim
 189 *>
 190 *> \param[in] LDGCOL
 191 *> \verbatim
 192 *>          LDGCOL is INTEGER, LDGCOL = > N.
 193 *>         The leading dimension of arrays GIVCOL and PERM.
 194 *> \endverbatim
 195 *>
 196 *> \param[in] PERM
 197 *> \verbatim
 198 *>          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
 199 *>         On entry, PERM(*, I) records permutations done on the I-th
 200 *>         level of the computation tree.
 201 *> \endverbatim
 202 *>
 203 *> \param[in] GIVNUM
 204 *> \verbatim
 205 *>          GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ).
 206 *>         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
 207 *>         values of Givens rotations performed on the I-th level on the
 208 *>         computation tree.
 209 *> \endverbatim
 210 *>
 211 *> \param[in] C
 212 *> \verbatim
 213 *>          C is REAL array, dimension ( N ).
 214 *>         On entry, if the I-th subproblem is not square,
 215 *>         C( I ) contains the C-value of a Givens rotation related to
 216 *>         the right null space of the I-th subproblem.
 217 *> \endverbatim
 218 *>
 219 *> \param[in] S
 220 *> \verbatim
 221 *>          S is REAL array, dimension ( N ).
 222 *>         On entry, if the I-th subproblem is not square,
 223 *>         S( I ) contains the S-value of a Givens rotation related to
 224 *>         the right null space of the I-th subproblem.
 225 *> \endverbatim
 226 *>
 227 *> \param[out] RWORK
 228 *> \verbatim
 229 *>          RWORK is REAL array, dimension at least
 230 *>         MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ).
 231 *> \endverbatim
 232 *>
 233 *> \param[out] IWORK
 234 *> \verbatim
 235 *>          IWORK is INTEGER array.
 236 *>         The dimension must be at least 3 * N
 237 *> \endverbatim
 238 *>
 239 *> \param[out] INFO
 240 *> \verbatim
 241 *>          INFO is INTEGER
 242 *>          = 0:  successful exit.
 243 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 244 *> \endverbatim
 245 *
 246 *  Authors:
 247 *  ========
 248 *
 249 *> \author Univ. of Tennessee
 250 *> \author Univ. of California Berkeley
 251 *> \author Univ. of Colorado Denver
 252 *> \author NAG Ltd.
 253 *
 254 *> \date September 2012
 255 *
 256 *> \ingroup complexOTHERcomputational
 257 *
 258 *> \par Contributors:
 259 *  ==================
 260 *>
 261 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
 262 *>       California at Berkeley, USA \n
 263 *>     Osni Marques, LBNL/NERSC, USA \n
 264 *
 265 *  =====================================================================
 266       SUBROUTINE CLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
 267      $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
 268      $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK,
 269      $                   IWORK, INFO )
 270 *
 271 *  -- LAPACK computational routine (version 3.4.2) --
 272 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 273 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 274 *     September 2012
 275 *
 276 *     .. Scalar Arguments ..
 277       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
 278      $                   SMLSIZ
 279 *     ..
 280 *     .. Array Arguments ..
 281       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
 282      $                   K( * ), PERM( LDGCOL, * )
 283       REAL               C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
 284      $                   GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
 285      $                   S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
 286       COMPLEX            B( LDB, * ), BX( LDBX, * )
 287 *     ..
 288 *
 289 *  =====================================================================
 290 *
 291 *     .. Parameters ..
 292       REAL               ZERO, ONE
 293       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 294 *     ..
 295 *     .. Local Scalars ..
 296       INTEGER            I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
 297      $                   JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
 298      $                   NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
 299 *     ..
 300 *     .. External Subroutines ..
 301       EXTERNAL           CCOPY, CLALS0, SGEMM, SLASDT, XERBLA
 302 *     ..
 303 *     .. Intrinsic Functions ..
 304       INTRINSIC          AIMAG, CMPLX, REAL
 305 *     ..
 306 *     .. Executable Statements ..
 307 *
 308 *     Test the input parameters.
 309 *
 310       INFO = 0
 311 *
 312       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 313          INFO = -1
 314       ELSE IF( SMLSIZ.LT.3 ) THEN
 315          INFO = -2
 316       ELSE IF( N.LT.SMLSIZ ) THEN
 317          INFO = -3
 318       ELSE IF( NRHS.LT.1 ) THEN
 319          INFO = -4
 320       ELSE IF( LDB.LT.N ) THEN
 321          INFO = -6
 322       ELSE IF( LDBX.LT.N ) THEN
 323          INFO = -8
 324       ELSE IF( LDU.LT.N ) THEN
 325          INFO = -10
 326       ELSE IF( LDGCOL.LT.N ) THEN
 327          INFO = -19
 328       END IF
 329       IF( INFO.NE.0 ) THEN
 330          CALL XERBLA( 'CLALSA', -INFO )
 331          RETURN
 332       END IF
 333 *
 334 *     Book-keeping and  setting up the computation tree.
 335 *
 336       INODE = 1
 337       NDIML = INODE + N
 338       NDIMR = NDIML + N
 339 *
 340       CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
 341      $             IWORK( NDIMR ), SMLSIZ )
 342 *
 343 *     The following code applies back the left singular vector factors.
 344 *     For applying back the right singular vector factors, go to 170.
 345 *
 346       IF( ICOMPQ.EQ.1 ) THEN
 347          GO TO 170
 348       END IF
 349 *
 350 *     The nodes on the bottom level of the tree were solved
 351 *     by SLASDQ. The corresponding left and right singular vector
 352 *     matrices are in explicit form. First apply back the left
 353 *     singular vector matrices.
 354 *
 355       NDB1 = ( ND+1 ) / 2
 356       DO 130 I = NDB1, ND
 357 *
 358 *        IC : center row of each node
 359 *        NL : number of rows of left  subproblem
 360 *        NR : number of rows of right subproblem
 361 *        NLF: starting row of the left   subproblem
 362 *        NRF: starting row of the right  subproblem
 363 *
 364          I1 = I - 1
 365          IC = IWORK( INODE+I1 )
 366          NL = IWORK( NDIML+I1 )
 367          NR = IWORK( NDIMR+I1 )
 368          NLF = IC - NL
 369          NRF = IC + 1
 370 *
 371 *        Since B and BX are complex, the following call to SGEMM
 372 *        is performed in two steps (real and imaginary parts).
 373 *
 374 *        CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
 375 *     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
 376 *
 377          J = NL*NRHS*2
 378          DO 20 JCOL = 1, NRHS
 379             DO 10 JROW = NLF, NLF + NL - 1
 380                J = J + 1
 381                RWORK( J ) = REAL( B( JROW, JCOL ) )
 382    10       CONTINUE
 383    20    CONTINUE
 384          CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
 385      $               RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1 ), NL )
 386          J = NL*NRHS*2
 387          DO 40 JCOL = 1, NRHS
 388             DO 30 JROW = NLF, NLF + NL - 1
 389                J = J + 1
 390                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 391    30       CONTINUE
 392    40    CONTINUE
 393          CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
 394      $               RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1+NL*NRHS ),
 395      $               NL )
 396          JREAL = 0
 397          JIMAG = NL*NRHS
 398          DO 60 JCOL = 1, NRHS
 399             DO 50 JROW = NLF, NLF + NL - 1
 400                JREAL = JREAL + 1
 401                JIMAG = JIMAG + 1
 402                BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 403      $                            RWORK( JIMAG ) )
 404    50       CONTINUE
 405    60    CONTINUE
 406 *
 407 *        Since B and BX are complex, the following call to SGEMM
 408 *        is performed in two steps (real and imaginary parts).
 409 *
 410 *        CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
 411 *    $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
 412 *
 413          J = NR*NRHS*2
 414          DO 80 JCOL = 1, NRHS
 415             DO 70 JROW = NRF, NRF + NR - 1
 416                J = J + 1
 417                RWORK( J ) = REAL( B( JROW, JCOL ) )
 418    70       CONTINUE
 419    80    CONTINUE
 420          CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
 421      $               RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1 ), NR )
 422          J = NR*NRHS*2
 423          DO 100 JCOL = 1, NRHS
 424             DO 90 JROW = NRF, NRF + NR - 1
 425                J = J + 1
 426                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 427    90       CONTINUE
 428   100    CONTINUE
 429          CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
 430      $               RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1+NR*NRHS ),
 431      $               NR )
 432          JREAL = 0
 433          JIMAG = NR*NRHS
 434          DO 120 JCOL = 1, NRHS
 435             DO 110 JROW = NRF, NRF + NR - 1
 436                JREAL = JREAL + 1
 437                JIMAG = JIMAG + 1
 438                BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 439      $                            RWORK( JIMAG ) )
 440   110       CONTINUE
 441   120    CONTINUE
 442 *
 443   130 CONTINUE
 444 *
 445 *     Next copy the rows of B that correspond to unchanged rows
 446 *     in the bidiagonal matrix to BX.
 447 *
 448       DO 140 I = 1, ND
 449          IC = IWORK( INODE+I-1 )
 450          CALL CCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
 451   140 CONTINUE
 452 *
 453 *     Finally go through the left singular vector matrices of all
 454 *     the other subproblems bottom-up on the tree.
 455 *
 456       J = 2**NLVL
 457       SQRE = 0
 458 *
 459       DO 160 LVL = NLVL, 1, -1
 460          LVL2 = 2*LVL - 1
 461 *
 462 *        find the first node LF and last node LL on
 463 *        the current level LVL
 464 *
 465          IF( LVL.EQ.1 ) THEN
 466             LF = 1
 467             LL = 1
 468          ELSE
 469             LF = 2**( LVL-1 )
 470             LL = 2*LF - 1
 471          END IF
 472          DO 150 I = LF, LL
 473             IM1 = I - 1
 474             IC = IWORK( INODE+IM1 )
 475             NL = IWORK( NDIML+IM1 )
 476             NR = IWORK( NDIMR+IM1 )
 477             NLF = IC - NL
 478             NRF = IC + 1
 479             J = J - 1
 480             CALL CLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
 481      $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
 482      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
 483      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
 484      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
 485      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
 486      $                   INFO )
 487   150    CONTINUE
 488   160 CONTINUE
 489       GO TO 330
 490 *
 491 *     ICOMPQ = 1: applying back the right singular vector factors.
 492 *
 493   170 CONTINUE
 494 *
 495 *     First now go through the right singular vector matrices of all
 496 *     the tree nodes top-down.
 497 *
 498       J = 0
 499       DO 190 LVL = 1, NLVL
 500          LVL2 = 2*LVL - 1
 501 *
 502 *        Find the first node LF and last node LL on
 503 *        the current level LVL.
 504 *
 505          IF( LVL.EQ.1 ) THEN
 506             LF = 1
 507             LL = 1
 508          ELSE
 509             LF = 2**( LVL-1 )
 510             LL = 2*LF - 1
 511          END IF
 512          DO 180 I = LL, LF, -1
 513             IM1 = I - 1
 514             IC = IWORK( INODE+IM1 )
 515             NL = IWORK( NDIML+IM1 )
 516             NR = IWORK( NDIMR+IM1 )
 517             NLF = IC - NL
 518             NRF = IC + 1
 519             IF( I.EQ.LL ) THEN
 520                SQRE = 0
 521             ELSE
 522                SQRE = 1
 523             END IF
 524             J = J + 1
 525             CALL CLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
 526      $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
 527      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
 528      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
 529      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
 530      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK,
 531      $                   INFO )
 532   180    CONTINUE
 533   190 CONTINUE
 534 *
 535 *     The nodes on the bottom level of the tree were solved
 536 *     by SLASDQ. The corresponding right singular vector
 537 *     matrices are in explicit form. Apply them back.
 538 *
 539       NDB1 = ( ND+1 ) / 2
 540       DO 320 I = NDB1, ND
 541          I1 = I - 1
 542          IC = IWORK( INODE+I1 )
 543          NL = IWORK( NDIML+I1 )
 544          NR = IWORK( NDIMR+I1 )
 545          NLP1 = NL + 1
 546          IF( I.EQ.ND ) THEN
 547             NRP1 = NR
 548          ELSE
 549             NRP1 = NR + 1
 550          END IF
 551          NLF = IC - NL
 552          NRF = IC + 1
 553 *
 554 *        Since B and BX are complex, the following call to SGEMM is
 555 *        performed in two steps (real and imaginary parts).
 556 *
 557 *        CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
 558 *    $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
 559 *
 560          J = NLP1*NRHS*2
 561          DO 210 JCOL = 1, NRHS
 562             DO 200 JROW = NLF, NLF + NLP1 - 1
 563                J = J + 1
 564                RWORK( J ) = REAL( B( JROW, JCOL ) )
 565   200       CONTINUE
 566   210    CONTINUE
 567          CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
 568      $               RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO, RWORK( 1 ),
 569      $               NLP1 )
 570          J = NLP1*NRHS*2
 571          DO 230 JCOL = 1, NRHS
 572             DO 220 JROW = NLF, NLF + NLP1 - 1
 573                J = J + 1
 574                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 575   220       CONTINUE
 576   230    CONTINUE
 577          CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
 578      $               RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO,
 579      $               RWORK( 1+NLP1*NRHS ), NLP1 )
 580          JREAL = 0
 581          JIMAG = NLP1*NRHS
 582          DO 250 JCOL = 1, NRHS
 583             DO 240 JROW = NLF, NLF + NLP1 - 1
 584                JREAL = JREAL + 1
 585                JIMAG = JIMAG + 1
 586                BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 587      $                            RWORK( JIMAG ) )
 588   240       CONTINUE
 589   250    CONTINUE
 590 *
 591 *        Since B and BX are complex, the following call to SGEMM is
 592 *        performed in two steps (real and imaginary parts).
 593 *
 594 *        CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
 595 *    $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
 596 *
 597          J = NRP1*NRHS*2
 598          DO 270 JCOL = 1, NRHS
 599             DO 260 JROW = NRF, NRF + NRP1 - 1
 600                J = J + 1
 601                RWORK( J ) = REAL( B( JROW, JCOL ) )
 602   260       CONTINUE
 603   270    CONTINUE
 604          CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
 605      $               RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO, RWORK( 1 ),
 606      $               NRP1 )
 607          J = NRP1*NRHS*2
 608          DO 290 JCOL = 1, NRHS
 609             DO 280 JROW = NRF, NRF + NRP1 - 1
 610                J = J + 1
 611                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 612   280       CONTINUE
 613   290    CONTINUE
 614          CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
 615      $               RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO,
 616      $               RWORK( 1+NRP1*NRHS ), NRP1 )
 617          JREAL = 0
 618          JIMAG = NRP1*NRHS
 619          DO 310 JCOL = 1, NRHS
 620             DO 300 JROW = NRF, NRF + NRP1 - 1
 621                JREAL = JREAL + 1
 622                JIMAG = JIMAG + 1
 623                BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 624      $                            RWORK( JIMAG ) )
 625   300       CONTINUE
 626   310    CONTINUE
 627 *
 628   320 CONTINUE
 629 *
 630   330 CONTINUE
 631 *
 632       RETURN
 633 *
 634 *     End of CLALSA
 635 *
 636       END