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