![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGEQP3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGEQP3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqp3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqp3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqp3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGEQP3( 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 * REAL A( LDA, * ), TAU( * ), WORK( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SGEQP3 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 REAL 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 REAL 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 REAL 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 realGEcomputational 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 SGEQP3( 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 REAL 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 SGEQRF, SLAQP2, SLAQPS, SORMQR, SSWAP, XERBLA 00180 * .. 00181 * .. External Functions .. 00182 INTEGER ILAENV 00183 REAL SNRM2 00184 EXTERNAL ILAENV, SNRM2 00185 * .. 00186 * .. Intrinsic Functions .. 00187 INTRINSIC INT, MAX, MIN 00188 * .. 00189 * .. Executable Statements .. 00190 * 00191 INFO = 0 00192 LQUERY = ( LWORK.EQ.-1 ) 00193 IF( M.LT.0 ) THEN 00194 INFO = -1 00195 ELSE IF( N.LT.0 ) THEN 00196 INFO = -2 00197 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00198 INFO = -4 00199 END IF 00200 * 00201 IF( INFO.EQ.0 ) THEN 00202 MINMN = MIN( M, N ) 00203 IF( MINMN.EQ.0 ) THEN 00204 IWS = 1 00205 LWKOPT = 1 00206 ELSE 00207 IWS = 3*N + 1 00208 NB = ILAENV( INB, 'SGEQRF', ' ', M, N, -1, -1 ) 00209 LWKOPT = 2*N + ( N + 1 )*NB 00210 END IF 00211 WORK( 1 ) = LWKOPT 00212 * 00213 IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN 00214 INFO = -8 00215 END IF 00216 END IF 00217 * 00218 IF( INFO.NE.0 ) THEN 00219 CALL XERBLA( 'SGEQP3', -INFO ) 00220 RETURN 00221 ELSE IF( LQUERY ) THEN 00222 RETURN 00223 END IF 00224 * 00225 * Quick return if possible. 00226 * 00227 IF( MINMN.EQ.0 ) THEN 00228 RETURN 00229 END IF 00230 * 00231 * Move initial columns up front. 00232 * 00233 NFXD = 1 00234 DO 10 J = 1, N 00235 IF( JPVT( J ).NE.0 ) THEN 00236 IF( J.NE.NFXD ) THEN 00237 CALL SSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 ) 00238 JPVT( J ) = JPVT( NFXD ) 00239 JPVT( NFXD ) = J 00240 ELSE 00241 JPVT( J ) = J 00242 END IF 00243 NFXD = NFXD + 1 00244 ELSE 00245 JPVT( J ) = J 00246 END IF 00247 10 CONTINUE 00248 NFXD = NFXD - 1 00249 * 00250 * Factorize fixed columns 00251 * ======================= 00252 * 00253 * Compute the QR factorization of fixed columns and update 00254 * remaining columns. 00255 * 00256 IF( NFXD.GT.0 ) THEN 00257 NA = MIN( M, NFXD ) 00258 *CC CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO ) 00259 CALL SGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO ) 00260 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 00261 IF( NA.LT.N ) THEN 00262 *CC CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA, 00263 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO ) 00264 CALL SORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU, 00265 $ A( 1, NA+1 ), LDA, WORK, LWORK, INFO ) 00266 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 00267 END IF 00268 END IF 00269 * 00270 * Factorize free columns 00271 * ====================== 00272 * 00273 IF( NFXD.LT.MINMN ) THEN 00274 * 00275 SM = M - NFXD 00276 SN = N - NFXD 00277 SMINMN = MINMN - NFXD 00278 * 00279 * Determine the block size. 00280 * 00281 NB = ILAENV( INB, 'SGEQRF', ' ', SM, SN, -1, -1 ) 00282 NBMIN = 2 00283 NX = 0 00284 * 00285 IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN 00286 * 00287 * Determine when to cross over from blocked to unblocked code. 00288 * 00289 NX = MAX( 0, ILAENV( IXOVER, 'SGEQRF', ' ', SM, SN, -1, 00290 $ -1 ) ) 00291 * 00292 * 00293 IF( NX.LT.SMINMN ) THEN 00294 * 00295 * Determine if workspace is large enough for blocked code. 00296 * 00297 MINWS = 2*SN + ( SN+1 )*NB 00298 IWS = MAX( IWS, MINWS ) 00299 IF( LWORK.LT.MINWS ) THEN 00300 * 00301 * Not enough workspace to use optimal NB: Reduce NB and 00302 * determine the minimum value of NB. 00303 * 00304 NB = ( LWORK-2*SN ) / ( SN+1 ) 00305 NBMIN = MAX( 2, ILAENV( INBMIN, 'SGEQRF', ' ', SM, SN, 00306 $ -1, -1 ) ) 00307 * 00308 * 00309 END IF 00310 END IF 00311 END IF 00312 * 00313 * Initialize partial column norms. The first N elements of work 00314 * store the exact column norms. 00315 * 00316 DO 20 J = NFXD + 1, N 00317 WORK( J ) = SNRM2( SM, A( NFXD+1, J ), 1 ) 00318 WORK( N+J ) = WORK( J ) 00319 20 CONTINUE 00320 * 00321 IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND. 00322 $ ( NX.LT.SMINMN ) ) THEN 00323 * 00324 * Use blocked code initially. 00325 * 00326 J = NFXD + 1 00327 * 00328 * Compute factorization: while loop. 00329 * 00330 * 00331 TOPBMN = MINMN - NX 00332 30 CONTINUE 00333 IF( J.LE.TOPBMN ) THEN 00334 JB = MIN( NB, TOPBMN-J+1 ) 00335 * 00336 * Factorize JB columns among columns J:N. 00337 * 00338 CALL SLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA, 00339 $ JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ), 00340 $ WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 ) 00341 * 00342 J = J + FJB 00343 GO TO 30 00344 END IF 00345 ELSE 00346 J = NFXD + 1 00347 END IF 00348 * 00349 * Use unblocked code to factor the last or only block. 00350 * 00351 * 00352 IF( J.LE.MINMN ) 00353 $ CALL SLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ), 00354 $ TAU( J ), WORK( J ), WORK( N+J ), 00355 $ WORK( 2*N+1 ) ) 00356 * 00357 END IF 00358 * 00359 WORK( 1 ) = IWS 00360 RETURN 00361 * 00362 * End of SGEQP3 00363 * 00364 END