LAPACK: SRC/zlalsd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
00002      $                   RANK, WORK, RWORK, IWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          UPLO
00011       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
00012       DOUBLE PRECISION   RCOND
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IWORK( * )
00016       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
00017       COMPLEX*16         B( LDB, * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  ZLALSD uses the singular value decomposition of A to solve the least
00024 *  squares problem of finding X to minimize the Euclidean norm of each
00025 *  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
00026 *  are N-by-NRHS. The solution X overwrites B.
00027 *
00028 *  The singular values of A smaller than RCOND times the largest
00029 *  singular value are treated as zero in solving the least squares
00030 *  problem; in this case a minimum norm solution is returned.
00031 *  The actual singular values are returned in D in ascending order.
00032 *
00033 *  This code makes very mild assumptions about floating point
00034 *  arithmetic. It will work on machines with a guard digit in
00035 *  add/subtract, or on those binary machines without guard digits
00036 *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
00037 *  It could conceivably fail on hexadecimal or decimal machines
00038 *  without guard digits, but we know of none.
00039 *
00040 *  Arguments
00041 *  =========
00042 *
00043 *  UPLO   (input) CHARACTER*1
00044 *         = 'U': D and E define an upper bidiagonal matrix.
00045 *         = 'L': D and E define a  lower bidiagonal matrix.
00046 *
00047 *  SMLSIZ (input) INTEGER
00048 *         The maximum size of the subproblems at the bottom of the
00049 *         computation tree.
00050 *
00051 *  N      (input) INTEGER
00052 *         The dimension of the  bidiagonal matrix.  N >= 0.
00053 *
00054 *  NRHS   (input) INTEGER
00055 *         The number of columns of B. NRHS must be at least 1.
00056 *
00057 *  D      (input/output) DOUBLE PRECISION array, dimension (N)
00058 *         On entry D contains the main diagonal of the bidiagonal
00059 *         matrix. On exit, if INFO = 0, D contains its singular values.
00060 *
00061 *  E      (input/output) DOUBLE PRECISION array, dimension (N-1)
00062 *         Contains the super-diagonal entries of the bidiagonal matrix.
00063 *         On exit, E has been destroyed.
00064 *
00065 *  B      (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
00066 *         On input, B contains the right hand sides of the least
00067 *         squares problem. On output, B contains the solution X.
00068 *
00069 *  LDB    (input) INTEGER
00070 *         The leading dimension of B in the calling subprogram.
00071 *         LDB must be at least max(1,N).
00072 *
00073 *  RCOND  (input) DOUBLE PRECISION
00074 *         The singular values of A less than or equal to RCOND times
00075 *         the largest singular value are treated as zero in solving
00076 *         the least squares problem. If RCOND is negative,
00077 *         machine precision is used instead.
00078 *         For example, if diag(S)*X=B were the least squares problem,
00079 *         where diag(S) is a diagonal matrix of singular values, the
00080 *         solution would be X(i) = B(i) / S(i) if S(i) is greater than
00081 *         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
00082 *         RCOND*max(S).
00083 *
00084 *  RANK   (output) INTEGER
00085 *         The number of singular values of A greater than RCOND times
00086 *         the largest singular value.
00087 *
00088 *  WORK   (workspace) COMPLEX*16 array, dimension at least
00089 *         (N * NRHS).
00090 *
00091 *  RWORK  (workspace) DOUBLE PRECISION array, dimension at least
00092 *         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
00093 *         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
00094 *         where
00095 *         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
00096 *
00097 *  IWORK  (workspace) INTEGER array, dimension at least
00098 *         (3*N*NLVL + 11*N).
00099 *
00100 *  INFO   (output) INTEGER
00101 *         = 0:  successful exit.
00102 *         < 0:  if INFO = -i, the i-th argument had an illegal value.
00103 *         > 0:  The algorithm failed to compute a singular value while
00104 *               working on the submatrix lying in rows and columns
00105 *               INFO/(N+1) through MOD(INFO,N+1).
00106 *
00107 *  Further Details
00108 *  ===============
00109 *
00110 *  Based on contributions by
00111 *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00112 *       California at Berkeley, USA
00113 *     Osni Marques, LBNL/NERSC, USA
00114 *
00115 *  =====================================================================
00116 *
00117 *     .. Parameters ..
00118       DOUBLE PRECISION   ZERO, ONE, TWO
00119       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00120       COMPLEX*16         CZERO
00121       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
00122 *     ..
00123 *     .. Local Scalars ..
00124       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
00125      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
00126      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
00127      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
00128      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
00129      $                   U, VT, Z
00130       DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
00131 *     ..
00132 *     .. External Functions ..
00133       INTEGER            IDAMAX
00134       DOUBLE PRECISION   DLAMCH, DLANST
00135       EXTERNAL           IDAMAX, DLAMCH, DLANST
00136 *     ..
00137 *     .. External Subroutines ..
00138       EXTERNAL           DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
00139      $                   DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
00140      $                   ZLASCL, ZLASET
00141 *     ..
00142 *     .. Intrinsic Functions ..
00143       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
00144 *     ..
00145 *     .. Executable Statements ..
00146 *
00147 *     Test the input parameters.
00148 *
00149       INFO = 0
00150 *
00151       IF( N.LT.0 ) THEN
00152          INFO = -3
00153       ELSE IF( NRHS.LT.1 ) THEN
00154          INFO = -4
00155       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
00156          INFO = -8
00157       END IF
00158       IF( INFO.NE.0 ) THEN
00159          CALL XERBLA( 'ZLALSD', -INFO )
00160          RETURN
00161       END IF
00162 *
00163       EPS = DLAMCH( 'Epsilon' )
00164 *
00165 *     Set up the tolerance.
00166 *
00167       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
00168          RCND = EPS
00169       ELSE
00170          RCND = RCOND
00171       END IF
00172 *
00173       RANK = 0
00174 *
00175 *     Quick return if possible.
00176 *
00177       IF( N.EQ.0 ) THEN
00178          RETURN
00179       ELSE IF( N.EQ.1 ) THEN
00180          IF( D( 1 ).EQ.ZERO ) THEN
00181             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
00182          ELSE
00183             RANK = 1
00184             CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
00185             D( 1 ) = ABS( D( 1 ) )
00186          END IF
00187          RETURN
00188       END IF
00189 *
00190 *     Rotate the matrix if it is lower bidiagonal.
00191 *
00192       IF( UPLO.EQ.'L' ) THEN
00193          DO 10 I = 1, N - 1
00194             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00195             D( I ) = R
00196             E( I ) = SN*D( I+1 )
00197             D( I+1 ) = CS*D( I+1 )
00198             IF( NRHS.EQ.1 ) THEN
00199                CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
00200             ELSE
00201                RWORK( I*2-1 ) = CS
00202                RWORK( I*2 ) = SN
00203             END IF
00204    10    CONTINUE
00205          IF( NRHS.GT.1 ) THEN
00206             DO 30 I = 1, NRHS
00207                DO 20 J = 1, N - 1
00208                   CS = RWORK( J*2-1 )
00209                   SN = RWORK( J*2 )
00210                   CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
00211    20          CONTINUE
00212    30       CONTINUE
00213          END IF
00214       END IF
00215 *
00216 *     Scale.
00217 *
00218       NM1 = N - 1
00219       ORGNRM = DLANST( 'M', N, D, E )
00220       IF( ORGNRM.EQ.ZERO ) THEN
00221          CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
00222          RETURN
00223       END IF
00224 *
00225       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00226       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
00227 *
00228 *     If N is smaller than the minimum divide size SMLSIZ, then solve
00229 *     the problem with another solver.
00230 *
00231       IF( N.LE.SMLSIZ ) THEN
00232          IRWU = 1
00233          IRWVT = IRWU + N*N
00234          IRWWRK = IRWVT + N*N
00235          IRWRB = IRWWRK
00236          IRWIB = IRWRB + N*NRHS
00237          IRWB = IRWIB + N*NRHS
00238          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
00239          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
00240          CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
00241      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
00242      $                RWORK( IRWWRK ), INFO )
00243          IF( INFO.NE.0 ) THEN
00244             RETURN
00245          END IF
00246 *
00247 *        In the real version, B is passed to DLASDQ and multiplied
00248 *        internally by Q'. Here B is complex and that product is
00249 *        computed below in two steps (real and imaginary parts).
00250 *
00251          J = IRWB - 1
00252          DO 50 JCOL = 1, NRHS
00253             DO 40 JROW = 1, N
00254                J = J + 1
00255                RWORK( J ) = DBLE( B( JROW, JCOL ) )
00256    40       CONTINUE
00257    50    CONTINUE
00258          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00259      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00260          J = IRWB - 1
00261          DO 70 JCOL = 1, NRHS
00262             DO 60 JROW = 1, N
00263                J = J + 1
00264                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00265    60       CONTINUE
00266    70    CONTINUE
00267          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00268      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00269          JREAL = IRWRB - 1
00270          JIMAG = IRWIB - 1
00271          DO 90 JCOL = 1, NRHS
00272             DO 80 JROW = 1, N
00273                JREAL = JREAL + 1
00274                JIMAG = JIMAG + 1
00275                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00276      $                           RWORK( JIMAG ) )
00277    80       CONTINUE
00278    90    CONTINUE
00279 *
00280          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00281          DO 100 I = 1, N
00282             IF( D( I ).LE.TOL ) THEN
00283                CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00284             ELSE
00285                CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
00286      $                      LDB, INFO )
00287                RANK = RANK + 1
00288             END IF
00289   100    CONTINUE
00290 *
00291 *        Since B is complex, the following call to DGEMM is performed
00292 *        in two steps (real and imaginary parts). That is for V * B
00293 *        (in the real version of the code V' is stored in WORK).
00294 *
00295 *        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
00296 *    $               WORK( NWORK ), N )
00297 *
00298          J = IRWB - 1
00299          DO 120 JCOL = 1, NRHS
00300             DO 110 JROW = 1, N
00301                J = J + 1
00302                RWORK( J ) = DBLE( B( JROW, JCOL ) )
00303   110       CONTINUE
00304   120    CONTINUE
00305          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00306      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00307          J = IRWB - 1
00308          DO 140 JCOL = 1, NRHS
00309             DO 130 JROW = 1, N
00310                J = J + 1
00311                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00312   130       CONTINUE
00313   140    CONTINUE
00314          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00315      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00316          JREAL = IRWRB - 1
00317          JIMAG = IRWIB - 1
00318          DO 160 JCOL = 1, NRHS
00319             DO 150 JROW = 1, N
00320                JREAL = JREAL + 1
00321                JIMAG = JIMAG + 1
00322                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00323      $                           RWORK( JIMAG ) )
00324   150       CONTINUE
00325   160    CONTINUE
00326 *
00327 *        Unscale.
00328 *
00329          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00330          CALL DLASRT( 'D', N, D, INFO )
00331          CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00332 *
00333          RETURN
00334       END IF
00335 *
00336 *     Book-keeping and setting up some constants.
00337 *
00338       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
00339 *
00340       SMLSZP = SMLSIZ + 1
00341 *
00342       U = 1
00343       VT = 1 + SMLSIZ*N
00344       DIFL = VT + SMLSZP*N
00345       DIFR = DIFL + NLVL*N
00346       Z = DIFR + NLVL*N*2
00347       C = Z + NLVL*N
00348       S = C + N
00349       POLES = S + N
00350       GIVNUM = POLES + 2*NLVL*N
00351       NRWORK = GIVNUM + 2*NLVL*N
00352       BX = 1
00353 *
00354       IRWRB = NRWORK
00355       IRWIB = IRWRB + SMLSIZ*NRHS
00356       IRWB = IRWIB + SMLSIZ*NRHS
00357 *
00358       SIZEI = 1 + N
00359       K = SIZEI + N
00360       GIVPTR = K + N
00361       PERM = GIVPTR + N
00362       GIVCOL = PERM + NLVL*N
00363       IWK = GIVCOL + NLVL*N*2
00364 *
00365       ST = 1
00366       SQRE = 0
00367       ICMPQ1 = 1
00368       ICMPQ2 = 0
00369       NSUB = 0
00370 *
00371       DO 170 I = 1, N
00372          IF( ABS( D( I ) ).LT.EPS ) THEN
00373             D( I ) = SIGN( EPS, D( I ) )
00374          END IF
00375   170 CONTINUE
00376 *
00377       DO 240 I = 1, NM1
00378          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
00379             NSUB = NSUB + 1
00380             IWORK( NSUB ) = ST
00381 *
00382 *           Subproblem found. First determine its size and then
00383 *           apply divide and conquer on it.
00384 *
00385             IF( I.LT.NM1 ) THEN
00386 *
00387 *              A subproblem with E(I) small for I < NM1.
00388 *
00389                NSIZE = I - ST + 1
00390                IWORK( SIZEI+NSUB-1 ) = NSIZE
00391             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
00392 *
00393 *              A subproblem with E(NM1) not too small but I = NM1.
00394 *
00395                NSIZE = N - ST + 1
00396                IWORK( SIZEI+NSUB-1 ) = NSIZE
00397             ELSE
00398 *
00399 *              A subproblem with E(NM1) small. This implies an
00400 *              1-by-1 subproblem at D(N), which is not solved
00401 *              explicitly.
00402 *
00403                NSIZE = I - ST + 1
00404                IWORK( SIZEI+NSUB-1 ) = NSIZE
00405                NSUB = NSUB + 1
00406                IWORK( NSUB ) = N
00407                IWORK( SIZEI+NSUB-1 ) = 1
00408                CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
00409             END IF
00410             ST1 = ST - 1
00411             IF( NSIZE.EQ.1 ) THEN
00412 *
00413 *              This is a 1-by-1 subproblem and is not solved
00414 *              explicitly.
00415 *
00416                CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
00417             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00418 *
00419 *              This is a small subproblem and is solved by DLASDQ.
00420 *
00421                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00422      $                      RWORK( VT+ST1 ), N )
00423                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00424      $                      RWORK( U+ST1 ), N )
00425                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
00426      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
00427      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
00428      $                      INFO )
00429                IF( INFO.NE.0 ) THEN
00430                   RETURN
00431                END IF
00432 *
00433 *              In the real version, B is passed to DLASDQ and multiplied
00434 *              internally by Q'. Here B is complex and that product is
00435 *              computed below in two steps (real and imaginary parts).
00436 *
00437                J = IRWB - 1
00438                DO 190 JCOL = 1, NRHS
00439                   DO 180 JROW = ST, ST + NSIZE - 1
00440                      J = J + 1
00441                      RWORK( J ) = DBLE( B( JROW, JCOL ) )
00442   180             CONTINUE
00443   190          CONTINUE
00444                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00445      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00446      $                     ZERO, RWORK( IRWRB ), NSIZE )
00447                J = IRWB - 1
00448                DO 210 JCOL = 1, NRHS
00449                   DO 200 JROW = ST, ST + NSIZE - 1
00450                      J = J + 1
00451                      RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00452   200             CONTINUE
00453   210          CONTINUE
00454                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00455      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00456      $                     ZERO, RWORK( IRWIB ), NSIZE )
00457                JREAL = IRWRB - 1
00458                JIMAG = IRWIB - 1
00459                DO 230 JCOL = 1, NRHS
00460                   DO 220 JROW = ST, ST + NSIZE - 1
00461                      JREAL = JREAL + 1
00462                      JIMAG = JIMAG + 1
00463                      B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00464      $                                 RWORK( JIMAG ) )
00465   220             CONTINUE
00466   230          CONTINUE
00467 *
00468                CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
00469      $                      WORK( BX+ST1 ), N )
00470             ELSE
00471 *
00472 *              A large problem. Solve it using divide and conquer.
00473 *
00474                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
00475      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
00476      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
00477      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
00478      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
00479      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
00480      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
00481      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
00482      $                      IWORK( IWK ), INFO )
00483                IF( INFO.NE.0 ) THEN
00484                   RETURN
00485                END IF
00486                BXST = BX + ST1
00487                CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
00488      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
00489      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
00490      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00491      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00492      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00493      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00494      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
00495      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
00496                IF( INFO.NE.0 ) THEN
00497                   RETURN
00498                END IF
00499             END IF
00500             ST = I + 1
00501          END IF
00502   240 CONTINUE
00503 *
00504 *     Apply the singular values and treat the tiny ones as zero.
00505 *
00506       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00507 *
00508       DO 250 I = 1, N
00509 *
00510 *        Some of the elements in D can be negative because 1-by-1
00511 *        subproblems were not solved explicitly.
00512 *
00513          IF( ABS( D( I ) ).LE.TOL ) THEN
00514             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
00515          ELSE
00516             RANK = RANK + 1
00517             CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
00518      $                   WORK( BX+I-1 ), N, INFO )
00519          END IF
00520          D( I ) = ABS( D( I ) )
00521   250 CONTINUE
00522 *
00523 *     Now apply back the right singular vectors.
00524 *
00525       ICMPQ2 = 1
00526       DO 320 I = 1, NSUB
00527          ST = IWORK( I )
00528          ST1 = ST - 1
00529          NSIZE = IWORK( SIZEI+I-1 )
00530          BXST = BX + ST1
00531          IF( NSIZE.EQ.1 ) THEN
00532             CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
00533          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00534 *
00535 *           Since B and BX are complex, the following call to DGEMM
00536 *           is performed in two steps (real and imaginary parts).
00537 *
00538 *           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00539 *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
00540 *    $                  B( ST, 1 ), LDB )
00541 *
00542             J = BXST - N - 1
00543             JREAL = IRWB - 1
00544             DO 270 JCOL = 1, NRHS
00545                J = J + N
00546                DO 260 JROW = 1, NSIZE
00547                   JREAL = JREAL + 1
00548                   RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
00549   260          CONTINUE
00550   270       CONTINUE
00551             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00552      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00553      $                  RWORK( IRWRB ), NSIZE )
00554             J = BXST - N - 1
00555             JIMAG = IRWB - 1
00556             DO 290 JCOL = 1, NRHS
00557                J = J + N
00558                DO 280 JROW = 1, NSIZE
00559                   JIMAG = JIMAG + 1
00560                   RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
00561   280          CONTINUE
00562   290       CONTINUE
00563             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00564      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00565      $                  RWORK( IRWIB ), NSIZE )
00566             JREAL = IRWRB - 1
00567             JIMAG = IRWIB - 1
00568             DO 310 JCOL = 1, NRHS
00569                DO 300 JROW = ST, ST + NSIZE - 1
00570                   JREAL = JREAL + 1
00571                   JIMAG = JIMAG + 1
00572                   B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00573      $                              RWORK( JIMAG ) )
00574   300          CONTINUE
00575   310       CONTINUE
00576          ELSE
00577             CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
00578      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
00579      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
00580      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00581      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00582      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00583      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00584      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
00585      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
00586             IF( INFO.NE.0 ) THEN
00587                RETURN
00588             END IF
00589          END IF
00590   320 CONTINUE
00591 *
00592 *     Unscale and sort the singular values.
00593 *
00594       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00595       CALL DLASRT( 'D', N, D, INFO )
00596       CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00597 *
00598       RETURN
00599 *
00600 *     End of ZLALSD
00601 *
00602       END