 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LWORK, M, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER JPVT( * ) 00014 REAL RWORK( * ) 00015 COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGEQP3 computes a QR factorization with column pivoting of a 00022 * matrix A: A*P = Q*R using Level 3 BLAS. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * M (input) INTEGER 00028 * The number of rows of the matrix A. M >= 0. 00029 * 00030 * N (input) INTEGER 00031 * The number of columns of the matrix A. N >= 0. 00032 * 00033 * A (input/output) COMPLEX array, dimension (LDA,N) 00034 * On entry, the M-by-N matrix A. 00035 * On exit, the upper triangle of the array contains the 00036 * min(M,N)-by-N upper trapezoidal matrix R; the elements below 00037 * the diagonal, together with the array TAU, represent the 00038 * unitary matrix Q as a product of min(M,N) elementary 00039 * reflectors. 00040 * 00041 * LDA (input) INTEGER 00042 * The leading dimension of the array A. LDA >= max(1,M). 00043 * 00044 * JPVT (input/output) INTEGER array, dimension (N) 00045 * On entry, if JPVT(J).ne.0, the J-th column of A is permuted 00046 * to the front of A*P (a leading column); if JPVT(J)=0, 00047 * the J-th column of A is a free column. 00048 * On exit, if JPVT(J)=K, then the J-th column of A*P was the 00049 * the K-th column of A. 00050 * 00051 * TAU (output) COMPLEX array, dimension (min(M,N)) 00052 * The scalar factors of the elementary reflectors. 00053 * 00054 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00055 * On exit, if INFO=0, WORK(1) returns the optimal LWORK. 00056 * 00057 * LWORK (input) INTEGER 00058 * The dimension of the array WORK. LWORK >= N+1. 00059 * For optimal performance LWORK >= ( N+1 )*NB, where NB 00060 * is the optimal blocksize. 00061 * 00062 * If LWORK = -1, then a workspace query is assumed; the routine 00063 * only calculates the optimal size of the WORK array, returns 00064 * this value as the first entry of the WORK array, and no error 00065 * message related to LWORK is issued by XERBLA. 00066 * 00067 * RWORK (workspace) REAL array, dimension (2*N) 00068 * 00069 * INFO (output) INTEGER 00070 * = 0: successful exit. 00071 * < 0: if INFO = -i, the i-th argument had an illegal value. 00072 * 00073 * Further Details 00074 * =============== 00075 * 00076 * The matrix Q is represented as a product of elementary reflectors 00077 * 00078 * Q = H(1) H(2) . . . H(k), where k = min(m,n). 00079 * 00080 * Each H(i) has the form 00081 * 00082 * H(i) = I - tau * v * v**H 00083 * 00084 * where tau is a real/complex scalar, and v is a real/complex vector 00085 * with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in 00086 * A(i+1:m,i), and tau in TAU(i). 00087 * 00088 * Based on contributions by 00089 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 00090 * X. Sun, Computer Science Dept., Duke University, USA 00091 * 00092 * ===================================================================== 00093 * 00094 * .. Parameters .. 00095 INTEGER INB, INBMIN, IXOVER 00096 PARAMETER ( INB = 1, INBMIN = 2, IXOVER = 3 ) 00097 * .. 00098 * .. Local Scalars .. 00099 LOGICAL LQUERY 00100 INTEGER FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB, 00101 $ NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN 00102 * .. 00103 * .. External Subroutines .. 00104 EXTERNAL CGEQRF, CLAQP2, CLAQPS, CSWAP, CUNMQR, XERBLA 00105 * .. 00106 * .. External Functions .. 00107 INTEGER ILAENV 00108 REAL SCNRM2 00109 EXTERNAL ILAENV, SCNRM2 00110 * .. 00111 * .. Intrinsic Functions .. 00112 INTRINSIC INT, MAX, MIN 00113 * .. 00114 * .. Executable Statements .. 00115 * 00116 * Test input arguments 00117 * ==================== 00118 * 00119 INFO = 0 00120 LQUERY = ( LWORK.EQ.-1 ) 00121 IF( M.LT.0 ) THEN 00122 INFO = -1 00123 ELSE IF( N.LT.0 ) THEN 00124 INFO = -2 00125 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00126 INFO = -4 00127 END IF 00128 * 00129 IF( INFO.EQ.0 ) THEN 00130 MINMN = MIN( M, N ) 00131 IF( MINMN.EQ.0 ) THEN 00132 IWS = 1 00133 LWKOPT = 1 00134 ELSE 00135 IWS = N + 1 00136 NB = ILAENV( INB, 'CGEQRF', ' ', M, N, -1, -1 ) 00137 LWKOPT = ( N + 1 )*NB 00138 END IF 00139 WORK( 1 ) = LWKOPT 00140 * 00141 IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN 00142 INFO = -8 00143 END IF 00144 END IF 00145 * 00146 IF( INFO.NE.0 ) THEN 00147 CALL XERBLA( 'CGEQP3', -INFO ) 00148 RETURN 00149 ELSE IF( LQUERY ) THEN 00150 RETURN 00151 END IF 00152 * 00153 * Quick return if possible. 00154 * 00155 IF( MINMN.EQ.0 ) THEN 00156 RETURN 00157 END IF 00158 * 00159 * Move initial columns up front. 00160 * 00161 NFXD = 1 00162 DO 10 J = 1, N 00163 IF( JPVT( J ).NE.0 ) THEN 00164 IF( J.NE.NFXD ) THEN 00165 CALL CSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 ) 00166 JPVT( J ) = JPVT( NFXD ) 00167 JPVT( NFXD ) = J 00168 ELSE 00169 JPVT( J ) = J 00170 END IF 00171 NFXD = NFXD + 1 00172 ELSE 00173 JPVT( J ) = J 00174 END IF 00175 10 CONTINUE 00176 NFXD = NFXD - 1 00177 * 00178 * Factorize fixed columns 00179 * ======================= 00180 * 00181 * Compute the QR factorization of fixed columns and update 00182 * remaining columns. 00183 * 00184 IF( NFXD.GT.0 ) THEN 00185 NA = MIN( M, NFXD ) 00186 *CC CALL CGEQR2( M, NA, A, LDA, TAU, WORK, INFO ) 00187 CALL CGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO ) 00188 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 00189 IF( NA.LT.N ) THEN 00190 *CC CALL CUNM2R( 'Left', 'Conjugate Transpose', M, N-NA, 00191 *CC $ NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK, 00192 *CC $ INFO ) 00193 CALL CUNMQR( 'Left', 'Conjugate Transpose', M, N-NA, NA, A, 00194 $ LDA, TAU, A( 1, NA+1 ), LDA, WORK, LWORK, 00195 $ INFO ) 00196 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 00197 END IF 00198 END IF 00199 * 00200 * Factorize free columns 00201 * ====================== 00202 * 00203 IF( NFXD.LT.MINMN ) THEN 00204 * 00205 SM = M - NFXD 00206 SN = N - NFXD 00207 SMINMN = MINMN - NFXD 00208 * 00209 * Determine the block size. 00210 * 00211 NB = ILAENV( INB, 'CGEQRF', ' ', SM, SN, -1, -1 ) 00212 NBMIN = 2 00213 NX = 0 00214 * 00215 IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN 00216 * 00217 * Determine when to cross over from blocked to unblocked code. 00218 * 00219 NX = MAX( 0, ILAENV( IXOVER, 'CGEQRF', ' ', SM, SN, -1, 00220 $ -1 ) ) 00221 * 00222 * 00223 IF( NX.LT.SMINMN ) THEN 00224 * 00225 * Determine if workspace is large enough for blocked code. 00226 * 00227 MINWS = ( SN+1 )*NB 00228 IWS = MAX( IWS, MINWS ) 00229 IF( LWORK.LT.MINWS ) THEN 00230 * 00231 * Not enough workspace to use optimal NB: Reduce NB and 00232 * determine the minimum value of NB. 00233 * 00234 NB = LWORK / ( SN+1 ) 00235 NBMIN = MAX( 2, ILAENV( INBMIN, 'CGEQRF', ' ', SM, SN, 00236 $ -1, -1 ) ) 00237 * 00238 * 00239 END IF 00240 END IF 00241 END IF 00242 * 00243 * Initialize partial column norms. The first N elements of work 00244 * store the exact column norms. 00245 * 00246 DO 20 J = NFXD + 1, N 00247 RWORK( J ) = SCNRM2( SM, A( NFXD+1, J ), 1 ) 00248 RWORK( N+J ) = RWORK( J ) 00249 20 CONTINUE 00250 * 00251 IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND. 00252 $ ( NX.LT.SMINMN ) ) THEN 00253 * 00254 * Use blocked code initially. 00255 * 00256 J = NFXD + 1 00257 * 00258 * Compute factorization: while loop. 00259 * 00260 * 00261 TOPBMN = MINMN - NX 00262 30 CONTINUE 00263 IF( J.LE.TOPBMN ) THEN 00264 JB = MIN( NB, TOPBMN-J+1 ) 00265 * 00266 * Factorize JB columns among columns J:N. 00267 * 00268 CALL CLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA, 00269 $ JPVT( J ), TAU( J ), RWORK( J ), 00270 $ RWORK( N+J ), WORK( 1 ), WORK( JB+1 ), 00271 $ N-J+1 ) 00272 * 00273 J = J + FJB 00274 GO TO 30 00275 END IF 00276 ELSE 00277 J = NFXD + 1 00278 END IF 00279 * 00280 * Use unblocked code to factor the last or only block. 00281 * 00282 * 00283 IF( J.LE.MINMN ) 00284 $ CALL CLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ), 00285 $ TAU( J ), RWORK( J ), RWORK( N+J ), WORK( 1 ) ) 00286 * 00287 END IF 00288 * 00289 WORK( 1 ) = IWS 00290 RETURN 00291 * 00292 * End of CGEQP3 00293 * 00294 END
1.7.3 
