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