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