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