lapack-netlib/SRC/sgelss.f

   1 *> \brief <b> SGELSS solves overdetermined or underdetermined systems for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGELSS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelss.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelss.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelss.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
  22 *                          WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
  26 *       REAL               RCOND
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SGELSS computes the minimum norm solution to a real linear least
  39 *> squares problem:
  40 *>
  41 *> Minimize 2-norm(| b - A*x |).
  42 *>
  43 *> using the singular value decomposition (SVD) of A. A is an M-by-N
  44 *> matrix which may be rank-deficient.
  45 *>
  46 *> Several right hand side vectors b and solution vectors x can be
  47 *> handled in a single call; they are stored as the columns of the
  48 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
  49 *> X.
  50 *>
  51 *> The effective rank of A is determined by treating as zero those
  52 *> singular values which are less than RCOND times the largest singular
  53 *> value.
  54 *> \endverbatim
  55 *
  56 *  Arguments:
  57 *  ==========
  58 *
  59 *> \param[in] M
  60 *> \verbatim
  61 *>          M is INTEGER
  62 *>          The number of rows of the matrix A. M >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] N
  66 *> \verbatim
  67 *>          N is INTEGER
  68 *>          The number of columns of the matrix A. N >= 0.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] NRHS
  72 *> \verbatim
  73 *>          NRHS is INTEGER
  74 *>          The number of right hand sides, i.e., the number of columns
  75 *>          of the matrices B and X. NRHS >= 0.
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] A
  79 *> \verbatim
  80 *>          A is REAL array, dimension (LDA,N)
  81 *>          On entry, the M-by-N matrix A.
  82 *>          On exit, the first min(m,n) rows of A are overwritten with
  83 *>          its right singular vectors, stored rowwise.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDA
  87 *> \verbatim
  88 *>          LDA is INTEGER
  89 *>          The leading dimension of the array A.  LDA >= max(1,M).
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] B
  93 *> \verbatim
  94 *>          B is REAL array, dimension (LDB,NRHS)
  95 *>          On entry, the M-by-NRHS right hand side matrix B.
  96 *>          On exit, B is overwritten by the N-by-NRHS solution
  97 *>          matrix X.  If m >= n and RANK = n, the residual
  98 *>          sum-of-squares for the solution in the i-th column is given
  99 *>          by the sum of squares of elements n+1:m in that column.
 100 *> \endverbatim
 101 *>
 102 *> \param[in] LDB
 103 *> \verbatim
 104 *>          LDB is INTEGER
 105 *>          The leading dimension of the array B. LDB >= max(1,max(M,N)).
 106 *> \endverbatim
 107 *>
 108 *> \param[out] S
 109 *> \verbatim
 110 *>          S is REAL array, dimension (min(M,N))
 111 *>          The singular values of A in decreasing order.
 112 *>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
 113 *> \endverbatim
 114 *>
 115 *> \param[in] RCOND
 116 *> \verbatim
 117 *>          RCOND is REAL
 118 *>          RCOND is used to determine the effective rank of A.
 119 *>          Singular values S(i) <= RCOND*S(1) are treated as zero.
 120 *>          If RCOND < 0, machine precision is used instead.
 121 *> \endverbatim
 122 *>
 123 *> \param[out] RANK
 124 *> \verbatim
 125 *>          RANK is INTEGER
 126 *>          The effective rank of A, i.e., the number of singular values
 127 *>          which are greater than RCOND*S(1).
 128 *> \endverbatim
 129 *>
 130 *> \param[out] WORK
 131 *> \verbatim
 132 *>          WORK is REAL array, dimension (MAX(1,LWORK))
 133 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] LWORK
 137 *> \verbatim
 138 *>          LWORK is INTEGER
 139 *>          The dimension of the array WORK. LWORK >= 1, and also:
 140 *>          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
 141 *>          For good performance, LWORK should generally be larger.
 142 *>
 143 *>          If LWORK = -1, then a workspace query is assumed; the routine
 144 *>          only calculates the optimal size of the WORK array, returns
 145 *>          this value as the first entry of the WORK array, and no error
 146 *>          message related to LWORK is issued by XERBLA.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] INFO
 150 *> \verbatim
 151 *>          INFO is INTEGER
 152 *>          = 0:  successful exit
 153 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 154 *>          > 0:  the algorithm for computing the SVD failed to converge;
 155 *>                if INFO = i, i off-diagonal elements of an intermediate
 156 *>                bidiagonal form did not converge to zero.
 157 *> \endverbatim
 158 *
 159 *  Authors:
 160 *  ========
 161 *
 162 *> \author Univ. of Tennessee
 163 *> \author Univ. of California Berkeley
 164 *> \author Univ. of Colorado Denver
 165 *> \author NAG Ltd.
 166 *
 167 *> \date December 2016
 168 *
 169 *> \ingroup realGEsolve
 170 *
 171 *  =====================================================================
 172       SUBROUTINE SGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
 173      $                   WORK, LWORK, INFO )
 174 *
 175 *  -- LAPACK driver routine (version 3.7.0) --
 176 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 177 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 178 *     December 2016
 179 *
 180 *     .. Scalar Arguments ..
 181       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
 182       REAL               RCOND
 183 *     ..
 184 *     .. Array Arguments ..
 185       REAL               A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
 186 *     ..
 187 *
 188 *  =====================================================================
 189 *
 190 *     .. Parameters ..
 191       REAL               ZERO, ONE
 192       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 193 *     ..
 194 *     .. Local Scalars ..
 195       LOGICAL            LQUERY
 196       INTEGER            BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
 197      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
 198      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
 199       INTEGER            LWORK_SGEQRF, LWORK_SORMQR, LWORK_SGEBRD,
 200      $                   LWORK_SORMBR, LWORK_SORGBR, LWORK_SORMLQ
 201       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
 202 *     ..
 203 *     .. Local Arrays ..
 204       REAL               DUM( 1 )
 205 *     ..
 206 *     .. External Subroutines ..
 207       EXTERNAL           SBDSQR, SCOPY, SGEBRD, SGELQF, SGEMM, SGEMV,
 208      $                   SGEQRF, SLABAD, SLACPY, SLASCL, SLASET, SORGBR,
 209      $                   SORMBR, SORMLQ, SORMQR, SRSCL, XERBLA
 210 *     ..
 211 *     .. External Functions ..
 212       INTEGER            ILAENV
 213       REAL               SLAMCH, SLANGE
 214       EXTERNAL           ILAENV, SLAMCH, SLANGE
 215 *     ..
 216 *     .. Intrinsic Functions ..
 217       INTRINSIC          MAX, MIN
 218 *     ..
 219 *     .. Executable Statements ..
 220 *
 221 *     Test the input arguments
 222 *
 223       INFO = 0
 224       MINMN = MIN( M, N )
 225       MAXMN = MAX( M, N )
 226       LQUERY = ( LWORK.EQ.-1 )
 227       IF( M.LT.0 ) THEN
 228          INFO = -1
 229       ELSE IF( N.LT.0 ) THEN
 230          INFO = -2
 231       ELSE IF( NRHS.LT.0 ) THEN
 232          INFO = -3
 233       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 234          INFO = -5
 235       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
 236          INFO = -7
 237       END IF
 238 *
 239 *     Compute workspace
 240 *      (Note: Comments in the code beginning "Workspace:" describe the
 241 *       minimal amount of workspace needed at that point in the code,
 242 *       as well as the preferred amount for good performance.
 243 *       NB refers to the optimal block size for the immediately
 244 *       following subroutine, as returned by ILAENV.)
 245 *
 246       IF( INFO.EQ.0 ) THEN
 247          MINWRK = 1
 248          MAXWRK = 1
 249          IF( MINMN.GT.0 ) THEN
 250             MM = M
 251             MNTHR = ILAENV( 6, 'SGELSS', ' ', M, N, NRHS, -1 )
 252             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
 253 *
 254 *              Path 1a - overdetermined, with many more rows than
 255 *                        columns
 256 *
 257 *              Compute space needed for SGEQRF
 258                CALL SGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
 259                LWORK_SGEQRF=DUM(1)
 260 *              Compute space needed for SORMQR
 261                CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
 262      $                   LDB, DUM(1), -1, INFO )
 263                LWORK_SORMQR=DUM(1)
 264                MM = N
 265                MAXWRK = MAX( MAXWRK, N + LWORK_SGEQRF )
 266                MAXWRK = MAX( MAXWRK, N + LWORK_SORMQR )
 267             END IF
 268             IF( M.GE.N ) THEN
 269 *
 270 *              Path 1 - overdetermined or exactly determined
 271 *
 272 *              Compute workspace needed for SBDSQR
 273 *
 274                BDSPAC = MAX( 1, 5*N )
 275 *              Compute space needed for SGEBRD
 276                CALL SGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
 277      $                      DUM(1), DUM(1), -1, INFO )
 278                LWORK_SGEBRD=DUM(1)
 279 *              Compute space needed for SORMBR
 280                CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
 281      $                B, LDB, DUM(1), -1, INFO )
 282                LWORK_SORMBR=DUM(1)
 283 *              Compute space needed for SORGBR
 284                CALL SORGBR( 'P', N, N, N, A, LDA, DUM(1),
 285      $                   DUM(1), -1, INFO )
 286                LWORK_SORGBR=DUM(1)
 287 *              Compute total workspace needed
 288                MAXWRK = MAX( MAXWRK, 3*N + LWORK_SGEBRD )
 289                MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORMBR )
 290                MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORGBR )
 291                MAXWRK = MAX( MAXWRK, BDSPAC )
 292                MAXWRK = MAX( MAXWRK, N*NRHS )
 293                MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC )
 294                MAXWRK = MAX( MINWRK, MAXWRK )
 295             END IF
 296             IF( N.GT.M ) THEN
 297 *
 298 *              Compute workspace needed for SBDSQR
 299 *
 300                BDSPAC = MAX( 1, 5*M )
 301                MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
 302                IF( N.GE.MNTHR ) THEN
 303 *
 304 *                 Path 2a - underdetermined, with many more columns
 305 *                 than rows
 306 *
 307 *                 Compute space needed for SGEBRD
 308                   CALL SGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
 309      $                      DUM(1), DUM(1), -1, INFO )
 310                   LWORK_SGEBRD=DUM(1)
 311 *                 Compute space needed for SORMBR
 312                   CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA,
 313      $                DUM(1), B, LDB, DUM(1), -1, INFO )
 314                   LWORK_SORMBR=DUM(1)
 315 *                 Compute space needed for SORGBR
 316                   CALL SORGBR( 'P', M, M, M, A, LDA, DUM(1),
 317      $                   DUM(1), -1, INFO )
 318                   LWORK_SORGBR=DUM(1)
 319 *                 Compute space needed for SORMLQ
 320                   CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
 321      $                 B, LDB, DUM(1), -1, INFO )
 322                   LWORK_SORMLQ=DUM(1)
 323 *                 Compute total workspace needed
 324                   MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1,
 325      $                                  -1 )
 326                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SGEBRD )
 327                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORMBR )
 328                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORGBR )
 329                   MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC )
 330                   IF( NRHS.GT.1 ) THEN
 331                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
 332                   ELSE
 333                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
 334                   END IF
 335                   MAXWRK = MAX( MAXWRK, M + LWORK_SORMLQ )
 336                ELSE
 337 *
 338 *                 Path 2 - underdetermined
 339 *
 340 *                 Compute space needed for SGEBRD
 341                   CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
 342      $                      DUM(1), DUM(1), -1, INFO )
 343                   LWORK_SGEBRD=DUM(1)
 344 *                 Compute space needed for SORMBR
 345                   CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA,
 346      $                DUM(1), B, LDB, DUM(1), -1, INFO )
 347                   LWORK_SORMBR=DUM(1)
 348 *                 Compute space needed for SORGBR
 349                   CALL SORGBR( 'P', M, N, M, A, LDA, DUM(1),
 350      $                   DUM(1), -1, INFO )
 351                   LWORK_SORGBR=DUM(1)
 352                   MAXWRK = 3*M + LWORK_SGEBRD
 353                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORMBR )
 354                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORGBR )
 355                   MAXWRK = MAX( MAXWRK, BDSPAC )
 356                   MAXWRK = MAX( MAXWRK, N*NRHS )
 357                END IF
 358             END IF
 359             MAXWRK = MAX( MINWRK, MAXWRK )
 360          END IF
 361          WORK( 1 ) = MAXWRK
 362 *
 363          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
 364      $      INFO = -12
 365       END IF
 366 *
 367       IF( INFO.NE.0 ) THEN
 368          CALL XERBLA( 'SGELSS', -INFO )
 369          RETURN
 370       ELSE IF( LQUERY ) THEN
 371          RETURN
 372       END IF
 373 *
 374 *     Quick return if possible
 375 *
 376       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 377          RANK = 0
 378          RETURN
 379       END IF
 380 *
 381 *     Get machine parameters
 382 *
 383       EPS = SLAMCH( 'P' )
 384       SFMIN = SLAMCH( 'S' )
 385       SMLNUM = SFMIN / EPS
 386       BIGNUM = ONE / SMLNUM
 387       CALL SLABAD( SMLNUM, BIGNUM )
 388 *
 389 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 390 *
 391       ANRM = SLANGE( 'M', M, N, A, LDA, WORK )
 392       IASCL = 0
 393       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 394 *
 395 *        Scale matrix norm up to SMLNUM
 396 *
 397          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
 398          IASCL = 1
 399       ELSE IF( ANRM.GT.BIGNUM ) THEN
 400 *
 401 *        Scale matrix norm down to BIGNUM
 402 *
 403          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
 404          IASCL = 2
 405       ELSE IF( ANRM.EQ.ZERO ) THEN
 406 *
 407 *        Matrix all zero. Return zero solution.
 408 *
 409          CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
 410          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
 411          RANK = 0
 412          GO TO 70
 413       END IF
 414 *
 415 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 416 *
 417       BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK )
 418       IBSCL = 0
 419       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 420 *
 421 *        Scale matrix norm up to SMLNUM
 422 *
 423          CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
 424          IBSCL = 1
 425       ELSE IF( BNRM.GT.BIGNUM ) THEN
 426 *
 427 *        Scale matrix norm down to BIGNUM
 428 *
 429          CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
 430          IBSCL = 2
 431       END IF
 432 *
 433 *     Overdetermined case
 434 *
 435       IF( M.GE.N ) THEN
 436 *
 437 *        Path 1 - overdetermined or exactly determined
 438 *
 439          MM = M
 440          IF( M.GE.MNTHR ) THEN
 441 *
 442 *           Path 1a - overdetermined, with many more rows than columns
 443 *
 444             MM = N
 445             ITAU = 1
 446             IWORK = ITAU + N
 447 *
 448 *           Compute A=Q*R
 449 *           (Workspace: need 2*N, prefer N+N*NB)
 450 *
 451             CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
 452      $                   LWORK-IWORK+1, INFO )
 453 *
 454 *           Multiply B by transpose(Q)
 455 *           (Workspace: need N+NRHS, prefer N+NRHS*NB)
 456 *
 457             CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
 458      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 459 *
 460 *           Zero out below R
 461 *
 462             IF( N.GT.1 )
 463      $         CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
 464          END IF
 465 *
 466          IE = 1
 467          ITAUQ = IE + N
 468          ITAUP = ITAUQ + N
 469          IWORK = ITAUP + N
 470 *
 471 *        Bidiagonalize R in A
 472 *        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
 473 *
 474          CALL SGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 475      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
 476      $                INFO )
 477 *
 478 *        Multiply B by transpose of left bidiagonalizing vectors of R
 479 *        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
 480 *
 481          CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
 482      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 483 *
 484 *        Generate right bidiagonalizing vectors of R in A
 485 *        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
 486 *
 487          CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
 488      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 489          IWORK = IE + N
 490 *
 491 *        Perform bidiagonal QR iteration
 492 *          multiply B by transpose of left singular vectors
 493 *          compute right singular vectors in A
 494 *        (Workspace: need BDSPAC)
 495 *
 496          CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
 497      $                1, B, LDB, WORK( IWORK ), INFO )
 498          IF( INFO.NE.0 )
 499      $      GO TO 70
 500 *
 501 *        Multiply B by reciprocals of singular values
 502 *
 503          THR = MAX( RCOND*S( 1 ), SFMIN )
 504          IF( RCOND.LT.ZERO )
 505      $      THR = MAX( EPS*S( 1 ), SFMIN )
 506          RANK = 0
 507          DO 10 I = 1, N
 508             IF( S( I ).GT.THR ) THEN
 509                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 510                RANK = RANK + 1
 511             ELSE
 512                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
 513             END IF
 514    10    CONTINUE
 515 *
 516 *        Multiply B by right singular vectors
 517 *        (Workspace: need N, prefer N*NRHS)
 518 *
 519          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
 520             CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
 521      $                  WORK, LDB )
 522             CALL SLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
 523          ELSE IF( NRHS.GT.1 ) THEN
 524             CHUNK = LWORK / N
 525             DO 20 I = 1, NRHS, CHUNK
 526                BL = MIN( NRHS-I+1, CHUNK )
 527                CALL SGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ),
 528      $                     LDB, ZERO, WORK, N )
 529                CALL SLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
 530    20       CONTINUE
 531          ELSE
 532             CALL SGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
 533             CALL SCOPY( N, WORK, 1, B, 1 )
 534          END IF
 535 *
 536       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
 537      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
 538 *
 539 *        Path 2a - underdetermined, with many more columns than rows
 540 *        and sufficient workspace for an efficient algorithm
 541 *
 542          LDWORK = M
 543          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
 544      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
 545          ITAU = 1
 546          IWORK = M + 1
 547 *
 548 *        Compute A=L*Q
 549 *        (Workspace: need 2*M, prefer M+M*NB)
 550 *
 551          CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
 552      $                LWORK-IWORK+1, INFO )
 553          IL = IWORK
 554 *
 555 *        Copy L to WORK(IL), zeroing out above it
 556 *
 557          CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
 558          CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
 559      $                LDWORK )
 560          IE = IL + LDWORK*M
 561          ITAUQ = IE + M
 562          ITAUP = ITAUQ + M
 563          IWORK = ITAUP + M
 564 *
 565 *        Bidiagonalize L in WORK(IL)
 566 *        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
 567 *
 568          CALL SGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
 569      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
 570      $                LWORK-IWORK+1, INFO )
 571 *
 572 *        Multiply B by transpose of left bidiagonalizing vectors of L
 573 *        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
 574 *
 575          CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
 576      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
 577      $                LWORK-IWORK+1, INFO )
 578 *
 579 *        Generate right bidiagonalizing vectors of R in WORK(IL)
 580 *        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
 581 *
 582          CALL SORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
 583      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 584          IWORK = IE + M
 585 *
 586 *        Perform bidiagonal QR iteration,
 587 *           computing right singular vectors of L in WORK(IL) and
 588 *           multiplying B by transpose of left singular vectors
 589 *        (Workspace: need M*M+M+BDSPAC)
 590 *
 591          CALL SBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ),
 592      $                LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO )
 593          IF( INFO.NE.0 )
 594      $      GO TO 70
 595 *
 596 *        Multiply B by reciprocals of singular values
 597 *
 598          THR = MAX( RCOND*S( 1 ), SFMIN )
 599          IF( RCOND.LT.ZERO )
 600      $      THR = MAX( EPS*S( 1 ), SFMIN )
 601          RANK = 0
 602          DO 30 I = 1, M
 603             IF( S( I ).GT.THR ) THEN
 604                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 605                RANK = RANK + 1
 606             ELSE
 607                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
 608             END IF
 609    30    CONTINUE
 610          IWORK = IE
 611 *
 612 *        Multiply B by right singular vectors of L in WORK(IL)
 613 *        (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
 614 *
 615          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
 616             CALL SGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK,
 617      $                  B, LDB, ZERO, WORK( IWORK ), LDB )
 618             CALL SLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
 619          ELSE IF( NRHS.GT.1 ) THEN
 620             CHUNK = ( LWORK-IWORK+1 ) / M
 621             DO 40 I = 1, NRHS, CHUNK
 622                BL = MIN( NRHS-I+1, CHUNK )
 623                CALL SGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK,
 624      $                     B( 1, I ), LDB, ZERO, WORK( IWORK ), M )
 625                CALL SLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
 626      $                      LDB )
 627    40       CONTINUE
 628          ELSE
 629             CALL SGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ),
 630      $                  1, ZERO, WORK( IWORK ), 1 )
 631             CALL SCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
 632          END IF
 633 *
 634 *        Zero out below first M rows of B
 635 *
 636          CALL SLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
 637          IWORK = ITAU + M
 638 *
 639 *        Multiply transpose(Q) by B
 640 *        (Workspace: need M+NRHS, prefer M+NRHS*NB)
 641 *
 642          CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
 643      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 644 *
 645       ELSE
 646 *
 647 *        Path 2 - remaining underdetermined cases
 648 *
 649          IE = 1
 650          ITAUQ = IE + M
 651          ITAUP = ITAUQ + M
 652          IWORK = ITAUP + M
 653 *
 654 *        Bidiagonalize A
 655 *        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
 656 *
 657          CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 658      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
 659      $                INFO )
 660 *
 661 *        Multiply B by transpose of left bidiagonalizing vectors
 662 *        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
 663 *
 664          CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
 665      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 666 *
 667 *        Generate right bidiagonalizing vectors in A
 668 *        (Workspace: need 4*M, prefer 3*M+M*NB)
 669 *
 670          CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
 671      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 672          IWORK = IE + M
 673 *
 674 *        Perform bidiagonal QR iteration,
 675 *           computing right singular vectors of A in A and
 676 *           multiplying B by transpose of left singular vectors
 677 *        (Workspace: need BDSPAC)
 678 *
 679          CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
 680      $                1, B, LDB, WORK( IWORK ), INFO )
 681          IF( INFO.NE.0 )
 682      $      GO TO 70
 683 *
 684 *        Multiply B by reciprocals of singular values
 685 *
 686          THR = MAX( RCOND*S( 1 ), SFMIN )
 687          IF( RCOND.LT.ZERO )
 688      $      THR = MAX( EPS*S( 1 ), SFMIN )
 689          RANK = 0
 690          DO 50 I = 1, M
 691             IF( S( I ).GT.THR ) THEN
 692                CALL SRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 693                RANK = RANK + 1
 694             ELSE
 695                CALL SLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
 696             END IF
 697    50    CONTINUE
 698 *
 699 *        Multiply B by right singular vectors of A
 700 *        (Workspace: need N, prefer N*NRHS)
 701 *
 702          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
 703             CALL SGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO,
 704      $                  WORK, LDB )
 705             CALL SLACPY( 'F', N, NRHS, WORK, LDB, B, LDB )
 706          ELSE IF( NRHS.GT.1 ) THEN
 707             CHUNK = LWORK / N
 708             DO 60 I = 1, NRHS, CHUNK
 709                BL = MIN( NRHS-I+1, CHUNK )
 710                CALL SGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ),
 711      $                     LDB, ZERO, WORK, N )
 712                CALL SLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
 713    60       CONTINUE
 714          ELSE
 715             CALL SGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
 716             CALL SCOPY( N, WORK, 1, B, 1 )
 717          END IF
 718       END IF
 719 *
 720 *     Undo scaling
 721 *
 722       IF( IASCL.EQ.1 ) THEN
 723          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
 724          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
 725      $                INFO )
 726       ELSE IF( IASCL.EQ.2 ) THEN
 727          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
 728          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
 729      $                INFO )
 730       END IF
 731       IF( IBSCL.EQ.1 ) THEN
 732          CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
 733       ELSE IF( IBSCL.EQ.2 ) THEN
 734          CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
 735       END IF
 736 *
 737    70 CONTINUE
 738       WORK( 1 ) = MAXWRK
 739       RETURN
 740 *
 741 *     End of SGELSS
 742 *
 743       END