LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgelsy.f
Go to the documentation of this file.
00001 *> \brief <b> ZGELSY 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 ZGELSY + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsy.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsy.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsy.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00022 *                          WORK, LWORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDA, LDB, LWORK, 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 *> ZGELSY computes the minimum-norm solution to a complex linear least
00041 *> squares problem:
00042 *>     minimize || A * X - B ||
00043 *> using a complete orthogonal factorization of A.  A is an M-by-N
00044 *> matrix which may be rank-deficient.
00045 *>
00046 *> Several right hand side vectors b and solution vectors x can be
00047 *> handled in a single call; they are stored as the columns of the
00048 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00049 *> matrix X.
00050 *>
00051 *> The routine first computes a QR factorization with column pivoting:
00052 *>     A * P = Q * [ R11 R12 ]
00053 *>                 [  0  R22 ]
00054 *> with R11 defined as the largest leading submatrix whose estimated
00055 *> condition number is less than 1/RCOND.  The order of R11, RANK,
00056 *> is the effective rank of A.
00057 *>
00058 *> Then, R22 is considered to be negligible, and R12 is annihilated
00059 *> by unitary transformations from the right, arriving at the
00060 *> complete orthogonal factorization:
00061 *>    A * P = Q * [ T11 0 ] * Z
00062 *>                [  0  0 ]
00063 *> The minimum-norm solution is then
00064 *>    X = P * Z**H [ inv(T11)*Q1**H*B ]
00065 *>                 [        0         ]
00066 *> where Q1 consists of the first RANK columns of Q.
00067 *>
00068 *> This routine is basically identical to the original xGELSX except
00069 *> three differences:
00070 *>   o The permutation of matrix B (the right hand side) is faster and
00071 *>     more simple.
00072 *>   o The call to the subroutine xGEQPF has been substituted by the
00073 *>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
00074 *>     version of the QR factorization with column pivoting.
00075 *>   o Matrix B (the right hand side) is updated with Blas-3.
00076 *> \endverbatim
00077 *
00078 *  Arguments:
00079 *  ==========
00080 *
00081 *> \param[in] M
00082 *> \verbatim
00083 *>          M is INTEGER
00084 *>          The number of rows of the matrix A.  M >= 0.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] N
00088 *> \verbatim
00089 *>          N is INTEGER
00090 *>          The number of columns of the matrix A.  N >= 0.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] NRHS
00094 *> \verbatim
00095 *>          NRHS is INTEGER
00096 *>          The number of right hand sides, i.e., the number of
00097 *>          columns of matrices B and X. NRHS >= 0.
00098 *> \endverbatim
00099 *>
00100 *> \param[in,out] A
00101 *> \verbatim
00102 *>          A is COMPLEX*16 array, dimension (LDA,N)
00103 *>          On entry, the M-by-N matrix A.
00104 *>          On exit, A has been overwritten by details of its
00105 *>          complete orthogonal factorization.
00106 *> \endverbatim
00107 *>
00108 *> \param[in] LDA
00109 *> \verbatim
00110 *>          LDA is INTEGER
00111 *>          The leading dimension of the array A.  LDA >= max(1,M).
00112 *> \endverbatim
00113 *>
00114 *> \param[in,out] B
00115 *> \verbatim
00116 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
00117 *>          On entry, the M-by-NRHS right hand side matrix B.
00118 *>          On exit, the N-by-NRHS solution matrix X.
00119 *> \endverbatim
00120 *>
00121 *> \param[in] LDB
00122 *> \verbatim
00123 *>          LDB is INTEGER
00124 *>          The leading dimension of the array B. LDB >= max(1,M,N).
00125 *> \endverbatim
00126 *>
00127 *> \param[in,out] JPVT
00128 *> \verbatim
00129 *>          JPVT is INTEGER array, dimension (N)
00130 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00131 *>          to the front of AP, otherwise column i is a free column.
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 (MAX(1,LWORK))
00156 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] LWORK
00160 *> \verbatim
00161 *>          LWORK is INTEGER
00162 *>          The dimension of the array WORK.
00163 *>          The unblocked strategy requires that:
00164 *>            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
00165 *>          where MN = min(M,N).
00166 *>          The block algorithm requires that:
00167 *>            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
00168 *>          where NB is an upper bound on the blocksize returned
00169 *>          by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
00170 *>          and ZUNMRZ.
00171 *>
00172 *>          If LWORK = -1, then a workspace query is assumed; the routine
00173 *>          only calculates the optimal size of the WORK array, returns
00174 *>          this value as the first entry of the WORK array, and no error
00175 *>          message related to LWORK is issued by XERBLA.
00176 *> \endverbatim
00177 *>
00178 *> \param[out] RWORK
00179 *> \verbatim
00180 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00181 *> \endverbatim
00182 *>
00183 *> \param[out] INFO
00184 *> \verbatim
00185 *>          INFO is INTEGER
00186 *>          = 0: successful exit
00187 *>          < 0: if INFO = -i, the i-th argument had an illegal value
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 complex16GEsolve
00201 *
00202 *> \par Contributors:
00203 *  ==================
00204 *>
00205 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n 
00206 *>    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
00207 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
00208 *>
00209 *  =====================================================================
00210       SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00211      $                   WORK, LWORK, RWORK, 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       DOUBLE PRECISION   RCOND
00221 *     ..
00222 *     .. Array Arguments ..
00223       INTEGER            JPVT( * )
00224       DOUBLE PRECISION   RWORK( * )
00225       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00226 *     ..
00227 *
00228 *  =====================================================================
00229 *
00230 *     .. Parameters ..
00231       INTEGER            IMAX, IMIN
00232       PARAMETER          ( IMAX = 1, IMIN = 2 )
00233       DOUBLE PRECISION   ZERO, ONE
00234       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00235       COMPLEX*16         CZERO, CONE
00236       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00237      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00238 *     ..
00239 *     .. Local Scalars ..
00240       LOGICAL            LQUERY
00241       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
00242      $                   NB, NB1, NB2, NB3, NB4
00243       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
00244      $                   SMLNUM, WSIZE
00245       COMPLEX*16         C1, C2, S1, S2
00246 *     ..
00247 *     .. External Subroutines ..
00248       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZGEQP3, ZLAIC1, ZLASCL,
00249      $                   ZLASET, ZTRSM, ZTZRZF, ZUNMQR, ZUNMRZ
00250 *     ..
00251 *     .. External Functions ..
00252       INTEGER            ILAENV
00253       DOUBLE PRECISION   DLAMCH, ZLANGE
00254       EXTERNAL           ILAENV, DLAMCH, ZLANGE
00255 *     ..
00256 *     .. Intrinsic Functions ..
00257       INTRINSIC          ABS, DBLE, DCMPLX, MAX, MIN
00258 *     ..
00259 *     .. Executable Statements ..
00260 *
00261       MN = MIN( M, N )
00262       ISMIN = MN + 1
00263       ISMAX = 2*MN + 1
00264 *
00265 *     Test the input arguments.
00266 *
00267       INFO = 0
00268       NB1 = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
00269       NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
00270       NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, NRHS, -1 )
00271       NB4 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, NRHS, -1 )
00272       NB = MAX( NB1, NB2, NB3, NB4 )
00273       LWKOPT = MAX( 1, MN+2*N+NB*( N+1 ), 2*MN+NB*NRHS )
00274       WORK( 1 ) = DCMPLX( LWKOPT )
00275       LQUERY = ( LWORK.EQ.-1 )
00276       IF( M.LT.0 ) THEN
00277          INFO = -1
00278       ELSE IF( N.LT.0 ) THEN
00279          INFO = -2
00280       ELSE IF( NRHS.LT.0 ) THEN
00281          INFO = -3
00282       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00283          INFO = -5
00284       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00285          INFO = -7
00286       ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND. .NOT.
00287      $         LQUERY ) THEN
00288          INFO = -12
00289       END IF
00290 *
00291       IF( INFO.NE.0 ) THEN
00292          CALL XERBLA( 'ZGELSY', -INFO )
00293          RETURN
00294       ELSE IF( LQUERY ) THEN
00295          RETURN
00296       END IF
00297 *
00298 *     Quick return if possible
00299 *
00300       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00301          RANK = 0
00302          RETURN
00303       END IF
00304 *
00305 *     Get machine parameters
00306 *
00307       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
00308       BIGNUM = ONE / SMLNUM
00309       CALL DLABAD( SMLNUM, BIGNUM )
00310 *
00311 *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
00312 *
00313       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
00314       IASCL = 0
00315       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00316 *
00317 *        Scale matrix norm up to SMLNUM
00318 *
00319          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00320          IASCL = 1
00321       ELSE IF( ANRM.GT.BIGNUM ) THEN
00322 *
00323 *        Scale matrix norm down to BIGNUM
00324 *
00325          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00326          IASCL = 2
00327       ELSE IF( ANRM.EQ.ZERO ) THEN
00328 *
00329 *        Matrix all zero. Return zero solution.
00330 *
00331          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00332          RANK = 0
00333          GO TO 70
00334       END IF
00335 *
00336       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
00337       IBSCL = 0
00338       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00339 *
00340 *        Scale matrix norm up to SMLNUM
00341 *
00342          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00343          IBSCL = 1
00344       ELSE IF( BNRM.GT.BIGNUM ) THEN
00345 *
00346 *        Scale matrix norm down to BIGNUM
00347 *
00348          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00349          IBSCL = 2
00350       END IF
00351 *
00352 *     Compute QR factorization with column pivoting of A:
00353 *        A * P = Q * R
00354 *
00355       CALL ZGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
00356      $             LWORK-MN, RWORK, INFO )
00357       WSIZE = MN + DBLE( WORK( MN+1 ) )
00358 *
00359 *     complex workspace: MN+NB*(N+1). real workspace 2*N.
00360 *     Details of Householder rotations stored in WORK(1:MN).
00361 *
00362 *     Determine RANK using incremental condition estimation
00363 *
00364       WORK( ISMIN ) = CONE
00365       WORK( ISMAX ) = CONE
00366       SMAX = ABS( A( 1, 1 ) )
00367       SMIN = SMAX
00368       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
00369          RANK = 0
00370          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00371          GO TO 70
00372       ELSE
00373          RANK = 1
00374       END IF
00375 *
00376    10 CONTINUE
00377       IF( RANK.LT.MN ) THEN
00378          I = RANK + 1
00379          CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
00380      $                A( I, I ), SMINPR, S1, C1 )
00381          CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
00382      $                A( I, I ), SMAXPR, S2, C2 )
00383 *
00384          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
00385             DO 20 I = 1, RANK
00386                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
00387                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
00388    20       CONTINUE
00389             WORK( ISMIN+RANK ) = C1
00390             WORK( ISMAX+RANK ) = C2
00391             SMIN = SMINPR
00392             SMAX = SMAXPR
00393             RANK = RANK + 1
00394             GO TO 10
00395          END IF
00396       END IF
00397 *
00398 *     complex workspace: 3*MN.
00399 *
00400 *     Logically partition R = [ R11 R12 ]
00401 *                             [  0  R22 ]
00402 *     where R11 = R(1:RANK,1:RANK)
00403 *
00404 *     [R11,R12] = [ T11, 0 ] * Y
00405 *
00406       IF( RANK.LT.N )
00407      $   CALL ZTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
00408      $                LWORK-2*MN, INFO )
00409 *
00410 *     complex workspace: 2*MN.
00411 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
00412 *
00413 *     B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
00414 *
00415       CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
00416      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
00417       WSIZE = MAX( WSIZE, 2*MN+DBLE( WORK( 2*MN+1 ) ) )
00418 *
00419 *     complex workspace: 2*MN+NB*NRHS.
00420 *
00421 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
00422 *
00423       CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
00424      $            NRHS, CONE, A, LDA, B, LDB )
00425 *
00426       DO 40 J = 1, NRHS
00427          DO 30 I = RANK + 1, N
00428             B( I, J ) = CZERO
00429    30    CONTINUE
00430    40 CONTINUE
00431 *
00432 *     B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
00433 *
00434       IF( RANK.LT.N ) THEN
00435          CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK,
00436      $                N-RANK, A, LDA, WORK( MN+1 ), B, LDB,
00437      $                WORK( 2*MN+1 ), LWORK-2*MN, INFO )
00438       END IF
00439 *
00440 *     complex workspace: 2*MN+NRHS.
00441 *
00442 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
00443 *
00444       DO 60 J = 1, NRHS
00445          DO 50 I = 1, N
00446             WORK( JPVT( I ) ) = B( I, J )
00447    50    CONTINUE
00448          CALL ZCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
00449    60 CONTINUE
00450 *
00451 *     complex workspace: N.
00452 *
00453 *     Undo scaling
00454 *
00455       IF( IASCL.EQ.1 ) THEN
00456          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00457          CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
00458      $                INFO )
00459       ELSE IF( IASCL.EQ.2 ) THEN
00460          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00461          CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
00462      $                INFO )
00463       END IF
00464       IF( IBSCL.EQ.1 ) THEN
00465          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00466       ELSE IF( IBSCL.EQ.2 ) THEN
00467          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00468       END IF
00469 *
00470    70 CONTINUE
00471       WORK( 1 ) = DCMPLX( LWKOPT )
00472 *
00473       RETURN
00474 *
00475 *     End of ZGELSY
00476 *
00477       END
 All Files Functions