LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgelss.f
Go to the documentation of this file.
00001 *> \brief <b> DGELSS 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 DGELSS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelss.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelss.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelss.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGELSS( 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 *       DOUBLE PRECISION   RCOND
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> DGELSS 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 doubleGEsolve
00170 *
00171 *  =====================================================================
00172       SUBROUTINE DGELSS( 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       DOUBLE PRECISION   RCOND
00183 *     ..
00184 *     .. Array Arguments ..
00185       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
00186 *     ..
00187 *
00188 *  =====================================================================
00189 *
00190 *     .. Parameters ..
00191       DOUBLE PRECISION   ZERO, ONE
00192       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD,
00200      $                   LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ,
00201      $                   LWORK_DGELQF
00202       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
00203 *     ..
00204 *     .. Local Arrays ..
00205       DOUBLE PRECISION   DUM( 1 )
00206 *     ..
00207 *     .. External Subroutines ..
00208       EXTERNAL           DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV,
00209      $                   DGEQRF, DLABAD, DLACPY, DLASCL, DLASET, DORGBR,
00210      $                   DORMBR, DORMLQ, DORMQR, DRSCL, XERBLA
00211 *     ..
00212 *     .. External Functions ..
00213       INTEGER            ILAENV
00214       DOUBLE PRECISION   DLAMCH, DLANGE
00215       EXTERNAL           ILAENV, DLAMCH, DLANGE
00216 *     ..
00217 *     .. Intrinsic Functions ..
00218       INTRINSIC          MAX, MIN
00219 *     ..
00220 *     .. Executable Statements ..
00221 *
00222 *     Test the input arguments
00223 *
00224       INFO = 0
00225       MINMN = MIN( M, N )
00226       MAXMN = MAX( M, N )
00227       LQUERY = ( LWORK.EQ.-1 )
00228       IF( M.LT.0 ) THEN
00229          INFO = -1
00230       ELSE IF( N.LT.0 ) THEN
00231          INFO = -2
00232       ELSE IF( NRHS.LT.0 ) THEN
00233          INFO = -3
00234       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00235          INFO = -5
00236       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
00237          INFO = -7
00238       END IF
00239 *
00240 *     Compute workspace
00241 *      (Note: Comments in the code beginning "Workspace:" describe the
00242 *       minimal amount of workspace needed at that point in the code,
00243 *       as well as the preferred amount for good performance.
00244 *       NB refers to the optimal block size for the immediately
00245 *       following subroutine, as returned by ILAENV.)
00246 *
00247       IF( INFO.EQ.0 ) THEN
00248          MINWRK = 1
00249          MAXWRK = 1
00250          IF( MINMN.GT.0 ) THEN
00251             MM = M
00252             MNTHR = ILAENV( 6, 'DGELSS', ' ', M, N, NRHS, -1 )
00253             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
00254 *
00255 *              Path 1a - overdetermined, with many more rows than
00256 *                        columns
00257 *
00258 *              Compute space needed for DGEQRF
00259                CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
00260                LWORK_DGEQRF=DUM(1)
00261 *              Compute space needed for DORMQR
00262                CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
00263      $                   LDB, DUM(1), -1, INFO )
00264                LWORK_DORMQR=DUM(1)
00265                MM = N
00266                MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF )
00267                MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR )
00268             END IF
00269             IF( M.GE.N ) THEN
00270 *
00271 *              Path 1 - overdetermined or exactly determined
00272 *
00273 *              Compute workspace needed for DBDSQR
00274 *
00275                BDSPAC = MAX( 1, 5*N )
00276 *              Compute space needed for DGEBRD
00277                CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
00278      $                      DUM(1), DUM(1), -1, INFO )
00279                LWORK_DGEBRD=DUM(1)
00280 *              Compute space needed for DORMBR
00281                CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
00282      $                B, LDB, DUM(1), -1, INFO )
00283                LWORK_DORMBR=DUM(1)
00284 *              Compute space needed for DORGBR
00285                CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
00286      $                   DUM(1), -1, INFO )
00287                LWORK_DORGBR=DUM(1)
00288 *              Compute total workspace needed 
00289                MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
00290                MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR )
00291                MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR )
00292                MAXWRK = MAX( MAXWRK, BDSPAC )
00293                MAXWRK = MAX( MAXWRK, N*NRHS )
00294                MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC )
00295                MAXWRK = MAX( MINWRK, MAXWRK )
00296             END IF
00297             IF( N.GT.M ) THEN
00298 *
00299 *              Compute workspace needed for DBDSQR
00300 *
00301                BDSPAC = MAX( 1, 5*M )
00302                MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
00303                IF( N.GE.MNTHR ) THEN
00304 *
00305 *                 Path 2a - underdetermined, with many more columns
00306 *                 than rows
00307 *
00308 *                 Compute space needed for DGELQF
00309                   CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1),
00310      $                -1, INFO )
00311                   LWORK_DGELQF=DUM(1)
00312 *                 Compute space needed for DGEBRD
00313                   CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
00314      $                      DUM(1), DUM(1), -1, INFO )
00315                   LWORK_DGEBRD=DUM(1)
00316 *                 Compute space needed for DORMBR
00317                   CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, 
00318      $                DUM(1), B, LDB, DUM(1), -1, INFO )
00319                   LWORK_DORMBR=DUM(1)
00320 *                 Compute space needed for DORGBR
00321                   CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1),
00322      $                   DUM(1), -1, INFO )
00323                   LWORK_DORGBR=DUM(1)
00324 *                 Compute space needed for DORMLQ
00325                   CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
00326      $                 B, LDB, DUM(1), -1, INFO )
00327                   LWORK_DORMLQ=DUM(1)
00328 *                 Compute total workspace needed 
00329                   MAXWRK = M + LWORK_DGELQF
00330                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD )
00331                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORMBR )
00332                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORGBR )
00333                   MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC )
00334                   IF( NRHS.GT.1 ) THEN
00335                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
00336                   ELSE
00337                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
00338                   END IF
00339                   MAXWRK = MAX( MAXWRK, M + LWORK_DORMLQ )
00340                ELSE
00341 *
00342 *                 Path 2 - underdetermined
00343 *
00344 *                 Compute space needed for DGEBRD
00345                   CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
00346      $                      DUM(1), DUM(1), -1, INFO )
00347                   LWORK_DGEBRD=DUM(1)
00348 *                 Compute space needed for DORMBR
00349                   CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, 
00350      $                DUM(1), B, LDB, DUM(1), -1, INFO )
00351                   LWORK_DORMBR=DUM(1)
00352 *                 Compute space needed for DORGBR
00353                   CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1),
00354      $                   DUM(1), -1, INFO )
00355                   LWORK_DORGBR=DUM(1)
00356                   MAXWRK = 3*M + LWORK_DGEBRD
00357                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR )
00358                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR )
00359                   MAXWRK = MAX( MAXWRK, BDSPAC )
00360                   MAXWRK = MAX( MAXWRK, N*NRHS )
00361                END IF
00362             END IF
00363             MAXWRK = MAX( MINWRK, MAXWRK )
00364          END IF
00365          WORK( 1 ) = MAXWRK
00366 *
00367          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
00368      $      INFO = -12
00369       END IF
00370 *
00371       IF( INFO.NE.0 ) THEN
00372          CALL XERBLA( 'DGELSS', -INFO )
00373          RETURN
00374       ELSE IF( LQUERY ) THEN
00375          RETURN
00376       END IF
00377 *
00378 *     Quick return if possible
00379 *
00380       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00381          RANK = 0
00382          RETURN
00383       END IF
00384 *
00385 *     Get machine parameters
00386 *
00387       EPS = DLAMCH( 'P' )
00388       SFMIN = DLAMCH( 'S' )
00389       SMLNUM = SFMIN / EPS
00390       BIGNUM = ONE / SMLNUM
00391       CALL DLABAD( SMLNUM, BIGNUM )
00392 *
00393 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00394 *
00395       ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
00396       IASCL = 0
00397       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00398 *
00399 *        Scale matrix norm up to SMLNUM
00400 *
00401          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00402          IASCL = 1
00403       ELSE IF( ANRM.GT.BIGNUM ) THEN
00404 *
00405 *        Scale matrix norm down to BIGNUM
00406 *
00407          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00408          IASCL = 2
00409       ELSE IF( ANRM.EQ.ZERO ) THEN
00410 *
00411 *        Matrix all zero. Return zero solution.
00412 *
00413          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00414          CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
00415          RANK = 0
00416          GO TO 70
00417       END IF
00418 *
00419 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00420 *
00421       BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
00422       IBSCL = 0
00423       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00424 *
00425 *        Scale matrix norm up to SMLNUM
00426 *
00427          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00428          IBSCL = 1
00429       ELSE IF( BNRM.GT.BIGNUM ) THEN
00430 *
00431 *        Scale matrix norm down to BIGNUM
00432 *
00433          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00434          IBSCL = 2
00435       END IF
00436 *
00437 *     Overdetermined case
00438 *
00439       IF( M.GE.N ) THEN
00440 *
00441 *        Path 1 - overdetermined or exactly determined
00442 *
00443          MM = M
00444          IF( M.GE.MNTHR ) THEN
00445 *
00446 *           Path 1a - overdetermined, with many more rows than columns
00447 *
00448             MM = N
00449             ITAU = 1
00450             IWORK = ITAU + N
00451 *
00452 *           Compute A=Q*R
00453 *           (Workspace: need 2*N, prefer N+N*NB)
00454 *
00455             CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00456      $                   LWORK-IWORK+1, INFO )
00457 *
00458 *           Multiply B by transpose(Q)
00459 *           (Workspace: need N+NRHS, prefer N+NRHS*NB)
00460 *
00461             CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
00462      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00463 *
00464 *           Zero out below R
00465 *
00466             IF( N.GT.1 )
00467      $         CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
00468          END IF
00469 *
00470          IE = 1
00471          ITAUQ = IE + N
00472          ITAUP = ITAUQ + N
00473          IWORK = ITAUP + N
00474 *
00475 *        Bidiagonalize R in A
00476 *        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
00477 *
00478          CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00479      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00480      $                INFO )
00481 *
00482 *        Multiply B by transpose of left bidiagonalizing vectors of R
00483 *        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
00484 *
00485          CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
00486      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00487 *
00488 *        Generate right bidiagonalizing vectors of R in A
00489 *        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
00490 *
00491          CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
00492      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00493          IWORK = IE + N
00494 *
00495 *        Perform bidiagonal QR iteration
00496 *          multiply B by transpose of left singular vectors
00497 *          compute right singular vectors in A
00498 *        (Workspace: need BDSPAC)
00499 *
00500          CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
00501      $                1, B, LDB, WORK( IWORK ), INFO )
00502          IF( INFO.NE.0 )
00503      $      GO TO 70
00504 *
00505 *        Multiply B by reciprocals of singular values
00506 *
00507          THR = MAX( RCOND*S( 1 ), SFMIN )
00508          IF( RCOND.LT.ZERO )
00509      $      THR = MAX( EPS*S( 1 ), SFMIN )
00510          RANK = 0
00511          DO 10 I = 1, N
00512             IF( S( I ).GT.THR ) THEN
00513                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00514                RANK = RANK + 1
00515             ELSE
00516                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00517             END IF
00518    10    CONTINUE
00519 *
00520 *        Multiply B by right singular vectors
00521 *        (Workspace: need N, prefer N*NRHS)
00522 *
00523          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00524             CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
00525      $                  WORK, LDB )
00526             CALL DLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
00527          ELSE IF( NRHS.GT.1 ) THEN
00528             CHUNK = LWORK / N
00529             DO 20 I = 1, NRHS, CHUNK
00530                BL = MIN( NRHS-I+1, CHUNK )
00531                CALL DGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ),
00532      $                     LDB, ZERO, WORK, N )
00533                CALL DLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
00534    20       CONTINUE
00535          ELSE
00536             CALL DGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
00537             CALL DCOPY( N, WORK, 1, B, 1 )
00538          END IF
00539 *
00540       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
00541      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
00542 *
00543 *        Path 2a - underdetermined, with many more columns than rows
00544 *        and sufficient workspace for an efficient algorithm
00545 *
00546          LDWORK = M
00547          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
00548      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
00549          ITAU = 1
00550          IWORK = M + 1
00551 *
00552 *        Compute A=L*Q
00553 *        (Workspace: need 2*M, prefer M+M*NB)
00554 *
00555          CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00556      $                LWORK-IWORK+1, INFO )
00557          IL = IWORK
00558 *
00559 *        Copy L to WORK(IL), zeroing out above it
00560 *
00561          CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
00562          CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
00563      $                LDWORK )
00564          IE = IL + LDWORK*M
00565          ITAUQ = IE + M
00566          ITAUP = ITAUQ + M
00567          IWORK = ITAUP + M
00568 *
00569 *        Bidiagonalize L in WORK(IL)
00570 *        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
00571 *
00572          CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
00573      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
00574      $                LWORK-IWORK+1, INFO )
00575 *
00576 *        Multiply B by transpose of left bidiagonalizing vectors of L
00577 *        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
00578 *
00579          CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
00580      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
00581      $                LWORK-IWORK+1, INFO )
00582 *
00583 *        Generate right bidiagonalizing vectors of R in WORK(IL)
00584 *        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
00585 *
00586          CALL DORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
00587      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00588          IWORK = IE + M
00589 *
00590 *        Perform bidiagonal QR iteration,
00591 *           computing right singular vectors of L in WORK(IL) and
00592 *           multiplying B by transpose of left singular vectors
00593 *        (Workspace: need M*M+M+BDSPAC)
00594 *
00595          CALL DBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ),
00596      $                LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO )
00597          IF( INFO.NE.0 )
00598      $      GO TO 70
00599 *
00600 *        Multiply B by reciprocals of singular values
00601 *
00602          THR = MAX( RCOND*S( 1 ), SFMIN )
00603          IF( RCOND.LT.ZERO )
00604      $      THR = MAX( EPS*S( 1 ), SFMIN )
00605          RANK = 0
00606          DO 30 I = 1, M
00607             IF( S( I ).GT.THR ) THEN
00608                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00609                RANK = RANK + 1
00610             ELSE
00611                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00612             END IF
00613    30    CONTINUE
00614          IWORK = IE
00615 *
00616 *        Multiply B by right singular vectors of L in WORK(IL)
00617 *        (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
00618 *
00619          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
00620             CALL DGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK,
00621      $                  B, LDB, ZERO, WORK( IWORK ), LDB )
00622             CALL DLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
00623          ELSE IF( NRHS.GT.1 ) THEN
00624             CHUNK = ( LWORK-IWORK+1 ) / M
00625             DO 40 I = 1, NRHS, CHUNK
00626                BL = MIN( NRHS-I+1, CHUNK )
00627                CALL DGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK,
00628      $                     B( 1, I ), LDB, ZERO, WORK( IWORK ), M )
00629                CALL DLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
00630      $                      LDB )
00631    40       CONTINUE
00632          ELSE
00633             CALL DGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ),
00634      $                  1, ZERO, WORK( IWORK ), 1 )
00635             CALL DCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
00636          END IF
00637 *
00638 *        Zero out below first M rows of B
00639 *
00640          CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
00641          IWORK = ITAU + M
00642 *
00643 *        Multiply transpose(Q) by B
00644 *        (Workspace: need M+NRHS, prefer M+NRHS*NB)
00645 *
00646          CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
00647      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00648 *
00649       ELSE
00650 *
00651 *        Path 2 - remaining underdetermined cases
00652 *
00653          IE = 1
00654          ITAUQ = IE + M
00655          ITAUP = ITAUQ + M
00656          IWORK = ITAUP + M
00657 *
00658 *        Bidiagonalize A
00659 *        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
00660 *
00661          CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00662      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00663      $                INFO )
00664 *
00665 *        Multiply B by transpose of left bidiagonalizing vectors
00666 *        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
00667 *
00668          CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
00669      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00670 *
00671 *        Generate right bidiagonalizing vectors in A
00672 *        (Workspace: need 4*M, prefer 3*M+M*NB)
00673 *
00674          CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
00675      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00676          IWORK = IE + M
00677 *
00678 *        Perform bidiagonal QR iteration,
00679 *           computing right singular vectors of A in A and
00680 *           multiplying B by transpose of left singular vectors
00681 *        (Workspace: need BDSPAC)
00682 *
00683          CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
00684      $                1, B, LDB, WORK( IWORK ), INFO )
00685          IF( INFO.NE.0 )
00686      $      GO TO 70
00687 *
00688 *        Multiply B by reciprocals of singular values
00689 *
00690          THR = MAX( RCOND*S( 1 ), SFMIN )
00691          IF( RCOND.LT.ZERO )
00692      $      THR = MAX( EPS*S( 1 ), SFMIN )
00693          RANK = 0
00694          DO 50 I = 1, M
00695             IF( S( I ).GT.THR ) THEN
00696                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00697                RANK = RANK + 1
00698             ELSE
00699                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00700             END IF
00701    50    CONTINUE
00702 *
00703 *        Multiply B by right singular vectors of A
00704 *        (Workspace: need N, prefer N*NRHS)
00705 *
00706          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00707             CALL DGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO,
00708      $                  WORK, LDB )
00709             CALL DLACPY( 'F', N, NRHS, WORK, LDB, B, LDB )
00710          ELSE IF( NRHS.GT.1 ) THEN
00711             CHUNK = LWORK / N
00712             DO 60 I = 1, NRHS, CHUNK
00713                BL = MIN( NRHS-I+1, CHUNK )
00714                CALL DGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ),
00715      $                     LDB, ZERO, WORK, N )
00716                CALL DLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
00717    60       CONTINUE
00718          ELSE
00719             CALL DGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
00720             CALL DCOPY( N, WORK, 1, B, 1 )
00721          END IF
00722       END IF
00723 *
00724 *     Undo scaling
00725 *
00726       IF( IASCL.EQ.1 ) THEN
00727          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00728          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
00729      $                INFO )
00730       ELSE IF( IASCL.EQ.2 ) THEN
00731          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00732          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
00733      $                INFO )
00734       END IF
00735       IF( IBSCL.EQ.1 ) THEN
00736          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00737       ELSE IF( IBSCL.EQ.2 ) THEN
00738          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00739       END IF
00740 *
00741    70 CONTINUE
00742       WORK( 1 ) = MAXWRK
00743       RETURN
00744 *
00745 *     End of DGELSS
00746 *
00747       END
 All Files Functions