LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgeqp3.f
Go to the documentation of this file.
00001 *> \brief \b DGEQP3
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DGEQP3 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqp3.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqp3.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqp3.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, LWORK, M, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       INTEGER            JPVT( * )
00028 *       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DGEQP3 computes a QR factorization with column pivoting of a
00038 *> matrix A:  A*P = Q*R  using Level 3 BLAS.
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] M
00045 *> \verbatim
00046 *>          M is INTEGER
00047 *>          The number of rows of the matrix A. M >= 0.
00048 *> \endverbatim
00049 *>
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The number of columns of the matrix A.  N >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in,out] A
00057 *> \verbatim
00058 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00059 *>          On entry, the M-by-N matrix A.
00060 *>          On exit, the upper triangle of the array contains the
00061 *>          min(M,N)-by-N upper trapezoidal matrix R; the elements below
00062 *>          the diagonal, together with the array TAU, represent the
00063 *>          orthogonal matrix Q as a product of min(M,N) elementary
00064 *>          reflectors.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] LDA
00068 *> \verbatim
00069 *>          LDA is INTEGER
00070 *>          The leading dimension of the array A. LDA >= max(1,M).
00071 *> \endverbatim
00072 *>
00073 *> \param[in,out] JPVT
00074 *> \verbatim
00075 *>          JPVT is INTEGER array, dimension (N)
00076 *>          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
00077 *>          to the front of A*P (a leading column); if JPVT(J)=0,
00078 *>          the J-th column of A is a free column.
00079 *>          On exit, if JPVT(J)=K, then the J-th column of A*P was the
00080 *>          the K-th column of A.
00081 *> \endverbatim
00082 *>
00083 *> \param[out] TAU
00084 *> \verbatim
00085 *>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
00086 *>          The scalar factors of the elementary reflectors.
00087 *> \endverbatim
00088 *>
00089 *> \param[out] WORK
00090 *> \verbatim
00091 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00092 *>          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LWORK
00096 *> \verbatim
00097 *>          LWORK is INTEGER
00098 *>          The dimension of the array WORK. LWORK >= 3*N+1.
00099 *>          For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
00100 *>          is the optimal blocksize.
00101 *>
00102 *>          If LWORK = -1, then a workspace query is assumed; the routine
00103 *>          only calculates the optimal size of the WORK array, returns
00104 *>          this value as the first entry of the WORK array, and no error
00105 *>          message related to LWORK is issued by XERBLA.
00106 *> \endverbatim
00107 *>
00108 *> \param[out] INFO
00109 *> \verbatim
00110 *>          INFO is INTEGER
00111 *>          = 0: successful exit.
00112 *>          < 0: if INFO = -i, the i-th argument had an illegal value.
00113 *> \endverbatim
00114 *
00115 *  Authors:
00116 *  ========
00117 *
00118 *> \author Univ. of Tennessee 
00119 *> \author Univ. of California Berkeley 
00120 *> \author Univ. of Colorado Denver 
00121 *> \author NAG Ltd. 
00122 *
00123 *> \date November 2011
00124 *
00125 *> \ingroup doubleGEcomputational
00126 *
00127 *> \par Further Details:
00128 *  =====================
00129 *>
00130 *> \verbatim
00131 *>
00132 *>  The matrix Q is represented as a product of elementary reflectors
00133 *>
00134 *>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
00135 *>
00136 *>  Each H(i) has the form
00137 *>
00138 *>     H(i) = I - tau * v * v**T
00139 *>
00140 *>  where tau is a real/complex scalar, and v is a real/complex vector
00141 *>  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
00142 *>  A(i+1:m,i), and tau in TAU(i).
00143 *> \endverbatim
00144 *
00145 *> \par Contributors:
00146 *  ==================
00147 *>
00148 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00149 *>    X. Sun, Computer Science Dept., Duke University, USA
00150 *>
00151 *  =====================================================================
00152       SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
00153 *
00154 *  -- LAPACK computational routine (version 3.4.0) --
00155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00157 *     November 2011
00158 *
00159 *     .. Scalar Arguments ..
00160       INTEGER            INFO, LDA, LWORK, M, N
00161 *     ..
00162 *     .. Array Arguments ..
00163       INTEGER            JPVT( * )
00164       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
00165 *     ..
00166 *
00167 *  =====================================================================
00168 *
00169 *     .. Parameters ..
00170       INTEGER            INB, INBMIN, IXOVER
00171       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
00172 *     ..
00173 *     .. Local Scalars ..
00174       LOGICAL            LQUERY
00175       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
00176      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
00177 *     ..
00178 *     .. External Subroutines ..
00179       EXTERNAL           DGEQRF, DLAQP2, DLAQPS, DORMQR, DSWAP, XERBLA
00180 *     ..
00181 *     .. External Functions ..
00182       INTEGER            ILAENV
00183       DOUBLE PRECISION   DNRM2
00184       EXTERNAL           ILAENV, DNRM2
00185 *     ..
00186 *     .. Intrinsic Functions ..
00187       INTRINSIC          INT, MAX, MIN
00188 *     ..
00189 *     .. Executable Statements ..
00190 *
00191 *     Test input arguments
00192 *  ====================
00193 *
00194       INFO = 0
00195       LQUERY = ( LWORK.EQ.-1 )
00196       IF( M.LT.0 ) THEN
00197          INFO = -1
00198       ELSE IF( N.LT.0 ) THEN
00199          INFO = -2
00200       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00201          INFO = -4
00202       END IF
00203 *
00204       IF( INFO.EQ.0 ) THEN
00205          MINMN = MIN( M, N )
00206          IF( MINMN.EQ.0 ) THEN
00207             IWS = 1
00208             LWKOPT = 1
00209          ELSE
00210             IWS = 3*N + 1
00211             NB = ILAENV( INB, 'DGEQRF', ' ', M, N, -1, -1 )
00212             LWKOPT = 2*N + ( N + 1 )*NB
00213          END IF
00214          WORK( 1 ) = LWKOPT
00215 *
00216          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
00217             INFO = -8
00218          END IF
00219       END IF
00220 *
00221       IF( INFO.NE.0 ) THEN
00222          CALL XERBLA( 'DGEQP3', -INFO )
00223          RETURN
00224       ELSE IF( LQUERY ) THEN
00225          RETURN
00226       END IF
00227 *
00228 *     Quick return if possible.
00229 *
00230       IF( MINMN.EQ.0 ) THEN
00231          RETURN
00232       END IF
00233 *
00234 *     Move initial columns up front.
00235 *
00236       NFXD = 1
00237       DO 10 J = 1, N
00238          IF( JPVT( J ).NE.0 ) THEN
00239             IF( J.NE.NFXD ) THEN
00240                CALL DSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
00241                JPVT( J ) = JPVT( NFXD )
00242                JPVT( NFXD ) = J
00243             ELSE
00244                JPVT( J ) = J
00245             END IF
00246             NFXD = NFXD + 1
00247          ELSE
00248             JPVT( J ) = J
00249          END IF
00250    10 CONTINUE
00251       NFXD = NFXD - 1
00252 *
00253 *     Factorize fixed columns
00254 *  =======================
00255 *
00256 *     Compute the QR factorization of fixed columns and update
00257 *     remaining columns.
00258 *
00259       IF( NFXD.GT.0 ) THEN
00260          NA = MIN( M, NFXD )
00261 *CC      CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
00262          CALL DGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
00263          IWS = MAX( IWS, INT( WORK( 1 ) ) )
00264          IF( NA.LT.N ) THEN
00265 *CC         CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
00266 *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
00267             CALL DORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
00268      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
00269             IWS = MAX( IWS, INT( WORK( 1 ) ) )
00270          END IF
00271       END IF
00272 *
00273 *     Factorize free columns
00274 *  ======================
00275 *
00276       IF( NFXD.LT.MINMN ) THEN
00277 *
00278          SM = M - NFXD
00279          SN = N - NFXD
00280          SMINMN = MINMN - NFXD
00281 *
00282 *        Determine the block size.
00283 *
00284          NB = ILAENV( INB, 'DGEQRF', ' ', SM, SN, -1, -1 )
00285          NBMIN = 2
00286          NX = 0
00287 *
00288          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
00289 *
00290 *           Determine when to cross over from blocked to unblocked code.
00291 *
00292             NX = MAX( 0, ILAENV( IXOVER, 'DGEQRF', ' ', SM, SN, -1,
00293      $           -1 ) )
00294 *
00295 *
00296             IF( NX.LT.SMINMN ) THEN
00297 *
00298 *              Determine if workspace is large enough for blocked code.
00299 *
00300                MINWS = 2*SN + ( SN+1 )*NB
00301                IWS = MAX( IWS, MINWS )
00302                IF( LWORK.LT.MINWS ) THEN
00303 *
00304 *                 Not enough workspace to use optimal NB: Reduce NB and
00305 *                 determine the minimum value of NB.
00306 *
00307                   NB = ( LWORK-2*SN ) / ( SN+1 )
00308                   NBMIN = MAX( 2, ILAENV( INBMIN, 'DGEQRF', ' ', SM, SN,
00309      $                    -1, -1 ) )
00310 *
00311 *
00312                END IF
00313             END IF
00314          END IF
00315 *
00316 *        Initialize partial column norms. The first N elements of work
00317 *        store the exact column norms.
00318 *
00319          DO 20 J = NFXD + 1, N
00320             WORK( J ) = DNRM2( SM, A( NFXD+1, J ), 1 )
00321             WORK( N+J ) = WORK( J )
00322    20    CONTINUE
00323 *
00324          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
00325      $       ( NX.LT.SMINMN ) ) THEN
00326 *
00327 *           Use blocked code initially.
00328 *
00329             J = NFXD + 1
00330 *
00331 *           Compute factorization: while loop.
00332 *
00333 *
00334             TOPBMN = MINMN - NX
00335    30       CONTINUE
00336             IF( J.LE.TOPBMN ) THEN
00337                JB = MIN( NB, TOPBMN-J+1 )
00338 *
00339 *              Factorize JB columns among columns J:N.
00340 *
00341                CALL DLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
00342      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
00343      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
00344 *
00345                J = J + FJB
00346                GO TO 30
00347             END IF
00348          ELSE
00349             J = NFXD + 1
00350          END IF
00351 *
00352 *        Use unblocked code to factor the last or only block.
00353 *
00354 *
00355          IF( J.LE.MINMN )
00356      $      CALL DLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
00357      $                   TAU( J ), WORK( J ), WORK( N+J ),
00358      $                   WORK( 2*N+1 ) )
00359 *
00360       END IF
00361 *
00362       WORK( 1 ) = IWS
00363       RETURN
00364 *
00365 *     End of DGEQP3
00366 *
00367       END
 All Files Functions