LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgelsx.f
Go to the documentation of this file.
00001 *> \brief <b> ZGELSX 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 ZGELSX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00022 *                          WORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
00026 *       DOUBLE PRECISION   RCOND
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            JPVT( * )
00030 *       DOUBLE PRECISION   RWORK( * )
00031 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> This routine is deprecated and has been replaced by routine ZGELSY.
00041 *>
00042 *> ZGELSX computes the minimum-norm solution to a complex linear least
00043 *> squares problem:
00044 *>     minimize || A * X - B ||
00045 *> using a complete orthogonal factorization of A.  A is an M-by-N
00046 *> matrix which may be rank-deficient.
00047 *>
00048 *> Several right hand side vectors b and solution vectors x can be
00049 *> handled in a single call; they are stored as the columns of the
00050 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00051 *> matrix X.
00052 *>
00053 *> The routine first computes a QR factorization with column pivoting:
00054 *>     A * P = Q * [ R11 R12 ]
00055 *>                 [  0  R22 ]
00056 *> with R11 defined as the largest leading submatrix whose estimated
00057 *> condition number is less than 1/RCOND.  The order of R11, RANK,
00058 *> is the effective rank of A.
00059 *>
00060 *> Then, R22 is considered to be negligible, and R12 is annihilated
00061 *> by unitary transformations from the right, arriving at the
00062 *> complete orthogonal factorization:
00063 *>    A * P = Q * [ T11 0 ] * Z
00064 *>                [  0  0 ]
00065 *> The minimum-norm solution is then
00066 *>    X = P * Z**H [ inv(T11)*Q1**H*B ]
00067 *>                 [        0         ]
00068 *> where Q1 consists of the first RANK columns of Q.
00069 *> \endverbatim
00070 *
00071 *  Arguments:
00072 *  ==========
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 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 *>          On exit, A has been overwritten by details of its
00098 *>          complete orthogonal factorization.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] LDA
00102 *> \verbatim
00103 *>          LDA is INTEGER
00104 *>          The leading dimension of the array A.  LDA >= max(1,M).
00105 *> \endverbatim
00106 *>
00107 *> \param[in,out] B
00108 *> \verbatim
00109 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
00110 *>          On entry, the M-by-NRHS right hand side matrix B.
00111 *>          On exit, the N-by-NRHS solution matrix X.
00112 *>          If m >= n and RANK = n, the residual sum-of-squares for
00113 *>          the solution in the i-th column is given by the sum of
00114 *>          squares of elements N+1:M in that column.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] LDB
00118 *> \verbatim
00119 *>          LDB is INTEGER
00120 *>          The leading dimension of the array B. LDB >= max(1,M,N).
00121 *> \endverbatim
00122 *>
00123 *> \param[in,out] JPVT
00124 *> \verbatim
00125 *>          JPVT is INTEGER array, dimension (N)
00126 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
00127 *>          initial column, otherwise it is a free column.  Before
00128 *>          the QR factorization of A, all initial columns are
00129 *>          permuted to the leading positions; only the remaining
00130 *>          free columns are moved as a result of column pivoting
00131 *>          during the factorization.
00132 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
00133 *>          was the k-th column of A.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] RCOND
00137 *> \verbatim
00138 *>          RCOND is DOUBLE PRECISION
00139 *>          RCOND is used to determine the effective rank of A, which
00140 *>          is defined as the order of the largest leading triangular
00141 *>          submatrix R11 in the QR factorization with pivoting of A,
00142 *>          whose estimated condition number < 1/RCOND.
00143 *> \endverbatim
00144 *>
00145 *> \param[out] RANK
00146 *> \verbatim
00147 *>          RANK is INTEGER
00148 *>          The effective rank of A, i.e., the order of the submatrix
00149 *>          R11.  This is the same as the order of the submatrix T11
00150 *>          in the complete orthogonal factorization of A.
00151 *> \endverbatim
00152 *>
00153 *> \param[out] WORK
00154 *> \verbatim
00155 *>          WORK is COMPLEX*16 array, dimension
00156 *>                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
00157 *> \endverbatim
00158 *>
00159 *> \param[out] RWORK
00160 *> \verbatim
00161 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00162 *> \endverbatim
00163 *>
00164 *> \param[out] INFO
00165 *> \verbatim
00166 *>          INFO is INTEGER
00167 *>          = 0:  successful exit
00168 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00169 *> \endverbatim
00170 *
00171 *  Authors:
00172 *  ========
00173 *
00174 *> \author Univ. of Tennessee 
00175 *> \author Univ. of California Berkeley 
00176 *> \author Univ. of Colorado Denver 
00177 *> \author NAG Ltd. 
00178 *
00179 *> \date November 2011
00180 *
00181 *> \ingroup complex16GEsolve
00182 *
00183 *  =====================================================================
00184       SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00185      $                   WORK, RWORK, INFO )
00186 *
00187 *  -- LAPACK driver routine (version 3.4.0) --
00188 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00189 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00190 *     November 2011
00191 *
00192 *     .. Scalar Arguments ..
00193       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
00194       DOUBLE PRECISION   RCOND
00195 *     ..
00196 *     .. Array Arguments ..
00197       INTEGER            JPVT( * )
00198       DOUBLE PRECISION   RWORK( * )
00199       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00200 *     ..
00201 *
00202 *  =====================================================================
00203 *
00204 *     .. Parameters ..
00205       INTEGER            IMAX, IMIN
00206       PARAMETER          ( IMAX = 1, IMIN = 2 )
00207       DOUBLE PRECISION   ZERO, ONE, DONE, NTDONE
00208       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, DONE = ZERO,
00209      $                   NTDONE = ONE )
00210       COMPLEX*16         CZERO, CONE
00211       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00212      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00213 *     ..
00214 *     .. Local Scalars ..
00215       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
00216       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
00217      $                   SMLNUM
00218       COMPLEX*16         C1, C2, S1, S2, T1, T2
00219 *     ..
00220 *     .. External Subroutines ..
00221       EXTERNAL           XERBLA, ZGEQPF, ZLAIC1, ZLASCL, ZLASET, ZLATZM,
00222      $                   ZTRSM, ZTZRQF, ZUNM2R
00223 *     ..
00224 *     .. External Functions ..
00225       DOUBLE PRECISION   DLAMCH, ZLANGE
00226       EXTERNAL           DLAMCH, ZLANGE
00227 *     ..
00228 *     .. Intrinsic Functions ..
00229       INTRINSIC          ABS, DCONJG, MAX, MIN
00230 *     ..
00231 *     .. Executable Statements ..
00232 *
00233       MN = MIN( M, N )
00234       ISMIN = MN + 1
00235       ISMAX = 2*MN + 1
00236 *
00237 *     Test the input arguments.
00238 *
00239       INFO = 0
00240       IF( M.LT.0 ) THEN
00241          INFO = -1
00242       ELSE IF( N.LT.0 ) THEN
00243          INFO = -2
00244       ELSE IF( NRHS.LT.0 ) THEN
00245          INFO = -3
00246       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00247          INFO = -5
00248       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00249          INFO = -7
00250       END IF
00251 *
00252       IF( INFO.NE.0 ) THEN
00253          CALL XERBLA( 'ZGELSX', -INFO )
00254          RETURN
00255       END IF
00256 *
00257 *     Quick return if possible
00258 *
00259       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00260          RANK = 0
00261          RETURN
00262       END IF
00263 *
00264 *     Get machine parameters
00265 *
00266       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
00267       BIGNUM = ONE / SMLNUM
00268       CALL DLABAD( SMLNUM, BIGNUM )
00269 *
00270 *     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
00271 *
00272       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
00273       IASCL = 0
00274       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00275 *
00276 *        Scale matrix norm up to SMLNUM
00277 *
00278          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00279          IASCL = 1
00280       ELSE IF( ANRM.GT.BIGNUM ) THEN
00281 *
00282 *        Scale matrix norm down to BIGNUM
00283 *
00284          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00285          IASCL = 2
00286       ELSE IF( ANRM.EQ.ZERO ) THEN
00287 *
00288 *        Matrix all zero. Return zero solution.
00289 *
00290          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00291          RANK = 0
00292          GO TO 100
00293       END IF
00294 *
00295       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
00296       IBSCL = 0
00297       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00298 *
00299 *        Scale matrix norm up to SMLNUM
00300 *
00301          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00302          IBSCL = 1
00303       ELSE IF( BNRM.GT.BIGNUM ) THEN
00304 *
00305 *        Scale matrix norm down to BIGNUM
00306 *
00307          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00308          IBSCL = 2
00309       END IF
00310 *
00311 *     Compute QR factorization with column pivoting of A:
00312 *        A * P = Q * R
00313 *
00314       CALL ZGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK,
00315      $             INFO )
00316 *
00317 *     complex workspace MN+N. Real workspace 2*N. Details of Householder
00318 *     rotations stored in WORK(1:MN).
00319 *
00320 *     Determine RANK using incremental condition estimation
00321 *
00322       WORK( ISMIN ) = CONE
00323       WORK( ISMAX ) = CONE
00324       SMAX = ABS( A( 1, 1 ) )
00325       SMIN = SMAX
00326       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
00327          RANK = 0
00328          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00329          GO TO 100
00330       ELSE
00331          RANK = 1
00332       END IF
00333 *
00334    10 CONTINUE
00335       IF( RANK.LT.MN ) THEN
00336          I = RANK + 1
00337          CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
00338      $                A( I, I ), SMINPR, S1, C1 )
00339          CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
00340      $                A( I, I ), SMAXPR, S2, C2 )
00341 *
00342          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
00343             DO 20 I = 1, RANK
00344                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
00345                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
00346    20       CONTINUE
00347             WORK( ISMIN+RANK ) = C1
00348             WORK( ISMAX+RANK ) = C2
00349             SMIN = SMINPR
00350             SMAX = SMAXPR
00351             RANK = RANK + 1
00352             GO TO 10
00353          END IF
00354       END IF
00355 *
00356 *     Logically partition R = [ R11 R12 ]
00357 *                             [  0  R22 ]
00358 *     where R11 = R(1:RANK,1:RANK)
00359 *
00360 *     [R11,R12] = [ T11, 0 ] * Y
00361 *
00362       IF( RANK.LT.N )
00363      $   CALL ZTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
00364 *
00365 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
00366 *
00367 *     B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
00368 *
00369       CALL ZUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
00370      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO )
00371 *
00372 *     workspace NRHS
00373 *
00374 *      B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
00375 *
00376       CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
00377      $            NRHS, CONE, A, LDA, B, LDB )
00378 *
00379       DO 40 I = RANK + 1, N
00380          DO 30 J = 1, NRHS
00381             B( I, J ) = CZERO
00382    30    CONTINUE
00383    40 CONTINUE
00384 *
00385 *     B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
00386 *
00387       IF( RANK.LT.N ) THEN
00388          DO 50 I = 1, RANK
00389             CALL ZLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
00390      $                   DCONJG( WORK( MN+I ) ), B( I, 1 ),
00391      $                   B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) )
00392    50    CONTINUE
00393       END IF
00394 *
00395 *     workspace NRHS
00396 *
00397 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
00398 *
00399       DO 90 J = 1, NRHS
00400          DO 60 I = 1, N
00401             WORK( 2*MN+I ) = NTDONE
00402    60    CONTINUE
00403          DO 80 I = 1, N
00404             IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
00405                IF( JPVT( I ).NE.I ) THEN
00406                   K = I
00407                   T1 = B( K, J )
00408                   T2 = B( JPVT( K ), J )
00409    70             CONTINUE
00410                   B( JPVT( K ), J ) = T1
00411                   WORK( 2*MN+K ) = DONE
00412                   T1 = T2
00413                   K = JPVT( K )
00414                   T2 = B( JPVT( K ), J )
00415                   IF( JPVT( K ).NE.I )
00416      $               GO TO 70
00417                   B( I, J ) = T1
00418                   WORK( 2*MN+K ) = DONE
00419                END IF
00420             END IF
00421    80    CONTINUE
00422    90 CONTINUE
00423 *
00424 *     Undo scaling
00425 *
00426       IF( IASCL.EQ.1 ) THEN
00427          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00428          CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
00429      $                INFO )
00430       ELSE IF( IASCL.EQ.2 ) THEN
00431          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00432          CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
00433      $                INFO )
00434       END IF
00435       IF( IBSCL.EQ.1 ) THEN
00436          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00437       ELSE IF( IBSCL.EQ.2 ) THEN
00438          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00439       END IF
00440 *
00441   100 CONTINUE
00442 *
00443       RETURN
00444 *
00445 *     End of ZGELSX
00446 *
00447       END
 All Files Functions