zlalsd.f

Go to the documentation of this file.

*> \brief \b ZLALSD uses the singular value decomposition of A to solve the least squares problem.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download ZLALSD + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
*                          RANK, WORK, RWORK, IWORK, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
*       DOUBLE PRECISION   RCOND
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
*       COMPLEX*16         B( LDB, * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZLALSD uses the singular value decomposition of A to solve the least
*> squares problem of finding X to minimize the Euclidean norm of each
*> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
*> are N-by-NRHS. The solution X overwrites B.
*>
*> The singular values of A smaller than RCOND times the largest
*> singular value are treated as zero in solving the least squares
*> problem; in this case a minimum norm solution is returned.
*> The actual singular values are returned in D in ascending order.
*>
*> This code makes very mild assumptions about floating point
*> arithmetic. It will work on machines with a guard digit in
*> add/subtract, or on those binary machines without guard digits
*> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
*> It could conceivably fail on hexadecimal or decimal machines
*> without guard digits, but we know of none.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>         = 'U': D and E define an upper bidiagonal matrix.
*>         = 'L': D and E define a  lower bidiagonal matrix.
*> \endverbatim
*>
*> \param[in] SMLSIZ
*> \verbatim
*>          SMLSIZ is INTEGER
*>         The maximum size of the subproblems at the bottom of the
*>         computation tree.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>         The dimension of the  bidiagonal matrix.  N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>         The number of columns of B. NRHS must be at least 1.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>         On entry D contains the main diagonal of the bidiagonal
*>         matrix. On exit, if INFO = 0, D contains its singular values.
*> \endverbatim
*>
*> \param[in,out] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (N-1)
*>         Contains the super-diagonal entries of the bidiagonal matrix.
*>         On exit, E has been destroyed.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
*>         On input, B contains the right hand sides of the least
*>         squares problem. On output, B contains the solution X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>         The leading dimension of B in the calling subprogram.
*>         LDB must be at least max(1,N).
*> \endverbatim
*>
*> \param[in] RCOND
*> \verbatim
*>          RCOND is DOUBLE PRECISION
*>         The singular values of A less than or equal to RCOND times
*>         the largest singular value are treated as zero in solving
*>         the least squares problem. If RCOND is negative,
*>         machine precision is used instead.
*>         For example, if diag(S)*X=B were the least squares problem,
*>         where diag(S) is a diagonal matrix of singular values, the
*>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
*>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
*>         RCOND*max(S).
*> \endverbatim
*>
*> \param[out] RANK
*> \verbatim
*>          RANK is INTEGER
*>         The number of singular values of A greater than RCOND times
*>         the largest singular value.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension at least
*>         (N * NRHS).
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension at least
*>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
*>         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
*>         where
*>         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension at least
*>         (3*N*NLVL + 11*N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>         = 0:  successful exit.
*>         < 0:  if INFO = -i, the i-th argument had an illegal value.
*>         > 0:  The algorithm failed to compute a singular value while
*>               working on the submatrix lying in rows and columns
*>               INFO/(N+1) through MOD(INFO,N+1).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup complex16OTHERcomputational
*
*> \par Contributors:
*  ==================
*>
*>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
*>       California at Berkeley, USA \n
*>     Osni Marques, LBNL/NERSC, USA \n
*
*  =====================================================================
      SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
     $                   rank, work, rwork, iwork, info )
*
*  -- LAPACK computational routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
      COMPLEX*16         B( ldb, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter                ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
      COMPLEX*16         CZERO
      parameter                ( czero = ( 0.0d0, 0.0d0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
     $                   givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
     $                   irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
     $                   jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
     $                   perm, poles, s, sizei, smlszp, sqre, st, st1,
     $                   u, vt, z
      DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           idamax, dlamch, dlanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemm, dlartg, dlascl, dlasda, dlasdq, dlaset,
     $                   dlasrt, xerbla, zcopy, zdrot, zlacpy, zlalsa,
     $                   zlascl, zlaset
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dcmplx, dimag, int, log, sign
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( nrhs.LT.1 ) THEN
         info = -4
      ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
         info = -8
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZLALSD', -info )
         RETURN
      END IF
*
      eps = dlamch( 'Epsilon' )
*
*     Set up the tolerance.
*
      IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
         rcnd = eps
      ELSE
         rcnd = rcond
      END IF
*
      rank = 0
*
*     Quick return if possible.
*
      IF( n.EQ.0 ) THEN
         RETURN
      ELSE IF( n.EQ.1 ) THEN
         IF( d( 1 ).EQ.zero ) THEN
            CALL zlaset( 'A', 1, nrhs, czero, czero, b, ldb )
         ELSE
            rank = 1
            CALL zlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
            d( 1 ) = abs( d( 1 ) )
         END IF
         RETURN
      END IF
*
*     Rotate the matrix if it is lower bidiagonal.
*
      IF( uplo.EQ.'L' ) THEN
         DO 10 i = 1, n - 1
            CALL dlartg( d( i ), e( i ), cs, sn, r )
            d( i ) = r
            e( i ) = sn*d( i+1 )
            d( i+1 ) = cs*d( i+1 )
            IF( nrhs.EQ.1 ) THEN
               CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
            ELSE
               rwork( i*2-1 ) = cs
               rwork( i*2 ) = sn
            END IF
   10    CONTINUE
         IF( nrhs.GT.1 ) THEN
            DO 30 i = 1, nrhs
               DO 20 j = 1, n - 1
                  cs = rwork( j*2-1 )
                  sn = rwork( j*2 )
                  CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
   20          CONTINUE
   30       CONTINUE
         END IF
      END IF
*
*     Scale.
*
      nm1 = n - 1
      orgnrm = dlanst( 'M', n, d, e )
      IF( orgnrm.EQ.zero ) THEN
         CALL zlaset( 'A', n, nrhs, czero, czero, b, ldb )
         RETURN
      END IF
*
      CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
      CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
*
*     If N is smaller than the minimum divide size SMLSIZ, then solve
*     the problem with another solver.
*
      IF( n.LE.smlsiz ) THEN
         irwu = 1
         irwvt = irwu + n*n
         irwwrk = irwvt + n*n
         irwrb = irwwrk
         irwib = irwrb + n*nrhs
         irwb = irwib + n*nrhs
         CALL dlaset( 'A', n, n, zero, one, rwork( irwu ), n )
         CALL dlaset( 'A', n, n, zero, one, rwork( irwvt ), n )
         CALL dlasdq( 'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
     $                rwork( irwu ), n, rwork( irwwrk ), 1,
     $                rwork( irwwrk ), info )
         IF( info.NE.0 ) THEN
            RETURN
         END IF
*
*        In the real version, B is passed to DLASDQ and multiplied
*        internally by Q**H. Here B is complex and that product is
*        computed below in two steps (real and imaginary parts).
*
         j = irwb - 1
         DO 50 jcol = 1, nrhs
            DO 40 jrow = 1, n
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
   40       CONTINUE
   50    CONTINUE
         CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
     $               rwork( irwb ), n, zero, rwork( irwrb ), n )
         j = irwb - 1
         DO 70 jcol = 1, nrhs
            DO 60 jrow = 1, n
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
   60       CONTINUE
   70    CONTINUE
         CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
     $               rwork( irwb ), n, zero, rwork( irwib ), n )
         jreal = irwrb - 1
         jimag = irwib - 1
         DO 90 jcol = 1, nrhs
            DO 80 jrow = 1, n
               jreal = jreal + 1
               jimag = jimag + 1
               b( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                           rwork( jimag ) )
   80       CONTINUE
   90    CONTINUE
*
         tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
         DO 100 i = 1, n
            IF( d( i ).LE.tol ) THEN
               CALL zlaset( 'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
            ELSE
               CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
     $                      ldb, info )
               rank = rank + 1
            END IF
  100    CONTINUE
*
*        Since B is complex, the following call to DGEMM is performed
*        in two steps (real and imaginary parts). That is for V * B
*        (in the real version of the code V**H is stored in WORK).
*
*        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
*    $               WORK( NWORK ), N )
*
         j = irwb - 1
         DO 120 jcol = 1, nrhs
            DO 110 jrow = 1, n
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
  110       CONTINUE
  120    CONTINUE
         CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
     $               rwork( irwb ), n, zero, rwork( irwrb ), n )
         j = irwb - 1
         DO 140 jcol = 1, nrhs
            DO 130 jrow = 1, n
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
  130       CONTINUE
  140    CONTINUE
         CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
     $               rwork( irwb ), n, zero, rwork( irwib ), n )
         jreal = irwrb - 1
         jimag = irwib - 1
         DO 160 jcol = 1, nrhs
            DO 150 jrow = 1, n
               jreal = jreal + 1
               jimag = jimag + 1
               b( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                           rwork( jimag ) )
  150       CONTINUE
  160    CONTINUE
