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