SRC/zlalsd.f

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