LAPACK: SRC/sgelss.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
00002      $                   WORK, LWORK, 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       REAL               A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SGELSS computes the minimum norm solution to a real linear least
00021 *  squares problem:
00022 *
00023 *  Minimize 2-norm(| b - A*x |).
00024 *
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 matrix
00031 *  X.
00032 *
00033 *  The effective rank of A is determined by treating as zero those
00034 *  singular values which are less than RCOND times the largest singular
00035 *  value.
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  M       (input) INTEGER
00041 *          The number of rows of the matrix A. M >= 0.
00042 *
00043 *  N       (input) INTEGER
00044 *          The number of columns of the matrix A. N >= 0.
00045 *
00046 *  NRHS    (input) INTEGER
00047 *          The number of right hand sides, i.e., the number of columns
00048 *          of the matrices B and X. NRHS >= 0.
00049 *
00050 *  A       (input/output) REAL array, dimension (LDA,N)
00051 *          On entry, the M-by-N matrix A.
00052 *          On exit, the first min(m,n) rows of A are overwritten with
00053 *          its right singular vectors, stored rowwise.
00054 *
00055 *  LDA     (input) INTEGER
00056 *          The leading dimension of the array A.  LDA >= max(1,M).
00057 *
00058 *  B       (input/output) REAL array, dimension (LDB,NRHS)
00059 *          On entry, the M-by-NRHS right hand side matrix B.
00060 *          On exit, B is overwritten by the N-by-NRHS solution
00061 *          matrix X.  If m >= n and RANK = n, the residual
00062 *          sum-of-squares for the solution in the i-th column is given
00063 *          by the sum of squares of elements n+1:m in that column.
00064 *
00065 *  LDB     (input) INTEGER
00066 *          The leading dimension of the array B. LDB >= max(1,max(M,N)).
00067 *
00068 *  S       (output) REAL array, dimension (min(M,N))
00069 *          The singular values of A in decreasing order.
00070 *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
00071 *
00072 *  RCOND   (input) REAL
00073 *          RCOND is used to determine the effective rank of A.
00074 *          Singular values S(i) <= RCOND*S(1) are treated as zero.
00075 *          If RCOND < 0, machine precision is used instead.
00076 *
00077 *  RANK    (output) INTEGER
00078 *          The effective rank of A, i.e., the number of singular values
00079 *          which are greater than RCOND*S(1).
00080 *
00081 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00082 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00083 *
00084 *  LWORK   (input) INTEGER
00085 *          The dimension of the array WORK. LWORK >= 1, and also:
00086 *          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
00087 *          For good performance, LWORK should generally be larger.
00088 *
00089 *          If LWORK = -1, then a workspace query is assumed; the routine
00090 *          only calculates the optimal size of the WORK array, returns
00091 *          this value as the first entry of the WORK array, and no error
00092 *          message related to LWORK is issued by XERBLA.
00093 *
00094 *  INFO    (output) INTEGER
00095 *          = 0:  successful exit
00096 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00097 *          > 0:  the algorithm for computing the SVD failed to converge;
00098 *                if INFO = i, i off-diagonal elements of an intermediate
00099 *                bidiagonal form did not converge to zero.
00100 *
00101 *  =====================================================================
00102 *
00103 *     .. Parameters ..
00104       REAL               ZERO, ONE
00105       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00106 *     ..
00107 *     .. Local Scalars ..
00108       LOGICAL            LQUERY
00109       INTEGER            BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
00110      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
00111      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
00112       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
00113 *     ..
00114 *     .. Local Arrays ..
00115       REAL               VDUM( 1 )
00116 *     ..
00117 *     .. External Subroutines ..
00118       EXTERNAL           SBDSQR, SCOPY, SGEBRD, SGELQF, SGEMM, SGEMV,
00119      $                   SGEQRF, SLABAD, SLACPY, SLASCL, SLASET, SORGBR,
00120      $                   SORMBR, SORMLQ, SORMQR, SRSCL, XERBLA
00121 *     ..
00122 *     .. External Functions ..
00123       INTEGER            ILAENV
00124       REAL               SLAMCH, SLANGE
00125       EXTERNAL           ILAENV, SLAMCH, SLANGE
00126 *     ..
00127 *     .. Intrinsic Functions ..
00128       INTRINSIC          MAX, MIN
00129 *     ..
00130 *     .. Executable Statements ..
00131 *
00132 *     Test the input arguments
00133 *
00134       INFO = 0
00135       MINMN = MIN( M, N )
00136       MAXMN = MAX( M, N )
00137       LQUERY = ( LWORK.EQ.-1 )
00138       IF( M.LT.0 ) THEN
00139          INFO = -1
00140       ELSE IF( N.LT.0 ) THEN
00141          INFO = -2
00142       ELSE IF( NRHS.LT.0 ) THEN
00143          INFO = -3
00144       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00145          INFO = -5
00146       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
00147          INFO = -7
00148       END IF
00149 *
00150 *     Compute workspace
00151 *      (Note: Comments in the code beginning "Workspace:" describe the
00152 *       minimal amount of workspace needed at that point in the code,
00153 *       as well as the preferred amount for good performance.
00154 *       NB refers to the optimal block size for the immediately
00155 *       following subroutine, as returned by ILAENV.)
00156 *
00157       IF( INFO.EQ.0 ) THEN
00158          MINWRK = 1
00159          MAXWRK = 1
00160          IF( MINMN.GT.0 ) THEN
00161             MM = M
00162             MNTHR = ILAENV( 6, 'SGELSS', ' ', M, N, NRHS, -1 )
00163             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
00164 *
00165 *              Path 1a - overdetermined, with many more rows than
00166 *                        columns
00167 *
00168                MM = N
00169                MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'SGEQRF', ' ', M,
00170      $                       N, -1, -1 ) )
00171                MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'SORMQR', 'LT',
00172      $                       M, NRHS, N, -1 ) )
00173             END IF
00174             IF( M.GE.N ) THEN
00175 *
00176 *              Path 1 - overdetermined or exactly determined
00177 *
00178 *              Compute workspace needed for SBDSQR
00179 *
00180                BDSPAC = MAX( 1, 5*N )
00181                MAXWRK = MAX( MAXWRK, 3*N + ( MM + N )*ILAENV( 1,
00182      $                       'SGEBRD', ' ', MM, N, -1, -1 ) )
00183                MAXWRK = MAX( MAXWRK, 3*N + NRHS*ILAENV( 1, 'SORMBR',
00184      $                       'QLT', MM, NRHS, N, -1 ) )
00185                MAXWRK = MAX( MAXWRK, 3*N + ( N - 1 )*ILAENV( 1,
00186      $                       'SORGBR', 'P', N, N, N, -1 ) )
00187                MAXWRK = MAX( MAXWRK, BDSPAC )
00188                MAXWRK = MAX( MAXWRK, N*NRHS )
00189                MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC )
00190                MAXWRK = MAX( MINWRK, MAXWRK )
00191             END IF
00192             IF( N.GT.M ) THEN
00193 *
00194 *              Compute workspace needed for SBDSQR
00195 *
00196                BDSPAC = MAX( 1, 5*M )
00197                MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
00198                IF( N.GE.MNTHR ) THEN
00199 *
00200 *                 Path 2a - underdetermined, with many more columns
00201 *                 than rows
00202 *
00203                   MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1,
00204      $                                  -1 )
00205                   MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
00206      $                          'SGEBRD', ' ', M, M, -1, -1 ) )
00207                   MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
00208      $                          'SORMBR', 'QLT', M, NRHS, M, -1 ) )
00209                   MAXWRK = MAX( MAXWRK, M*M + 4*M +
00210      $                          ( M - 1 )*ILAENV( 1, 'SORGBR', 'P', M,
00211      $                          M, M, -1 ) )
00212                   MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC )
00213                   IF( NRHS.GT.1 ) THEN
00214                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
00215                   ELSE
00216                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
00217                   END IF
00218                   MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'SORMLQ',
00219      $                          'LT', N, NRHS, M, -1 ) )
00220                ELSE
00221 *
00222 *                 Path 2 - underdetermined
00223 *
00224                   MAXWRK = 3*M + ( N + M )*ILAENV( 1, 'SGEBRD', ' ', M,
00225      $                     N, -1, -1 )
00226                   MAXWRK = MAX( MAXWRK, 3*M + NRHS*ILAENV( 1, 'SORMBR',
00227      $                          'QLT', M, NRHS, M, -1 ) )
00228                   MAXWRK = MAX( MAXWRK, 3*M + M*ILAENV( 1, 'SORGBR',
00229      $                          'P', M, N, M, -1 ) )
00230                   MAXWRK = MAX( MAXWRK, BDSPAC )
00231                   MAXWRK = MAX( MAXWRK, N*NRHS )
00232                END IF
00233             END IF
00234             MAXWRK = MAX( MINWRK, MAXWRK )
00235          END IF
00236          WORK( 1 ) = MAXWRK
00237 *
00238          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
00239      $      INFO = -12
00240       END IF
00241 *
00242       IF( INFO.NE.0 ) THEN
00243          CALL XERBLA( 'SGELSS', -INFO )
00244          RETURN
00245       ELSE IF( LQUERY ) THEN
00246          RETURN
00247       END IF
00248 *
00249 *     Quick return if possible
00250 *
00251       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00252          RANK = 0
00253          RETURN
00254       END IF
00255 *
00256 *     Get machine parameters
00257 *
00258       EPS = SLAMCH( 'P' )
00259       SFMIN = SLAMCH( 'S' )
00260       SMLNUM = SFMIN / EPS
00261       BIGNUM = ONE / SMLNUM
00262       CALL SLABAD( SMLNUM, BIGNUM )
00263 *
00264 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00265 *
00266       ANRM = SLANGE( 'M', M, N, A, LDA, WORK )
00267       IASCL = 0
00268       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00269 *
00270 *        Scale matrix norm up to SMLNUM
00271 *
00272          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00273          IASCL = 1
00274       ELSE IF( ANRM.GT.BIGNUM ) THEN
00275 *
00276 *        Scale matrix norm down to BIGNUM
00277 *
00278          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00279          IASCL = 2
00280       ELSE IF( ANRM.EQ.ZERO ) THEN
00281 *
00282 *        Matrix all zero. Return zero solution.
00283 *
00284          CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00285          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
00286          RANK = 0
00287          GO TO 70
00288       END IF
00289 *
00290 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00291 *
00292       BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK )
00293       IBSCL = 0
00294       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00295 *
00296 *        Scale matrix norm up to SMLNUM
00297 *
00298          CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00299          IBSCL = 1
00300       ELSE IF( BNRM.GT.BIGNUM ) THEN
00301 *
00302 *        Scale matrix norm down to BIGNUM
00303 *
00304          CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00305          IBSCL = 2
00306       END IF
00307 *
00308 *     Overdetermined case
00309 *
00310       IF( M.GE.N ) THEN
00311 *
00312 *        Path 1 - overdetermined or exactly determined
00313 *
00314          MM = M
00315          IF( M.GE.MNTHR ) THEN
00316 *
00317 *           Path 1a - overdetermined, with many more rows than columns
00318 *
00319             MM = N
00320             ITAU = 1
00321             IWORK = ITAU + N
00322 *
00323 *           Compute A=Q*R
00324 *           (Workspace: need 2*N, prefer N+N*NB)
00325 *
00326             CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00327      $                   LWORK-IWORK+1, INFO )
00328 *
00329 *           Multiply B by transpose(Q)
00330 *           (Workspace: need N+NRHS, prefer N+NRHS*NB)
00331 *
00332             CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
00333      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00334 *
00335 *           Zero out below R
00336 *
00337             IF( N.GT.1 )
00338      $         CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
00339          END IF
00340 *
00341          IE = 1
00342          ITAUQ = IE + N
00343          ITAUP = ITAUQ + N
00344          IWORK = ITAUP + N
00345 *
00346 *        Bidiagonalize R in A
00347 *        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
00348 *
00349          CALL SGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00350      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00351      $                INFO )
00352 *
00353 *        Multiply B by transpose of left bidiagonalizing vectors of R
00354 *        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
00355 *
00356          CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
00357      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00358 *
00359 *        Generate right bidiagonalizing vectors of R in A
00360 *        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
00361 *
00362          CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
00363      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00364          IWORK = IE + N
00365 *
00366 *        Perform bidiagonal QR iteration
00367 *          multiply B by transpose of left singular vectors
00368 *          compute right singular vectors in A
00369 *        (Workspace: need BDSPAC)
00370 *
00371          CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM,
00372      $                1, B, LDB, WORK( IWORK ), INFO )
00373          IF( INFO.NE.0 )
00374      $      GO TO 70
00375 *
00376 *        Multiply B by reciprocals of singular values
00377 *
00378          THR = MAX( RCOND*S( 1 ), SFMIN )
00379          IF( RCOND.LT.ZERO )
00380      $      THR = MAX( EPS*S( 1 ), SFMIN )
00381          RANK = 0
00382          DO 10 I = 1, N
00383             IF( S( I ).GT.THR ) THEN
00384                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00385                RANK = RANK + 1
00386             ELSE
00387                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00388             END IF
00389    10    CONTINUE
00390 *
00391 *        Multiply B by right singular vectors
00392 *        (Workspace: need N, prefer N*NRHS)
00393 *
00394          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00395             CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
00396      $                  WORK, LDB )
00397             CALL SLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
00398          ELSE IF( NRHS.GT.1 ) THEN
00399             CHUNK = LWORK / N
00400             DO 20 I = 1, NRHS, CHUNK
00401                BL = MIN( NRHS-I+1, CHUNK )
00402                CALL SGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ),
00403      $                     LDB, ZERO, WORK, N )
00404                CALL SLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
00405    20       CONTINUE
00406          ELSE
00407             CALL SGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
00408             CALL SCOPY( N, WORK, 1, B, 1 )
00409          END IF
00410 *
00411       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
00412      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
00413 *
00414 *        Path 2a - underdetermined, with many more columns than rows
00415 *        and sufficient workspace for an efficient algorithm
00416 *
00417          LDWORK = M
00418          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
00419      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
00420          ITAU = 1
00421          IWORK = M + 1
00422 *
00423 *        Compute A=L*Q
00424 *        (Workspace: need 2*M, prefer M+M*NB)
00425 *
00426          CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00427      $                LWORK-IWORK+1, INFO )
00428          IL = IWORK
00429 *
00430 *        Copy L to WORK(IL), zeroing out above it
00431 *
00432          CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
00433          CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
00434      $                LDWORK )
00435          IE = IL + LDWORK*M
00436          ITAUQ = IE + M
00437          ITAUP = ITAUQ + M
00438          IWORK = ITAUP + M
00439 *
00440 *        Bidiagonalize L in WORK(IL)
00441 *        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
00442 *
00443          CALL SGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
00444      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
00445      $                LWORK-IWORK+1, INFO )
00446 *
00447 *        Multiply B by transpose of left bidiagonalizing vectors of L
00448 *        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
00449 *
00450          CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
00451      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
00452      $                LWORK-IWORK+1, INFO )
00453 *
00454 *        Generate right bidiagonalizing vectors of R in WORK(IL)
00455 *        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
00456 *
00457          CALL SORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
00458      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00459          IWORK = IE + M
00460 *
00461 *        Perform bidiagonal QR iteration,
00462 *           computing right singular vectors of L in WORK(IL) and
00463 *           multiplying B by transpose of left singular vectors
00464 *        (Workspace: need M*M+M+BDSPAC)
00465 *
00466          CALL SBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ),
00467      $                LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO )
00468          IF( INFO.NE.0 )
00469      $      GO TO 70
00470 *
00471 *        Multiply B by reciprocals of singular values
00472 *
00473          THR = MAX( RCOND*S( 1 ), SFMIN )
00474          IF( RCOND.LT.ZERO )
00475      $      THR = MAX( EPS*S( 1 ), SFMIN )
00476          RANK = 0
00477          DO 30 I = 1, M
00478             IF( S( I ).GT.THR ) THEN
00479                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00480                RANK = RANK + 1
00481             ELSE
00482                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00483             END IF
00484    30    CONTINUE
00485          IWORK = IE
00486 *
00487 *        Multiply B by right singular vectors of L in WORK(IL)
00488 *        (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
00489 *
00490          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
00491             CALL SGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK,
00492      $                  B, LDB, ZERO, WORK( IWORK ), LDB )
00493             CALL SLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
00494          ELSE IF( NRHS.GT.1 ) THEN
00495             CHUNK = ( LWORK-IWORK+1 ) / M
00496             DO 40 I = 1, NRHS, CHUNK
00497                BL = MIN( NRHS-I+1, CHUNK )
00498                CALL SGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK,
00499      $                     B( 1, I ), LDB, ZERO, WORK( IWORK ), M )
00500                CALL SLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
00501      $                      LDB )
00502    40       CONTINUE
00503          ELSE
00504             CALL SGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ),
00505      $                  1, ZERO, WORK( IWORK ), 1 )
00506             CALL SCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
00507          END IF
00508 *
00509 *        Zero out below first M rows of B
00510 *
00511          CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
00512          IWORK = ITAU + M
00513 *
00514 *        Multiply transpose(Q) by B
00515 *        (Workspace: need M+NRHS, prefer M+NRHS*NB)
00516 *
00517          CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
00518      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00519 *
00520       ELSE
00521 *
00522 *        Path 2 - remaining underdetermined cases
00523 *
00524          IE = 1
00525          ITAUQ = IE + M
00526          ITAUP = ITAUQ + M
00527          IWORK = ITAUP + M
00528 *
00529 *        Bidiagonalize A
00530 *        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
00531 *
00532          CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00533      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00534      $                INFO )
00535 *
00536 *        Multiply B by transpose of left bidiagonalizing vectors
00537 *        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
00538 *
00539          CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
00540      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00541 *
00542 *        Generate right bidiagonalizing vectors in A
00543 *        (Workspace: need 4*M, prefer 3*M+M*NB)
00544 *
00545          CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
00546      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00547          IWORK = IE + M
00548 *
00549 *        Perform bidiagonal QR iteration,
00550 *           computing right singular vectors of A in A and
00551 *           multiplying B by transpose of left singular vectors
00552 *        (Workspace: need BDSPAC)
00553 *
00554          CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM,
00555      $                1, B, LDB, WORK( IWORK ), INFO )
00556          IF( INFO.NE.0 )
00557      $      GO TO 70
00558 *
00559 *        Multiply B by reciprocals of singular values
00560 *
00561          THR = MAX( RCOND*S( 1 ), SFMIN )
00562          IF( RCOND.LT.ZERO )
00563      $      THR = MAX( EPS*S( 1 ), SFMIN )
00564          RANK = 0
00565          DO 50 I = 1, M
00566             IF( S( I ).GT.THR ) THEN
00567                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00568                RANK = RANK + 1
00569             ELSE
00570                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00571             END IF
00572    50    CONTINUE
00573 *
00574 *        Multiply B by right singular vectors of A
00575 *        (Workspace: need N, prefer N*NRHS)
00576 *
00577          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00578             CALL SGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO,
00579      $                  WORK, LDB )
00580             CALL SLACPY( 'F', N, NRHS, WORK, LDB, B, LDB )
00581          ELSE IF( NRHS.GT.1 ) THEN
00582             CHUNK = LWORK / N
00583             DO 60 I = 1, NRHS, CHUNK
00584                BL = MIN( NRHS-I+1, CHUNK )
00585                CALL SGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ),
00586      $                     LDB, ZERO, WORK, N )
00587                CALL SLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
00588    60       CONTINUE
00589          ELSE
00590             CALL SGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
00591             CALL SCOPY( N, WORK, 1, B, 1 )
00592          END IF
00593       END IF
00594 *
00595 *     Undo scaling
00596 *
00597       IF( IASCL.EQ.1 ) THEN
00598          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00599          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
00600      $                INFO )
00601       ELSE IF( IASCL.EQ.2 ) THEN
00602          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00603          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
00604      $                INFO )
00605       END IF
00606       IF( IBSCL.EQ.1 ) THEN
00607          CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00608       ELSE IF( IBSCL.EQ.2 ) THEN
00609          CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00610       END IF
00611 *
00612    70 CONTINUE
00613       WORK( 1 ) = MAXWRK
00614       RETURN
00615 *
00616 *     End of SGELSS
00617 *
00618       END