LAPACK: SRC/dgels.f Source File

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