*
*        Unscale.
*
         CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
         CALL dlasrt( 'D', n, d, info )
         CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
*
         RETURN
      END IF
*
*     Book-keeping and setting up some constants.
*
      nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
*
      smlszp = smlsiz + 1
*
      u = 1
      vt = 1 + smlsiz*n
      difl = vt + smlszp*n
      difr = difl + nlvl*n
      z = difr + nlvl*n*2
      c = z + nlvl*n
      s = c + n
      poles = s + n
      givnum = poles + 2*nlvl*n
      nrwork = givnum + 2*nlvl*n
      bx = 1
*
      irwrb = nrwork
      irwib = irwrb + smlsiz*nrhs
      irwb = irwib + smlsiz*nrhs
*
      sizei = 1 + n
      k = sizei + n
      givptr = k + n
      perm = givptr + n
      givcol = perm + nlvl*n
      iwk = givcol + nlvl*n*2
*
      st = 1
      sqre = 0
      icmpq1 = 1
      icmpq2 = 0
      nsub = 0
*
      DO 170 i = 1, n
         IF( abs( d( i ) ).LT.eps ) THEN
            d( i ) = sign( eps, d( i ) )
         END IF
  170 CONTINUE
*
      DO 240 i = 1, nm1
         IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
            nsub = nsub + 1
            iwork( nsub ) = st
