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