LAPACK: SRC/cgelsd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
00002      $                   WORK, LWORK, RWORK, IWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.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 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
00011       REAL               RCOND
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IWORK( * )
00015       REAL               RWORK( * ), S( * )
00016       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  CGELSD computes the minimum-norm solution to a real linear least
00023 *  squares problem:
00024 *      minimize 2-norm(| b - A*x |)
00025 *  using the singular value decomposition (SVD) of A. A is an M-by-N
00026 *  matrix which may be rank-deficient.
00027 *
00028 *  Several right hand side vectors b and solution vectors x can be
00029 *  handled in a single call; they are stored as the columns of the
00030 *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00031 *  matrix X.
00032 *
00033 *  The problem is solved in three steps:
00034 *  (1) Reduce the coefficient matrix A to bidiagonal form with
00035 *      Householder tranformations, reducing the original problem
00036 *      into a "bidiagonal least squares problem" (BLS)
00037 *  (2) Solve the BLS using a divide and conquer approach.
00038 *  (3) Apply back all the Householder tranformations to solve
00039 *      the original least squares problem.
00040 *
00041 *  The effective rank of A is determined by treating as zero those
00042 *  singular values which are less than RCOND times the largest singular
00043 *  value.
00044 *
00045 *  The divide and conquer algorithm makes very mild assumptions about
00046 *  floating point arithmetic. It will work on machines with a guard
00047 *  digit in add/subtract, or on those binary machines without guard
00048 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00049 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
00050 *  without guard digits, but we know of none.
00051 *
00052 *  Arguments
00053 *  =========
00054 *
00055 *  M       (input) INTEGER
00056 *          The number of rows of the matrix A. M >= 0.
00057 *
00058 *  N       (input) INTEGER
00059 *          The number of columns of the matrix A. N >= 0.
00060 *
00061 *  NRHS    (input) INTEGER
00062 *          The number of right hand sides, i.e., the number of columns
00063 *          of the matrices B and X. NRHS >= 0.
00064 *
00065 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00066 *          On entry, the M-by-N matrix A.
00067 *          On exit, A has been destroyed.
00068 *
00069 *  LDA     (input) INTEGER
00070 *          The leading dimension of the array A. LDA >= max(1,M).
00071 *
00072 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00073 *          On entry, the M-by-NRHS right hand side matrix B.
00074 *          On exit, B is overwritten by the N-by-NRHS solution matrix X.
00075 *          If m >= n and RANK = n, the residual sum-of-squares for
00076 *          the solution in the i-th column is given by the sum of
00077 *          squares of the modulus of elements n+1:m in that column.
00078 *
00079 *  LDB     (input) INTEGER
00080 *          The leading dimension of the array B.  LDB >= max(1,M,N).
00081 *
00082 *  S       (output) REAL array, dimension (min(M,N))
00083 *          The singular values of A in decreasing order.
00084 *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
00085 *
00086 *  RCOND   (input) REAL
00087 *          RCOND is used to determine the effective rank of A.
00088 *          Singular values S(i) <= RCOND*S(1) are treated as zero.
00089 *          If RCOND < 0, machine precision is used instead.
00090 *
00091 *  RANK    (output) INTEGER
00092 *          The effective rank of A, i.e., the number of singular values
00093 *          which are greater than RCOND*S(1).
00094 *
00095 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00096 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00097 *
00098 *  LWORK   (input) INTEGER
00099 *          The dimension of the array WORK. LWORK must be at least 1.
00100 *          The exact minimum amount of workspace needed depends on M,
00101 *          N and NRHS. As long as LWORK is at least
00102 *              2 * N + N * NRHS
00103 *          if M is greater than or equal to N or
00104 *              2 * M + M * NRHS
00105 *          if M is less than N, the code will execute correctly.
00106 *          For good performance, LWORK should generally be larger.
00107 *
00108 *          If LWORK = -1, then a workspace query is assumed; the routine
00109 *          only calculates the optimal size of the array WORK and the
00110 *          minimum sizes of the arrays RWORK and IWORK, and returns
00111 *          these values as the first entries of the WORK, RWORK and
00112 *          IWORK arrays, and no error message related to LWORK is issued
00113 *          by XERBLA.
00114 *
00115 *  RWORK   (workspace) REAL array, dimension (MAX(1,LRWORK))
00116 *          LRWORK >=
00117 *             10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
00118 *             MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
00119 *          if M is greater than or equal to N or
00120 *             10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
00121 *             MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
00122 *          if M is less than N, the code will execute correctly.
00123 *          SMLSIZ is returned by ILAENV and is equal to the maximum
00124 *          size of the subproblems at the bottom of the computation
00125 *          tree (usually about 25), and
00126 *             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
00127 *          On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
00128 *
00129 *  IWORK   (workspace) INTEGER array, dimension (MAX(1,LIWORK))
00130 *          LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
00131 *          where MINMN = MIN( M,N ).
00132 *          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
00133 *
00134 *  INFO    (output) INTEGER
00135 *          = 0: successful exit
00136 *          < 0: if INFO = -i, the i-th argument had an illegal value.
00137 *          > 0:  the algorithm for computing the SVD failed to converge;
00138 *                if INFO = i, i off-diagonal elements of an intermediate
00139 *                bidiagonal form did not converge to zero.
00140 *
00141 *  Further Details
00142 *  ===============
00143 *
00144 *  Based on contributions by
00145 *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00146 *       California at Berkeley, USA
00147 *     Osni Marques, LBNL/NERSC, USA
00148 *
00149 *  =====================================================================
00150 *
00151 *     .. Parameters ..
00152       REAL               ZERO, ONE, TWO
00153       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
00154       COMPLEX            CZERO
00155       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
00156 *     ..
00157 *     .. Local Scalars ..
00158       LOGICAL            LQUERY
00159       INTEGER            IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
00160      $                   LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
00161      $                   MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
00162       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
00163 *     ..
00164 *     .. External Subroutines ..
00165       EXTERNAL           CGEBRD, CGELQF, CGEQRF, CLACPY,
00166      $                   CLALSD, CLASCL, CLASET, CUNMBR,
00167      $                   CUNMLQ, CUNMQR, SLABAD, SLASCL,
00168      $                   SLASET, XERBLA
00169 *     ..
00170 *     .. External Functions ..
00171       INTEGER            ILAENV
00172       REAL               CLANGE, SLAMCH
00173       EXTERNAL           CLANGE, SLAMCH, ILAENV
00174 *     ..
00175 *     .. Intrinsic Functions ..
00176       INTRINSIC          INT, LOG, MAX, MIN, REAL
00177 *     ..
00178 *     .. Executable Statements ..
00179 *
00180 *     Test the input arguments.
00181 *
00182       INFO = 0
00183       MINMN = MIN( M, N )
00184       MAXMN = MAX( M, N )
00185       LQUERY = ( LWORK.EQ.-1 )
00186       IF( M.LT.0 ) THEN
00187          INFO = -1
00188       ELSE IF( N.LT.0 ) THEN
00189          INFO = -2
00190       ELSE IF( NRHS.LT.0 ) THEN
00191          INFO = -3
00192       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00193          INFO = -5
00194       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
00195          INFO = -7
00196       END IF
00197 *
00198 *     Compute workspace.
00199 *     (Note: Comments in the code beginning "Workspace:" describe the
00200 *     minimal amount of workspace needed at that point in the code,
00201 *     as well as the preferred amount for good performance.
00202 *     NB refers to the optimal block size for the immediately
00203 *     following subroutine, as returned by ILAENV.)
00204 *
00205       IF( INFO.EQ.0 ) THEN
00206          MINWRK = 1
00207          MAXWRK = 1
00208          LIWORK = 1
00209          LRWORK = 1
00210          IF( MINMN.GT.0 ) THEN
00211             SMLSIZ = ILAENV( 9, 'CGELSD', ' ', 0, 0, 0, 0 )
00212             MNTHR = ILAENV( 6, 'CGELSD', ' ', M, N, NRHS, -1 )
00213             NLVL = MAX( INT( LOG( REAL( MINMN ) / REAL( SMLSIZ + 1 ) ) /
00214      $                  LOG( TWO ) ) + 1, 0 )
00215             LIWORK = 3*MINMN*NLVL + 11*MINMN
00216             MM = M
00217             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
00218 *
00219 *              Path 1a - overdetermined, with many more rows than
00220 *                        columns.
00221 *
00222                MM = N
00223                MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'CGEQRF', ' ', M, N,
00224      $                       -1, -1 ) )
00225                MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'CUNMQR', 'LC', M,
00226      $                       NRHS, N, -1 ) )
00227             END IF
00228             IF( M.GE.N ) THEN
00229 *
00230 *              Path 1 - overdetermined or exactly determined.
00231 *
00232                LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
00233      $                  MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
00234                MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
00235      $                       'CGEBRD', ' ', MM, N, -1, -1 ) )
00236                MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR',
00237      $                       'QLC', MM, NRHS, N, -1 ) )
00238                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
00239      $                       'CUNMBR', 'PLN', N, NRHS, N, -1 ) )
00240                MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
00241                MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
00242             END IF
00243             IF( N.GT.M ) THEN
00244                LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
00245      $                  MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
00246                IF( N.GE.MNTHR ) THEN
00247 *
00248 *                 Path 2a - underdetermined, with many more columns
00249 *                           than rows.
00250 *
00251                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
00252      $                     -1 )
00253                   MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
00254      $                          'CGEBRD', ' ', M, M, -1, -1 ) )
00255                   MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
00256      $                          'CUNMBR', 'QLC', M, NRHS, M, -1 ) )
00257                   MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
00258      $                          'CUNMLQ', 'LC', N, NRHS, M, -1 ) )
00259                   IF( NRHS.GT.1 ) THEN
00260                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
00261                   ELSE
00262                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
00263                   END IF
00264                   MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
00265 !     XXX: Ensure the Path 2a case below is triggered.  The workspace
00266 !     calculation should use queries for all routines eventually.
00267                   MAXWRK = MAX( MAXWRK,
00268      $                 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
00269                ELSE
00270 *
00271 *                 Path 2 - underdetermined.
00272 *
00273                   MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M,
00274      $                     N, -1, -1 )
00275                   MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR',
00276      $                          'QLC', M, NRHS, M, -1 ) )
00277                   MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNMBR',
00278      $                          'PLN', N, NRHS, M, -1 ) )
00279                   MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
00280                END IF
00281                MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
00282             END IF
00283          END IF
00284          MINWRK = MIN( MINWRK, MAXWRK )
00285          WORK( 1 ) = MAXWRK
00286          IWORK( 1 ) = LIWORK
00287          RWORK( 1 ) = LRWORK
00288 *
00289          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00290             INFO = -12
00291          END IF
00292       END IF
00293 *
00294       IF( INFO.NE.0 ) THEN
00295          CALL XERBLA( 'CGELSD', -INFO )
00296          RETURN
00297       ELSE IF( LQUERY ) THEN
00298          RETURN
00299       END IF
00300 *
00301 *     Quick return if possible.
00302 *
00303       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00304          RANK = 0
00305          RETURN
00306       END IF
00307 *
00308 *     Get machine parameters.
00309 *
00310       EPS = SLAMCH( 'P' )
00311       SFMIN = SLAMCH( 'S' )
00312       SMLNUM = SFMIN / EPS
00313       BIGNUM = ONE / SMLNUM
00314       CALL SLABAD( SMLNUM, BIGNUM )
00315 *
00316 *     Scale A if max entry outside range [SMLNUM,BIGNUM].
00317 *
00318       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
00319       IASCL = 0
00320       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00321 *
00322 *        Scale matrix norm up to SMLNUM
00323 *
00324          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00325          IASCL = 1
00326       ELSE IF( ANRM.GT.BIGNUM ) THEN
00327 *
00328 *        Scale matrix norm down to BIGNUM.
00329 *
00330          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00331          IASCL = 2
00332       ELSE IF( ANRM.EQ.ZERO ) THEN
00333 *
00334 *        Matrix all zero. Return zero solution.
00335 *
00336          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00337          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
00338          RANK = 0
00339          GO TO 10
00340       END IF
00341 *
00342 *     Scale B if max entry outside range [SMLNUM,BIGNUM].
00343 *
00344       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
00345       IBSCL = 0
00346       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00347 *
00348 *        Scale matrix norm up to SMLNUM.
00349 *
00350          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00351          IBSCL = 1
00352       ELSE IF( BNRM.GT.BIGNUM ) THEN
00353 *
00354 *        Scale matrix norm down to BIGNUM.
00355 *
00356          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00357          IBSCL = 2
00358       END IF
00359 *
00360 *     If M < N make sure B(M+1:N,:) = 0
00361 *
00362       IF( M.LT.N )
00363      $   CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
00364 *
00365 *     Overdetermined case.
00366 *
00367       IF( M.GE.N ) THEN
00368 *
00369 *        Path 1 - overdetermined or exactly determined.
00370 *
00371          MM = M
00372          IF( M.GE.MNTHR ) THEN
00373 *
00374 *           Path 1a - overdetermined, with many more rows than columns
00375 *
00376             MM = N
00377             ITAU = 1
00378             NWORK = ITAU + N
00379 *
00380 *           Compute A=Q*R.
00381 *           (RWorkspace: need N)
00382 *           (CWorkspace: need N, prefer N*NB)
00383 *
00384             CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00385      $                   LWORK-NWORK+1, INFO )
00386 *
00387 *           Multiply B by transpose(Q).
00388 *           (RWorkspace: need N)
00389 *           (CWorkspace: need NRHS, prefer NRHS*NB)
00390 *
00391             CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
00392      $                   LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00393 *
00394 *           Zero out below R.
00395 *
00396             IF( N.GT.1 ) THEN
00397                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
00398      $                      LDA )
00399             END IF
00400          END IF
00401 *
00402          ITAUQ = 1
00403          ITAUP = ITAUQ + N
00404          NWORK = ITAUP + N
00405          IE = 1
00406          NRWORK = IE + N
00407 *
00408 *        Bidiagonalize R in A.
00409 *        (RWorkspace: need N)
00410 *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
00411 *
00412          CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00413      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00414      $                INFO )
00415 *
00416 *        Multiply B by transpose of left bidiagonalizing vectors of R.
00417 *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
00418 *
00419          CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
00420      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00421 *
00422 *        Solve the bidiagonal least squares problem.
00423 *
00424          CALL CLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
00425      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
00426      $                IWORK, INFO )
00427          IF( INFO.NE.0 ) THEN
00428             GO TO 10
00429          END IF
00430 *
00431 *        Multiply B by right bidiagonalizing vectors of R.
00432 *
00433          CALL CUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
00434      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00435 *
00436       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
00437      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
00438 *
00439 *        Path 2a - underdetermined, with many more columns than rows
00440 *        and sufficient workspace for an efficient algorithm.
00441 *
00442          LDWORK = M
00443          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
00444      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
00445          ITAU = 1
00446          NWORK = M + 1
00447 *
00448 *        Compute A=L*Q.
00449 *        (CWorkspace: need 2*M, prefer M+M*NB)
00450 *
00451          CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00452      $                LWORK-NWORK+1, INFO )
00453          IL = NWORK
00454 *
00455 *        Copy L to WORK(IL), zeroing out above its diagonal.
00456 *
00457          CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
00458          CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
00459      $                LDWORK )
00460          ITAUQ = IL + LDWORK*M
00461          ITAUP = ITAUQ + M
00462          NWORK = ITAUP + M
00463          IE = 1
00464          NRWORK = IE + M
00465 *
00466 *        Bidiagonalize L in WORK(IL).
00467 *        (RWorkspace: need M)
00468 *        (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
00469 *
00470          CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
00471      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
00472      $                LWORK-NWORK+1, INFO )
00473 *
00474 *        Multiply B by transpose of left bidiagonalizing vectors of L.
00475 *        (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
00476 *
00477          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
00478      $                WORK( ITAUQ ), B, LDB, WORK( NWORK ),
00479      $                LWORK-NWORK+1, INFO )
00480 *
00481 *        Solve the bidiagonal least squares problem.
00482 *
00483          CALL CLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
00484      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
00485      $                IWORK, INFO )
00486          IF( INFO.NE.0 ) THEN
00487             GO TO 10
00488          END IF
00489 *
00490 *        Multiply B by right bidiagonalizing vectors of L.
00491 *
00492          CALL CUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
00493      $                WORK( ITAUP ), B, LDB, WORK( NWORK ),
00494      $                LWORK-NWORK+1, INFO )
00495 *
00496 *        Zero out below first M rows of B.
00497 *
00498          CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
00499          NWORK = ITAU + M
00500 *
00501 *        Multiply transpose(Q) by B.
00502 *        (CWorkspace: need NRHS, prefer NRHS*NB)
00503 *
00504          CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
00505      $                LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00506 *
00507       ELSE
00508 *
00509 *        Path 2 - remaining underdetermined cases.
00510 *
00511          ITAUQ = 1
00512          ITAUP = ITAUQ + M
00513          NWORK = ITAUP + M
00514          IE = 1
00515          NRWORK = IE + M
00516 *
00517 *        Bidiagonalize A.
00518 *        (RWorkspace: need M)
00519 *        (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
00520 *
00521          CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00522      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00523      $                INFO )
00524 *
00525 *        Multiply B by transpose of left bidiagonalizing vectors.
00526 *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
00527 *
00528          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
00529      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00530 *
00531 *        Solve the bidiagonal least squares problem.
00532 *
00533          CALL CLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
00534      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
00535      $                IWORK, INFO )
00536          IF( INFO.NE.0 ) THEN
00537             GO TO 10
00538          END IF
00539 *
00540 *        Multiply B by right bidiagonalizing vectors of A.
00541 *
00542          CALL CUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
00543      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
00544 *
00545       END IF
00546 *
00547 *     Undo scaling.
00548 *
00549       IF( IASCL.EQ.1 ) THEN
00550          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00551          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
00552      $                INFO )
00553       ELSE IF( IASCL.EQ.2 ) THEN
00554          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00555          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
00556      $                INFO )
00557       END IF
00558       IF( IBSCL.EQ.1 ) THEN
00559          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00560       ELSE IF( IBSCL.EQ.2 ) THEN
00561          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00562       END IF
00563 *
00564    10 CONTINUE
00565       WORK( 1 ) = MAXWRK
00566       IWORK( 1 ) = LIWORK
00567       RWORK( 1 ) = LRWORK
00568       RETURN
00569 *
00570 *     End of CGELSD
00571 *
00572       END