*
*           Subproblem found. First determine its size and then
*           apply divide and conquer on it.
*
            IF( i.LT.nm1 ) THEN
*
*              A subproblem with E(I) small for I < NM1.
*
               nsize = i - st + 1
               iwork( sizei+nsub-1 ) = nsize
            ELSE IF( abs( e( i ) ).GE.eps ) THEN
*
*              A subproblem with E(NM1) not too small but I = NM1.
*
               nsize = n - st + 1
               iwork( sizei+nsub-1 ) = nsize
            ELSE
*
*              A subproblem with E(NM1) small. This implies an
*              1-by-1 subproblem at D(N), which is not solved
*              explicitly.
*
               nsize = i - st + 1
               iwork( sizei+nsub-1 ) = nsize
               nsub = nsub + 1
               iwork( nsub ) = n
               iwork( sizei+nsub-1 ) = 1
               CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
            END IF
            st1 = st - 1
            IF( nsize.EQ.1 ) THEN
*
*              This is a 1-by-1 subproblem and is not solved
*              explicitly.
*
               CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
            ELSE IF( nsize.LE.smlsiz ) THEN
*
*              This is a small subproblem and is solved by DLASDQ.
*
               CALL dlaset( 'A', nsize, nsize, zero, one,
     $                      rwork( vt+st1 ), n )
               CALL dlaset( 'A', nsize, nsize, zero, one,
     $                      rwork( u+st1 ), n )
               CALL dlasdq( 'U', 0, nsize, nsize, nsize, 0, d( st ),
     $                      e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
     $                      n, rwork( nrwork ), 1, rwork( nrwork ),
     $                      info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
*
*              In the real version, B is passed to DLASDQ and multiplied
*              internally by Q**H. Here B is complex and that product is
*              computed below in two steps (real and imaginary parts).
*
               j = irwb - 1
               DO 190 jcol = 1, nrhs
                  DO 180 jrow = st, st + nsize - 1
                     j = j + 1
                     rwork( j ) = dble( b( jrow, jcol ) )
  180             CONTINUE
  190          CONTINUE
               CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
     $                     rwork( u+st1 ), n, rwork( irwb ), nsize,
     $                     zero, rwork( irwrb ), nsize )
               j = irwb - 1
               DO 210 jcol = 1, nrhs
                  DO 200 jrow = st, st + nsize - 1
                     j = j + 1
                     rwork( j ) = dimag( b( jrow, jcol ) )
  200             CONTINUE
  210          CONTINUE
               CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
     $                     rwork( u+st1 ), n, rwork( irwb ), nsize,
     $                     zero, rwork( irwib ), nsize )
               jreal = irwrb - 1
               jimag = irwib - 1
               DO 230 jcol = 1, nrhs
                  DO 220 jrow = st, st + nsize - 1
                     jreal = jreal + 1
                     jimag = jimag + 1
                     b( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                                 rwork( jimag ) )
  220             CONTINUE
  230          CONTINUE
