![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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