LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgels.f
Go to the documentation of this file.
00001 *> \brief <b> DGELS 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 DGELS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgels.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgels.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgels.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
00022 *                         INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          TRANS
00026 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> DGELS solves overdetermined or underdetermined real linear systems
00039 *> involving an M-by-N matrix A, or its transpose, using a QR or LQ
00040 *> factorization of A.  It is assumed that A has full rank.
00041 *>
00042 *> The following options are provided:
00043 *>
00044 *> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
00045 *>    an overdetermined system, i.e., solve the least squares problem
00046 *>                 minimize || B - A*X ||.
00047 *>
00048 *> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
00049 *>    an underdetermined system A * X = B.
00050 *>
00051 *> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
00052 *>    an undetermined system A**T * X = B.
00053 *>
00054 *> 4. If TRANS = 'T' and m < n:  find the least squares solution of
00055 *>    an overdetermined system, i.e., solve the least squares problem
00056 *>                 minimize || B - A**T * X ||.
00057 *>
00058 *> Several right hand side vectors b and solution vectors x can be
00059 *> handled in a single call; they are stored as the columns of the
00060 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00061 *> matrix X.
00062 *> \endverbatim
00063 *
00064 *  Arguments:
00065 *  ==========
00066 *
00067 *> \param[in] TRANS
00068 *> \verbatim
00069 *>          TRANS is CHARACTER*1
00070 *>          = 'N': the linear system involves A;
00071 *>          = 'T': the linear system involves A**T.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] M
00075 *> \verbatim
00076 *>          M is INTEGER
00077 *>          The number of rows of the matrix A.  M >= 0.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] N
00081 *> \verbatim
00082 *>          N is INTEGER
00083 *>          The number of columns of the matrix A.  N >= 0.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] NRHS
00087 *> \verbatim
00088 *>          NRHS is INTEGER
00089 *>          The number of right hand sides, i.e., the number of
00090 *>          columns of the matrices B and X. NRHS >=0.
00091 *> \endverbatim
00092 *>
00093 *> \param[in,out] A
00094 *> \verbatim
00095 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00096 *>          On entry, the M-by-N matrix A.
00097 *>          On exit,
00098 *>            if M >= N, A is overwritten by details of its QR
00099 *>                       factorization as returned by DGEQRF;
00100 *>            if M <  N, A is overwritten by details of its LQ
00101 *>                       factorization as returned by DGELQF.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] LDA
00105 *> \verbatim
00106 *>          LDA is INTEGER
00107 *>          The leading dimension of the array A.  LDA >= max(1,M).
00108 *> \endverbatim
00109 *>
00110 *> \param[in,out] B
00111 *> \verbatim
00112 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
00113 *>          On entry, the matrix B of right hand side vectors, stored
00114 *>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
00115 *>          if TRANS = 'T'.
00116 *>          On exit, if INFO = 0, B is overwritten by the solution
00117 *>          vectors, stored columnwise:
00118 *>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
00119 *>          squares solution vectors; the residual sum of squares for the
00120 *>          solution in each column is given by the sum of squares of
00121 *>          elements N+1 to M in that column;
00122 *>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
00123 *>          minimum norm solution vectors;
00124 *>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
00125 *>          minimum norm solution vectors;
00126 *>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
00127 *>          least squares solution vectors; the residual sum of squares
00128 *>          for the solution in each column is given by the sum of
00129 *>          squares of elements M+1 to N in that column.
00130 *> \endverbatim
00131 *>
00132 *> \param[in] LDB
00133 *> \verbatim
00134 *>          LDB is INTEGER
00135 *>          The leading dimension of the array B. LDB >= MAX(1,M,N).
00136 *> \endverbatim
00137 *>
00138 *> \param[out] WORK
00139 *> \verbatim
00140 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00141 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00142 *> \endverbatim
00143 *>
00144 *> \param[in] LWORK
00145 *> \verbatim
00146 *>          LWORK is INTEGER
00147 *>          The dimension of the array WORK.
00148 *>          LWORK >= max( 1, MN + max( MN, NRHS ) ).
00149 *>          For optimal performance,
00150 *>          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
00151 *>          where MN = min(M,N) and NB is the optimum block size.
00152 *>
00153 *>          If LWORK = -1, then a workspace query is assumed; the routine
00154 *>          only calculates the optimal size of the WORK array, returns
00155 *>          this value as the first entry of the WORK array, and no error
00156 *>          message related to LWORK is issued by XERBLA.
00157 *> \endverbatim
00158 *>
00159 *> \param[out] INFO
00160 *> \verbatim
00161 *>          INFO is INTEGER
00162 *>          = 0:  successful exit
00163 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00164 *>          > 0:  if INFO =  i, the i-th diagonal element of the
00165 *>                triangular factor of A is zero, so that A does not have
00166 *>                full rank; the least squares solution could not be
00167 *>                computed.
00168 *> \endverbatim
00169 *
00170 *  Authors:
00171 *  ========
00172 *
00173 *> \author Univ. of Tennessee 
00174 *> \author Univ. of California Berkeley 
00175 *> \author Univ. of Colorado Denver 
00176 *> \author NAG Ltd. 
00177 *
00178 *> \date November 2011
00179 *
00180 *> \ingroup doubleGEsolve
00181 *
00182 *  =====================================================================
00183       SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
00184      $                  INFO )
00185 *
00186 *  -- LAPACK driver routine (version 3.4.0) --
00187 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00188 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00189 *     November 2011
00190 *
00191 *     .. Scalar Arguments ..
00192       CHARACTER          TRANS
00193       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
00194 *     ..
00195 *     .. Array Arguments ..
00196       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
00197 *     ..
00198 *
00199 *  =====================================================================
00200 *
00201 *     .. Parameters ..
00202       DOUBLE PRECISION   ZERO, ONE
00203       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00204 *     ..
00205 *     .. Local Scalars ..
00206       LOGICAL            LQUERY, TPSD
00207       INTEGER            BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
00208       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM
00209 *     ..
00210 *     .. Local Arrays ..
00211       DOUBLE PRECISION   RWORK( 1 )
00212 *     ..
00213 *     .. External Functions ..
00214       LOGICAL            LSAME
00215       INTEGER            ILAENV
00216       DOUBLE PRECISION   DLAMCH, DLANGE
00217       EXTERNAL           LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
00218 *     ..
00219 *     .. External Subroutines ..
00220       EXTERNAL           DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR,
00221      $                   DTRTRS, XERBLA
00222 *     ..
00223 *     .. Intrinsic Functions ..
00224       INTRINSIC          DBLE, MAX, MIN
00225 *     ..
00226 *     .. Executable Statements ..
00227 *
00228 *     Test the input arguments.
00229 *
00230       INFO = 0
00231       MN = MIN( M, N )
00232       LQUERY = ( LWORK.EQ.-1 )
00233       IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
00234          INFO = -1
00235       ELSE IF( M.LT.0 ) THEN
00236          INFO = -2
00237       ELSE IF( N.LT.0 ) THEN
00238          INFO = -3
00239       ELSE IF( NRHS.LT.0 ) THEN
00240          INFO = -4
00241       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00242          INFO = -6
00243       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00244          INFO = -8
00245       ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
00246      $          THEN
00247          INFO = -10
00248       END IF
00249 *
00250 *     Figure out optimal block size
00251 *
00252       IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
00253 *
00254          TPSD = .TRUE.
00255          IF( LSAME( TRANS, 'N' ) )
00256      $      TPSD = .FALSE.
00257 *
00258          IF( M.GE.N ) THEN
00259             NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00260             IF( TPSD ) THEN
00261                NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N,
00262      $              -1 ) )
00263             ELSE
00264                NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N,
00265      $              -1 ) )
00266             END IF
00267          ELSE
00268             NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
00269             IF( TPSD ) THEN
00270                NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M,
00271      $              -1 ) )
00272             ELSE
00273                NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M,
00274      $              -1 ) )
00275             END IF
00276          END IF
00277 *
00278          WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
00279          WORK( 1 ) = DBLE( WSIZE )
00280 *
00281       END IF
00282 *
00283       IF( INFO.NE.0 ) THEN
00284          CALL XERBLA( 'DGELS ', -INFO )
00285          RETURN
00286       ELSE IF( LQUERY ) THEN
00287          RETURN
00288       END IF
00289 *
00290 *     Quick return if possible
00291 *
00292       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00293          CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00294          RETURN
00295       END IF
00296 *
00297 *     Get machine parameters
00298 *
00299       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
00300       BIGNUM = ONE / SMLNUM
00301       CALL DLABAD( SMLNUM, BIGNUM )
00302 *
00303 *     Scale A, B if max element outside range [SMLNUM,BIGNUM]
00304 *
00305       ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
00306       IASCL = 0
00307       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00308 *
00309 *        Scale matrix norm up to SMLNUM
00310 *
00311          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00312          IASCL = 1
00313       ELSE IF( ANRM.GT.BIGNUM ) THEN
00314 *
00315 *        Scale matrix norm down to BIGNUM
00316 *
00317          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00318          IASCL = 2
00319       ELSE IF( ANRM.EQ.ZERO ) THEN
00320 *
00321 *        Matrix all zero. Return zero solution.
00322 *
00323          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00324          GO TO 50
00325       END IF
00326 *
00327       BROW = M
00328       IF( TPSD )
00329      $   BROW = N
00330       BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
00331       IBSCL = 0
00332       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00333 *
00334 *        Scale matrix norm up to SMLNUM
00335 *
00336          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
00337      $                INFO )
00338          IBSCL = 1
00339       ELSE IF( BNRM.GT.BIGNUM ) THEN
00340 *
00341 *        Scale matrix norm down to BIGNUM
00342 *
00343          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
00344      $                INFO )
00345          IBSCL = 2
00346       END IF
00347 *
00348       IF( M.GE.N ) THEN
00349 *
00350 *        compute QR factorization of A
00351 *
00352          CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
00353      $                INFO )
00354 *
00355 *        workspace at least N, optimally N*NB
00356 *
00357          IF( .NOT.TPSD ) THEN
00358 *
00359 *           Least-Squares Problem min || A * X - B ||
00360 *
00361 *           B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
00362 *
00363             CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
00364      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
00365      $                   INFO )
00366 *
00367 *           workspace at least NRHS, optimally NRHS*NB
00368 *
00369 *           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
00370 *
00371             CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
00372      $                   A, LDA, B, LDB, INFO )
00373 *
00374             IF( INFO.GT.0 ) THEN
00375                RETURN
00376             END IF
00377 *
00378             SCLLEN = N
00379 *
00380          ELSE
00381 *
00382 *           Overdetermined system of equations A**T * X = B
00383 *
00384 *           B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
00385 *
00386             CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
00387      $                   A, LDA, B, LDB, INFO )
00388 *
00389             IF( INFO.GT.0 ) THEN
00390                RETURN
00391             END IF
00392 *
00393 *           B(N+1:M,1:NRHS) = ZERO
00394 *
00395             DO 20 J = 1, NRHS
00396                DO 10 I = N + 1, M
00397                   B( I, J ) = ZERO
00398    10          CONTINUE
00399    20       CONTINUE
00400 *
00401 *           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
00402 *
00403             CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
00404      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
00405      $                   INFO )
00406 *
00407 *           workspace at least NRHS, optimally NRHS*NB
00408 *
00409             SCLLEN = M
00410 *
00411          END IF
00412 *
00413       ELSE
00414 *
00415 *        Compute LQ factorization of A
00416 *
00417          CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
00418      $                INFO )
00419 *
00420 *        workspace at least M, optimally M*NB.
00421 *
00422          IF( .NOT.TPSD ) THEN
00423 *
00424 *           underdetermined system of equations A * X = B
00425 *
00426 *           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
00427 *
00428             CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
00429      $                   A, LDA, B, LDB, INFO )
00430 *
00431             IF( INFO.GT.0 ) THEN
00432                RETURN
00433             END IF
00434 *
00435 *           B(M+1:N,1:NRHS) = 0
00436 *
00437             DO 40 J = 1, NRHS
00438                DO 30 I = M + 1, N
00439                   B( I, J ) = ZERO
00440    30          CONTINUE
00441    40       CONTINUE
00442 *
00443 *           B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
00444 *
00445             CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
00446      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
00447      $                   INFO )
00448 *
00449 *           workspace at least NRHS, optimally NRHS*NB
00450 *
00451             SCLLEN = N
00452 *
00453          ELSE
00454 *
00455 *           overdetermined system min || A**T * X - B ||
00456 *
00457 *           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
00458 *
00459             CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
00460      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
00461      $                   INFO )
00462 *
00463 *           workspace at least NRHS, optimally NRHS*NB
00464 *
00465 *           B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
00466 *
00467             CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
00468      $                   A, LDA, B, LDB, INFO )
00469 *
00470             IF( INFO.GT.0 ) THEN
00471                RETURN
00472             END IF
00473 *
00474             SCLLEN = M
00475 *
00476          END IF
00477 *
00478       END IF
00479 *
00480 *     Undo scaling
00481 *
00482       IF( IASCL.EQ.1 ) THEN
00483          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
00484      $                INFO )
00485       ELSE IF( IASCL.EQ.2 ) THEN
00486          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
00487      $                INFO )
00488       END IF
00489       IF( IBSCL.EQ.1 ) THEN
00490          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
00491      $                INFO )
00492       ELSE IF( IBSCL.EQ.2 ) THEN
00493          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
00494      $                INFO )
00495       END IF
00496 *
00497    50 CONTINUE
00498       WORK( 1 ) = DBLE( WSIZE )
00499 *
00500       RETURN
00501 *
00502 *     End of DGELS
00503 *
00504       END
 All Files Functions