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