LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgeqpf.f
Go to the documentation of this file.
00001 *> \brief \b ZGEQPF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZGEQPF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqpf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqpf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqpf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, M, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       INTEGER            JPVT( * )
00028 *       DOUBLE PRECISION   RWORK( * )
00029 *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> This routine is deprecated and has been replaced by routine ZGEQP3.
00039 *>
00040 *> ZGEQPF computes a QR factorization with column pivoting of a
00041 *> complex M-by-N matrix A: A*P = Q*R.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] M
00048 *> \verbatim
00049 *>          M is INTEGER
00050 *>          The number of rows of the matrix A. M >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>          The number of columns of the matrix A. N >= 0
00057 *> \endverbatim
00058 *>
00059 *> \param[in,out] A
00060 *> \verbatim
00061 *>          A is COMPLEX*16 array, dimension (LDA,N)
00062 *>          On entry, the M-by-N matrix A.
00063 *>          On exit, the upper triangle of the array contains the
00064 *>          min(M,N)-by-N upper triangular matrix R; the elements
00065 *>          below the diagonal, together with the array TAU,
00066 *>          represent the unitary matrix Q as a product of
00067 *>          min(m,n) elementary reflectors.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] LDA
00071 *> \verbatim
00072 *>          LDA is INTEGER
00073 *>          The leading dimension of the array A. LDA >= max(1,M).
00074 *> \endverbatim
00075 *>
00076 *> \param[in,out] JPVT
00077 *> \verbatim
00078 *>          JPVT is INTEGER array, dimension (N)
00079 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00080 *>          to the front of A*P (a leading column); if JPVT(i) = 0,
00081 *>          the i-th column of A is a free column.
00082 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
00083 *>          was the k-th column of A.
00084 *> \endverbatim
00085 *>
00086 *> \param[out] TAU
00087 *> \verbatim
00088 *>          TAU is COMPLEX*16 array, dimension (min(M,N))
00089 *>          The scalar factors of the elementary reflectors.
00090 *> \endverbatim
00091 *>
00092 *> \param[out] WORK
00093 *> \verbatim
00094 *>          WORK is COMPLEX*16 array, dimension (N)
00095 *> \endverbatim
00096 *>
00097 *> \param[out] RWORK
00098 *> \verbatim
00099 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00100 *> \endverbatim
00101 *>
00102 *> \param[out] INFO
00103 *> \verbatim
00104 *>          INFO is INTEGER
00105 *>          = 0:  successful exit
00106 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00107 *> \endverbatim
00108 *
00109 *  Authors:
00110 *  ========
00111 *
00112 *> \author Univ. of Tennessee 
00113 *> \author Univ. of California Berkeley 
00114 *> \author Univ. of Colorado Denver 
00115 *> \author NAG Ltd. 
00116 *
00117 *> \date November 2011
00118 *
00119 *> \ingroup complex16GEcomputational
00120 *
00121 *> \par Further Details:
00122 *  =====================
00123 *>
00124 *> \verbatim
00125 *>
00126 *>  The matrix Q is represented as a product of elementary reflectors
00127 *>
00128 *>     Q = H(1) H(2) . . . H(n)
00129 *>
00130 *>  Each H(i) has the form
00131 *>
00132 *>     H = I - tau * v * v**H
00133 *>
00134 *>  where tau is a complex scalar, and v is a complex vector with
00135 *>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
00136 *>
00137 *>  The matrix P is represented in jpvt as follows: If
00138 *>     jpvt(j) = i
00139 *>  then the jth column of P is the ith canonical unit vector.
00140 *>
00141 *>  Partial column norm updating strategy modified by
00142 *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00143 *>    University of Zagreb, Croatia.
00144 *>  -- April 2011                                                      --
00145 *>  For more details see LAPACK Working Note 176.
00146 *> \endverbatim
00147 *>
00148 *  =====================================================================
00149       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
00150 *
00151 *  -- LAPACK computational routine (version 3.4.0) --
00152 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00153 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00154 *     November 2011
00155 *
00156 *     .. Scalar Arguments ..
00157       INTEGER            INFO, LDA, M, N
00158 *     ..
00159 *     .. Array Arguments ..
00160       INTEGER            JPVT( * )
00161       DOUBLE PRECISION   RWORK( * )
00162       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
00163 *     ..
00164 *
00165 *  =====================================================================
00166 *
00167 *     .. Parameters ..
00168       DOUBLE PRECISION   ZERO, ONE
00169       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00170 *     ..
00171 *     .. Local Scalars ..
00172       INTEGER            I, ITEMP, J, MA, MN, PVT
00173       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
00174       COMPLEX*16         AII
00175 *     ..
00176 *     .. External Subroutines ..
00177       EXTERNAL           XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
00178 *     ..
00179 *     .. Intrinsic Functions ..
00180       INTRINSIC          ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
00181 *     ..
00182 *     .. External Functions ..
00183       INTEGER            IDAMAX
00184       DOUBLE PRECISION   DLAMCH, DZNRM2
00185       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
00186 *     ..
00187 *     .. Executable Statements ..
00188 *
00189 *     Test the input arguments
00190 *
00191       INFO = 0
00192       IF( M.LT.0 ) THEN
00193          INFO = -1
00194       ELSE IF( N.LT.0 ) THEN
00195          INFO = -2
00196       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00197          INFO = -4
00198       END IF
00199       IF( INFO.NE.0 ) THEN
00200          CALL XERBLA( 'ZGEQPF', -INFO )
00201          RETURN
00202       END IF
00203 *
00204       MN = MIN( M, N )
00205       TOL3Z = SQRT(DLAMCH('Epsilon'))
00206 *
00207 *     Move initial columns up front
00208 *
00209       ITEMP = 1
00210       DO 10 I = 1, N
00211          IF( JPVT( I ).NE.0 ) THEN
00212             IF( I.NE.ITEMP ) THEN
00213                CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
00214                JPVT( I ) = JPVT( ITEMP )
00215                JPVT( ITEMP ) = I
00216             ELSE
00217                JPVT( I ) = I
00218             END IF
00219             ITEMP = ITEMP + 1
00220          ELSE
00221             JPVT( I ) = I
00222          END IF
00223    10 CONTINUE
00224       ITEMP = ITEMP - 1
00225 *
00226 *     Compute the QR factorization and update remaining columns
00227 *
00228       IF( ITEMP.GT.0 ) THEN
00229          MA = MIN( ITEMP, M )
00230          CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
00231          IF( MA.LT.N ) THEN
00232             CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
00233      $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
00234          END IF
00235       END IF
00236 *
00237       IF( ITEMP.LT.MN ) THEN
00238 *
00239 *        Initialize partial column norms. The first n elements of
00240 *        work store the exact column norms.
00241 *
00242          DO 20 I = ITEMP + 1, N
00243             RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
00244             RWORK( N+I ) = RWORK( I )
00245    20    CONTINUE
00246 *
00247 *        Compute factorization
00248 *
00249          DO 40 I = ITEMP + 1, MN
00250 *
00251 *           Determine ith pivot column and swap if necessary
00252 *
00253             PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
00254 *
00255             IF( PVT.NE.I ) THEN
00256                CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00257                ITEMP = JPVT( PVT )
00258                JPVT( PVT ) = JPVT( I )
00259                JPVT( I ) = ITEMP
00260                RWORK( PVT ) = RWORK( I )
00261                RWORK( N+PVT ) = RWORK( N+I )
00262             END IF
00263 *
00264 *           Generate elementary reflector H(i)
00265 *
00266             AII = A( I, I )
00267             CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
00268      $                   TAU( I ) )
00269             A( I, I ) = AII
00270 *
00271             IF( I.LT.N ) THEN
00272 *
00273 *              Apply H(i) to A(i:m,i+1:n) from the left
00274 *
00275                AII = A( I, I )
00276                A( I, I ) = DCMPLX( ONE )
00277                CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
00278      $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
00279                A( I, I ) = AII
00280             END IF
00281 *
00282 *           Update partial column norms
00283 *
00284             DO 30 J = I + 1, N
00285                IF( RWORK( J ).NE.ZERO ) THEN
00286 *
00287 *                 NOTE: The following 4 lines follow from the analysis in
00288 *                 Lapack Working Note 176.
00289 *                 
00290                   TEMP = ABS( A( I, J ) ) / RWORK( J )
00291                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
00292                   TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
00293                   IF( TEMP2 .LE. TOL3Z ) THEN 
00294                      IF( M-I.GT.0 ) THEN
00295                         RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
00296                         RWORK( N+J ) = RWORK( J )
00297                      ELSE
00298                         RWORK( J ) = ZERO
00299                         RWORK( N+J ) = ZERO
00300                      END IF
00301                   ELSE
00302                      RWORK( J ) = RWORK( J )*SQRT( TEMP )
00303                   END IF
00304                END IF
00305    30       CONTINUE
00306 *
00307    40    CONTINUE
00308       END IF
00309       RETURN
00310 *
00311 *     End of ZGEQPF
00312 *
00313       END
 All Files Functions