SRC/clalsd.f

   1 *> \brief \b CLALSD uses the singular value decomposition of A to solve the least squares problem.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CLALSD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
  22 *                          RANK, WORK, RWORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
  27 *       REAL               RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * )
  31 *       REAL               D( * ), E( * ), RWORK( * )
  32 *       COMPLEX            B( LDB, * ), WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> CLALSD uses the singular value decomposition of A to solve the least
  42 *> squares problem of finding X to minimize the Euclidean norm of each
  43 *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
  44 *> are N-by-NRHS. The solution X overwrites B.
  45 *>
  46 *> The singular values of A smaller than RCOND times the largest
  47 *> singular value are treated as zero in solving the least squares
  48 *> problem; in this case a minimum norm solution is returned.
  49 *> The actual singular values are returned in D in ascending order.
  50 *>
  51 *> This code makes very mild assumptions about floating point
  52 *> arithmetic. It will work on machines with a guard digit in
  53 *> add/subtract, or on those binary machines without guard digits
  54 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
  55 *> It could conceivably fail on hexadecimal or decimal machines
  56 *> without guard digits, but we know of none.
  57 *> \endverbatim
  58 *
  59 *  Arguments:
  60 *  ==========
  61 *
  62 *> \param[in] UPLO
  63 *> \verbatim
  64 *>          UPLO is CHARACTER*1
  65 *>         = 'U': D and E define an upper bidiagonal matrix.
  66 *>         = 'L': D and E define a  lower bidiagonal 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 dimension of the  bidiagonal matrix.  N >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] NRHS
  83 *> \verbatim
  84 *>          NRHS is INTEGER
  85 *>         The number of columns of B. NRHS must be at least 1.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] D
  89 *> \verbatim
  90 *>          D is REAL array, dimension (N)
  91 *>         On entry D contains the main diagonal of the bidiagonal
  92 *>         matrix. On exit, if INFO = 0, D contains its singular values.
  93 *> \endverbatim
  94 *>
  95 *> \param[in,out] E
  96 *> \verbatim
  97 *>          E is REAL array, dimension (N-1)
  98 *>         Contains the super-diagonal entries of the bidiagonal matrix.
  99 *>         On exit, E has been destroyed.
 100 *> \endverbatim
 101 *>
 102 *> \param[in,out] B
 103 *> \verbatim
 104 *>          B is COMPLEX array, dimension (LDB,NRHS)
 105 *>         On input, B contains the right hand sides of the least
 106 *>         squares problem. On output, B contains the solution X.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LDB
 110 *> \verbatim
 111 *>          LDB is INTEGER
 112 *>         The leading dimension of B in the calling subprogram.
 113 *>         LDB must be at least max(1,N).
 114 *> \endverbatim
 115 *>
 116 *> \param[in] RCOND
 117 *> \verbatim
 118 *>          RCOND is REAL
 119 *>         The singular values of A less than or equal to RCOND times
 120 *>         the largest singular value are treated as zero in solving
 121 *>         the least squares problem. If RCOND is negative,
 122 *>         machine precision is used instead.
 123 *>         For example, if diag(S)*X=B were the least squares problem,
 124 *>         where diag(S) is a diagonal matrix of singular values, the
 125 *>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
 126 *>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
 127 *>         RCOND*max(S).
 128 *> \endverbatim
 129 *>
 130 *> \param[out] RANK
 131 *> \verbatim
 132 *>          RANK is INTEGER
 133 *>         The number of singular values of A greater than RCOND times
 134 *>         the largest singular value.
 135 *> \endverbatim
 136 *>
 137 *> \param[out] WORK
 138 *> \verbatim
 139 *>          WORK is COMPLEX array, dimension (N * NRHS).
 140 *> \endverbatim
 141 *>
 142 *> \param[out] RWORK
 143 *> \verbatim
 144 *>          RWORK is REAL array, dimension at least
 145 *>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
 146 *>         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
 147 *>         where
 148 *>         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
 149 *> \endverbatim
 150 *>
 151 *> \param[out] IWORK
 152 *> \verbatim
 153 *>          IWORK is INTEGER array, dimension (3*N*NLVL + 11*N).
 154 *> \endverbatim
 155 *>
 156 *> \param[out] INFO
 157 *> \verbatim
 158 *>          INFO is INTEGER
 159 *>         = 0:  successful exit.
 160 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
 161 *>         > 0:  The algorithm failed to compute a singular value while
 162 *>               working on the submatrix lying in rows and columns
 163 *>               INFO/(N+1) through MOD(INFO,N+1).
 164 *> \endverbatim
 165 *
 166 *  Authors:
 167 *  ========
 168 *
 169 *> \author Univ. of Tennessee
 170 *> \author Univ. of California Berkeley
 171 *> \author Univ. of Colorado Denver
 172 *> \author NAG Ltd.
 173 *
 174 *> \date September 2012
 175 *
 176 *> \ingroup complexOTHERcomputational
 177 *
 178 *> \par Contributors:
 179 *  ==================
 180 *>
 181 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
 182 *>       California at Berkeley, USA \n
 183 *>     Osni Marques, LBNL/NERSC, USA \n
 184 *
 185 *  =====================================================================
 186       SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
 187      $                   RANK, WORK, RWORK, IWORK, INFO )
 188 *
 189 *  -- LAPACK computational routine (version 3.4.2) --
 190 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 191 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 192 *     September 2012
 193 *
 194 *     .. Scalar Arguments ..
 195       CHARACTER          UPLO
 196       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
 197       REAL               RCOND
 198 *     ..
 199 *     .. Array Arguments ..
 200       INTEGER            IWORK( * )
 201       REAL               D( * ), E( * ), RWORK( * )
 202       COMPLEX            B( LDB, * ), WORK( * )
 203 *     ..
 204 *
 205 *  =====================================================================
 206 *
 207 *     .. Parameters ..
 208       REAL               ZERO, ONE, TWO
 209       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
 210       COMPLEX            CZERO
 211       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
 212 *     ..
 213 *     .. Local Scalars ..
 214       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
 215      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
 216      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
 217      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
 218      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
 219      $                   U, VT, Z
 220       REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
 221 *     ..
 222 *     .. External Functions ..
 223       INTEGER            ISAMAX
 224       REAL               SLAMCH, SLANST
 225       EXTERNAL           ISAMAX, SLAMCH, SLANST
 226 *     ..
 227 *     .. External Subroutines ..
 228       EXTERNAL           CCOPY, CLACPY, CLALSA, CLASCL, CLASET, CSROT,
 229      $                   SGEMM, SLARTG, SLASCL, SLASDA, SLASDQ, SLASET,
 230      $                   SLASRT, XERBLA
 231 *     ..
 232 *     .. Intrinsic Functions ..
 233       INTRINSIC          ABS, AIMAG, CMPLX, INT, LOG, REAL, SIGN
 234 *     ..
 235 *     .. Executable Statements ..
 236 *
 237 *     Test the input parameters.
 238 *
 239       INFO = 0
 240 *
 241       IF( N.LT.0 ) THEN
 242          INFO = -3
 243       ELSE IF( NRHS.LT.1 ) THEN
 244          INFO = -4
 245       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
 246          INFO = -8
 247       END IF
 248       IF( INFO.NE.0 ) THEN
 249          CALL XERBLA( 'CLALSD', -INFO )
 250          RETURN
 251       END IF
 252 *
 253       EPS = SLAMCH( 'Epsilon' )
 254 *
 255 *     Set up the tolerance.
 256 *
 257       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
 258          RCND = EPS
 259       ELSE
 260          RCND = RCOND
 261       END IF
 262 *
 263       RANK = 0
 264 *
 265 *     Quick return if possible.
 266 *
 267       IF( N.EQ.0 ) THEN
 268          RETURN
 269       ELSE IF( N.EQ.1 ) THEN
 270          IF( D( 1 ).EQ.ZERO ) THEN
 271             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
 272          ELSE
 273             RANK = 1
 274             CALL CLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
 275             D( 1 ) = ABS( D( 1 ) )
 276          END IF
 277          RETURN
 278       END IF
 279 *
 280 *     Rotate the matrix if it is lower bidiagonal.
 281 *
 282       IF( UPLO.EQ.'L' ) THEN
 283          DO 10 I = 1, N - 1
 284             CALL SLARTG( D( I ), E( I ), CS, SN, R )
 285             D( I ) = R
 286             E( I ) = SN*D( I+1 )
 287             D( I+1 ) = CS*D( I+1 )
 288             IF( NRHS.EQ.1 ) THEN
 289                CALL CSROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
 290             ELSE
 291                RWORK( I*2-1 ) = CS
 292                RWORK( I*2 ) = SN
 293             END IF
 294    10    CONTINUE
 295          IF( NRHS.GT.1 ) THEN
 296             DO 30 I = 1, NRHS
 297                DO 20 J = 1, N - 1
 298                   CS = RWORK( J*2-1 )
 299                   SN = RWORK( J*2 )
 300                   CALL CSROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
 301    20          CONTINUE
 302    30       CONTINUE
 303          END IF
 304       END IF
 305 *
 306 *     Scale.
 307 *
 308       NM1 = N - 1
 309       ORGNRM = SLANST( 'M', N, D, E )
 310       IF( ORGNRM.EQ.ZERO ) THEN
 311          CALL CLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
 312          RETURN
 313       END IF
 314 *
 315       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
 316       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
 317 *
 318 *     If N is smaller than the minimum divide size SMLSIZ, then solve
 319 *     the problem with another solver.
 320 *
 321       IF( N.LE.SMLSIZ ) THEN
 322          IRWU = 1
 323          IRWVT = IRWU + N*N
 324          IRWWRK = IRWVT + N*N
 325          IRWRB = IRWWRK
 326          IRWIB = IRWRB + N*NRHS
 327          IRWB = IRWIB + N*NRHS
 328          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
 329          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
 330          CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
 331      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
 332      $                RWORK( IRWWRK ), INFO )
 333          IF( INFO.NE.0 ) THEN
 334             RETURN
 335          END IF
 336 *
 337 *        In the real version, B is passed to SLASDQ and multiplied
 338 *        internally by Q**H. Here B is complex and that product is
 339 *        computed below in two steps (real and imaginary parts).
 340 *
 341          J = IRWB - 1
 342          DO 50 JCOL = 1, NRHS
 343             DO 40 JROW = 1, N
 344                J = J + 1
 345                RWORK( J ) = REAL( B( JROW, JCOL ) )
 346    40       CONTINUE
 347    50    CONTINUE
 348          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
 349      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
 350          J = IRWB - 1
 351          DO 70 JCOL = 1, NRHS
 352             DO 60 JROW = 1, N
 353                J = J + 1
 354                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 355    60       CONTINUE
 356    70    CONTINUE
 357          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
 358      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
 359          JREAL = IRWRB - 1
 360          JIMAG = IRWIB - 1
 361          DO 90 JCOL = 1, NRHS
 362             DO 80 JROW = 1, N
 363                JREAL = JREAL + 1
 364                JIMAG = JIMAG + 1
 365                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
 366    80       CONTINUE
 367    90    CONTINUE
 368 *
 369          TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
 370          DO 100 I = 1, N
 371             IF( D( I ).LE.TOL ) THEN
 372                CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
 373             ELSE
 374                CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
 375      $                      LDB, INFO )
 376                RANK = RANK + 1
 377             END IF
 378   100    CONTINUE
 379 *
 380 *        Since B is complex, the following call to SGEMM is performed
 381 *        in two steps (real and imaginary parts). That is for V * B
 382 *        (in the real version of the code V**H is stored in WORK).
 383 *
 384 *        CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
 385 *    $               WORK( NWORK ), N )
 386 *
 387          J = IRWB - 1
 388          DO 120 JCOL = 1, NRHS
 389             DO 110 JROW = 1, N
 390                J = J + 1
 391                RWORK( J ) = REAL( B( JROW, JCOL ) )
 392   110       CONTINUE
 393   120    CONTINUE
 394          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
 395      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
 396          J = IRWB - 1
 397          DO 140 JCOL = 1, NRHS
 398             DO 130 JROW = 1, N
 399                J = J + 1
 400                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 401   130       CONTINUE
 402   140    CONTINUE
 403          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
 404      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
 405          JREAL = IRWRB - 1
 406          JIMAG = IRWIB - 1
 407          DO 160 JCOL = 1, NRHS
 408             DO 150 JROW = 1, N
 409                JREAL = JREAL + 1
 410                JIMAG = JIMAG + 1
 411                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
 412   150       CONTINUE
 413   160    CONTINUE
 414 *
 415 *        Unscale.
 416 *
 417          CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 418          CALL SLASRT( 'D', N, D, INFO )
 419          CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
 420 *
 421          RETURN
 422       END IF
 423 *
 424 *     Book-keeping and setting up some constants.
 425 *
 426       NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
 427 *
 428       SMLSZP = SMLSIZ + 1
 429 *
 430       U = 1
 431       VT = 1 + SMLSIZ*N
 432       DIFL = VT + SMLSZP*N
 433       DIFR = DIFL + NLVL*N
 434       Z = DIFR + NLVL*N*2
 435       C = Z + NLVL*N
 436       S = C + N
 437       POLES = S + N
 438       GIVNUM = POLES + 2*NLVL*N
 439       NRWORK = GIVNUM + 2*NLVL*N
 440       BX = 1
 441 *
 442       IRWRB = NRWORK
 443       IRWIB = IRWRB + SMLSIZ*NRHS
 444       IRWB = IRWIB + SMLSIZ*NRHS
 445 *
 446       SIZEI = 1 + N
 447       K = SIZEI + N
 448       GIVPTR = K + N
 449       PERM = GIVPTR + N
 450       GIVCOL = PERM + NLVL*N
 451       IWK = GIVCOL + NLVL*N*2
 452 *
 453       ST = 1
 454       SQRE = 0
 455       ICMPQ1 = 1
 456       ICMPQ2 = 0
 457       NSUB = 0
 458 *
 459       DO 170 I = 1, N
 460          IF( ABS( D( I ) ).LT.EPS ) THEN
 461             D( I ) = SIGN( EPS, D( I ) )
 462          END IF
 463   170 CONTINUE
 464 *
 465       DO 240 I = 1, NM1
 466          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
 467             NSUB = NSUB + 1
 468             IWORK( NSUB ) = ST
 469 *
 470 *           Subproblem found. First determine its size and then
 471 *           apply divide and conquer on it.
 472 *
 473             IF( I.LT.NM1 ) THEN
 474 *
 475 *              A subproblem with E(I) small for I < NM1.
 476 *
 477                NSIZE = I - ST + 1
 478                IWORK( SIZEI+NSUB-1 ) = NSIZE
 479             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
 480 *
 481 *              A subproblem with E(NM1) not too small but I = NM1.
 482 *
 483                NSIZE = N - ST + 1
 484                IWORK( SIZEI+NSUB-1 ) = NSIZE
 485             ELSE
 486 *
 487 *              A subproblem with E(NM1) small. This implies an
 488 *              1-by-1 subproblem at D(N), which is not solved
 489 *              explicitly.
 490 *
 491                NSIZE = I - ST + 1
 492                IWORK( SIZEI+NSUB-1 ) = NSIZE
 493                NSUB = NSUB + 1
 494                IWORK( NSUB ) = N
 495                IWORK( SIZEI+NSUB-1 ) = 1
 496                CALL CCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
 497             END IF
 498             ST1 = ST - 1
 499             IF( NSIZE.EQ.1 ) THEN
 500 *
 501 *              This is a 1-by-1 subproblem and is not solved
 502 *              explicitly.
 503 *
 504                CALL CCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
 505             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
 506 *
 507 *              This is a small subproblem and is solved by SLASDQ.
 508 *
 509                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
 510      $                      RWORK( VT+ST1 ), N )
 511                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
 512      $                      RWORK( U+ST1 ), N )
 513                CALL SLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
 514      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
 515      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
 516      $                      INFO )
 517                IF( INFO.NE.0 ) THEN
 518                   RETURN
 519                END IF
 520 *
 521 *              In the real version, B is passed to SLASDQ and multiplied
 522 *              internally by Q**H. Here B is complex and that product is
 523 *              computed below in two steps (real and imaginary parts).
 524 *
 525                J = IRWB - 1
 526                DO 190 JCOL = 1, NRHS
 527                   DO 180 JROW = ST, ST + NSIZE - 1
 528                      J = J + 1
 529                      RWORK( J ) = REAL( B( JROW, JCOL ) )
 530   180             CONTINUE
 531   190          CONTINUE
 532                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 533      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
 534      $                     ZERO, RWORK( IRWRB ), NSIZE )
 535                J = IRWB - 1
 536                DO 210 JCOL = 1, NRHS
 537                   DO 200 JROW = ST, ST + NSIZE - 1
 538                      J = J + 1
 539                      RWORK( J ) = AIMAG( B( JROW, JCOL ) )
 540   200             CONTINUE
 541   210          CONTINUE
 542                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 543      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
 544      $                     ZERO, RWORK( IRWIB ), NSIZE )
 545                JREAL = IRWRB - 1
 546                JIMAG = IRWIB - 1
 547                DO 230 JCOL = 1, NRHS
 548                   DO 220 JROW = ST, ST + NSIZE - 1
 549                      JREAL = JREAL + 1
 550                      JIMAG = JIMAG + 1
 551                      B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 552      $                                 RWORK( JIMAG ) )
 553   220             CONTINUE
 554   230          CONTINUE
 555 *
 556                CALL CLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
 557      $                      WORK( BX+ST1 ), N )
 558             ELSE
 559 *
 560 *              A large problem. Solve it using divide and conquer.
 561 *
 562                CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
 563      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
 564      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
 565      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
 566      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
 567      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
 568      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
 569      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
 570      $                      IWORK( IWK ), INFO )
 571                IF( INFO.NE.0 ) THEN
 572                   RETURN
 573                END IF
 574                BXST = BX + ST1
 575                CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
 576      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
 577      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
 578      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
 579      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
 580      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
 581      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
 582      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
 583      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
 584                IF( INFO.NE.0 ) THEN
 585                   RETURN
 586                END IF
 587             END IF
 588             ST = I + 1
 589          END IF
 590   240 CONTINUE
 591 *
 592 *     Apply the singular values and treat the tiny ones as zero.
 593 *
 594       TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
 595 *
 596       DO 250 I = 1, N
 597 *
 598 *        Some of the elements in D can be negative because 1-by-1
 599 *        subproblems were not solved explicitly.
 600 *
 601          IF( ABS( D( I ) ).LE.TOL ) THEN
 602             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
 603          ELSE
 604             RANK = RANK + 1
 605             CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
 606      $                   WORK( BX+I-1 ), N, INFO )
 607          END IF
 608          D( I ) = ABS( D( I ) )
 609   250 CONTINUE
 610 *
 611 *     Now apply back the right singular vectors.
 612 *
 613       ICMPQ2 = 1
 614       DO 320 I = 1, NSUB
 615          ST = IWORK( I )
 616          ST1 = ST - 1
 617          NSIZE = IWORK( SIZEI+I-1 )
 618          BXST = BX + ST1
 619          IF( NSIZE.EQ.1 ) THEN
 620             CALL CCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
 621          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
 622 *
 623 *           Since B and BX are complex, the following call to SGEMM
 624 *           is performed in two steps (real and imaginary parts).
 625 *
 626 *           CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 627 *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
 628 *    $                  B( ST, 1 ), LDB )
 629 *
 630             J = BXST - N - 1
 631             JREAL = IRWB - 1
 632             DO 270 JCOL = 1, NRHS
 633                J = J + N
 634                DO 260 JROW = 1, NSIZE
 635                   JREAL = JREAL + 1
 636                   RWORK( JREAL ) = REAL( WORK( J+JROW ) )
 637   260          CONTINUE
 638   270       CONTINUE
 639             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 640      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
 641      $                  RWORK( IRWRB ), NSIZE )
 642             J = BXST - N - 1
 643             JIMAG = IRWB - 1
 644             DO 290 JCOL = 1, NRHS
 645                J = J + N
 646                DO 280 JROW = 1, NSIZE
 647                   JIMAG = JIMAG + 1
 648                   RWORK( JIMAG ) = AIMAG( WORK( J+JROW ) )
 649   280          CONTINUE
 650   290       CONTINUE
 651             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 652      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
 653      $                  RWORK( IRWIB ), NSIZE )
 654             JREAL = IRWRB - 1
 655             JIMAG = IRWIB - 1
 656             DO 310 JCOL = 1, NRHS
 657                DO 300 JROW = ST, ST + NSIZE - 1
 658                   JREAL = JREAL + 1
 659                   JIMAG = JIMAG + 1
 660                   B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
 661      $                              RWORK( JIMAG ) )
 662   300          CONTINUE
 663   310       CONTINUE
 664          ELSE
 665             CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
 666      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
 667      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
 668      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
 669      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
 670      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
 671      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
 672      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
 673      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
 674             IF( INFO.NE.0 ) THEN
 675                RETURN
 676             END IF
 677          END IF
 678   320 CONTINUE
 679 *
 680 *     Unscale and sort the singular values.
 681 *
 682       CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 683       CALL SLASRT( 'D', N, D, INFO )
 684       CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
 685 *
 686       RETURN
 687 *
 688 *     End of CLALSD
 689 *
 690       END