*
               CALL zlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
     $                      work( bx+st1 ), n )
            ELSE
*
*              A large problem. Solve it using divide and conquer.
*
               CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
     $                      e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
     $                      iwork( k+st1 ), rwork( difl+st1 ),
     $                      rwork( difr+st1 ), rwork( z+st1 ),
     $                      rwork( poles+st1 ), iwork( givptr+st1 ),
     $                      iwork( givcol+st1 ), n, iwork( perm+st1 ),
     $                      rwork( givnum+st1 ), rwork( c+st1 ),
     $                      rwork( s+st1 ), rwork( nrwork ),
     $                      iwork( iwk ), info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
               bxst = bx + st1
               CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
     $                      ldb, work( bxst ), n, rwork( u+st1 ), n,
     $                      rwork( vt+st1 ), iwork( k+st1 ),
     $                      rwork( difl+st1 ), rwork( difr+st1 ),
     $                      rwork( z+st1 ), rwork( poles+st1 ),
     $                      iwork( givptr+st1 ), iwork( givcol+st1 ), n,
     $                      iwork( perm+st1 ), rwork( givnum+st1 ),
     $                      rwork( c+st1 ), rwork( s+st1 ),
     $                      rwork( nrwork ), iwork( iwk ), info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            st = i + 1
         END IF
  240 CONTINUE
*
*     Apply the singular values and treat the tiny ones as zero.
*
      tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
*
      DO 250 i = 1, n
*
*        Some of the elements in D can be negative because 1-by-1
*        subproblems were not solved explicitly.
*
         IF( abs( d( i ) ).LE.tol ) THEN
            CALL zlaset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
         ELSE
            rank = rank + 1
            CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
     $                   work( bx+i-1 ), n, info )
         END IF
         d( i ) = abs( d( i ) )
  250 CONTINUE
*
*     Now apply back the right singular vectors.
*
      icmpq2 = 1
      DO 320 i = 1, nsub
         st = iwork( i )
         st1 = st - 1
         nsize = iwork( sizei+i-1 )
         bxst = bx + st1
         IF( nsize.EQ.1 ) THEN
            CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
         ELSE IF( nsize.LE.smlsiz ) THEN
*
*           Since B and BX are complex, the following call to DGEMM
*           is performed in two steps (real and imaginary parts).
*
*           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
*    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
*    $                  B( ST, 1 ), LDB )
*
            j = bxst - n - 1
            jreal = irwb - 1
            DO 270 jcol = 1, nrhs
               j = j + n
               DO 260 jrow = 1, nsize
                  jreal = jreal + 1
                  rwork( jreal ) = dble( work( j+jrow ) )
  260          CONTINUE
  270       CONTINUE
            CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
     $                  rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
     $                  rwork( irwrb ), nsize )
            j = bxst - n - 1
            jimag = irwb - 1
            DO 290 jcol = 1, nrhs
               j = j + n
               DO 280 jrow = 1, nsize
                  jimag = jimag + 1
                  rwork( jimag ) = dimag( work( j+jrow ) )
  280          CONTINUE
  290       CONTINUE
            CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
     $                  rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
     $                  rwork( irwib ), nsize )
            jreal = irwrb - 1
            jimag = irwib - 1
            DO 310 jcol = 1, nrhs
               DO 300 jrow = st, st + nsize - 1
                  jreal = jreal + 1
                  jimag = jimag + 1
                  b( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                              rwork( jimag ) )
  300          CONTINUE
  310       CONTINUE
         ELSE
            CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
     $                   b( st, 1 ), ldb, rwork( u+st1 ), n,
     $                   rwork( vt+st1 ), iwork( k+st1 ),
     $                   rwork( difl+st1 ), rwork( difr+st1 ),
     $                   rwork( z+st1 ), rwork( poles+st1 ),
     $                   iwork( givptr+st1 ), iwork( givcol+st1 ), n,
     $                   iwork( perm+st1 ), rwork( givnum+st1 ),
     $                   rwork( c+st1 ), rwork( s+st1 ),
     $                   rwork( nrwork ), iwork( iwk ), info )
            IF( info.NE.0 ) THEN
               RETURN
            END IF
         END IF
  320 CONTINUE
*
*     Unscale and sort the singular values.
*
      CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
      CALL dlasrt( 'D', n, d, info )
      CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
*
      RETURN
*
*     End of ZLALSD
*
      END