1 SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, 2 $ INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER INFO, LDA, LWORK, M, N 11* .. 12* .. Array Arguments .. 13 INTEGER JPVT( * ) 14 REAL RWORK( * ) 15 COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* CGEQP3 computes a QR factorization with column pivoting of a 22* matrix A: A*P = Q*R using Level 3 BLAS. 23* 24* Arguments 25* ========= 26* 27* M (input) INTEGER 28* The number of rows of the matrix A. M >= 0. 29* 30* N (input) INTEGER 31* The number of columns of the matrix A. N >= 0. 32* 33* A (input/output) COMPLEX array, dimension (LDA,N) 34* On entry, the M-by-N matrix A. 35* On exit, the upper triangle of the array contains the 36* min(M,N)-by-N upper trapezoidal matrix R; the elements below 37* the diagonal, together with the array TAU, represent the 38* unitary matrix Q as a product of min(M,N) elementary 39* reflectors. 40* 41* LDA (input) INTEGER 42* The leading dimension of the array A. LDA >= max(1,M). 43* 44* JPVT (input/output) INTEGER array, dimension (N) 45* On entry, if JPVT(J).ne.0, the J-th column of A is permuted 46* to the front of A*P (a leading column); if JPVT(J)=0, 47* the J-th column of A is a free column. 48* On exit, if JPVT(J)=K, then the J-th column of A*P was the 49* the K-th column of A. 50* 51* TAU (output) COMPLEX array, dimension (min(M,N)) 52* The scalar factors of the elementary reflectors. 53* 54* WORK (workspace/output) COMPLEX array, dimension (LWORK) 55* On exit, if INFO=0, WORK(1) returns the optimal LWORK. 56* 57* LWORK (input) INTEGER 58* The dimension of the array WORK. LWORK >= N+1. 59* For optimal performance LWORK >= ( N+1 )*NB, where NB 60* is the optimal blocksize. 61* 62* If LWORK = -1, then a workspace query is assumed; the routine 63* only calculates the optimal size of the WORK array, returns 64* this value as the first entry of the WORK array, and no error 65* message related to LWORK is issued by XERBLA. 66* 67* RWORK (workspace) REAL array, dimension (2*N) 68* 69* INFO (output) INTEGER 70* = 0: successful exit. 71* < 0: if INFO = -i, the i-th argument had an illegal value. 72* 73* Further Details 74* =============== 75* 76* The matrix Q is represented as a product of elementary reflectors 77* 78* Q = H(1) H(2) . . . H(k), where k = min(m,n). 79* 80* Each H(i) has the form 81* 82* H(i) = I - tau * v * v' 83* 84* where tau is a real/complex scalar, and v is a real/complex vector 85* with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in 86* A(i+1:m,i), and tau in TAU(i). 87* 88* Based on contributions by 89* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 90* X. Sun, Computer Science Dept., Duke University, USA 91* 92* ===================================================================== 93* 94* .. Parameters .. 95 INTEGER INB, INBMIN, IXOVER 96 PARAMETER ( INB = 1, INBMIN = 2, IXOVER = 3 ) 97* .. 98* .. Local Scalars .. 99 LOGICAL LQUERY 100 INTEGER FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB, 101 $ NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN 102* .. 103* .. External Subroutines .. 104 EXTERNAL CGEQRF, CLAQP2, CLAQPS, CSWAP, CUNMQR, XERBLA 105* .. 106* .. External Functions .. 107 INTEGER ILAENV 108 REAL SCNRM2 109 EXTERNAL ILAENV, SCNRM2 110* .. 111* .. Intrinsic Functions .. 112 INTRINSIC INT, MAX, MIN 113* .. 114* .. Executable Statements .. 115* 116 IWS = N + 1 117 MINMN = MIN( M, N ) 118* 119* Test input arguments 120* ==================== 121* 122 INFO = 0 123 NB = ILAENV( INB, 'CGEQRF', ' ', M, N, -1, -1 ) 124 LWKOPT = ( N+1 )*NB 125 WORK( 1 ) = LWKOPT 126 LQUERY = ( LWORK.EQ.-1 ) 127 IF( M.LT.0 ) THEN 128 INFO = -1 129 ELSE IF( N.LT.0 ) THEN 130 INFO = -2 131 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 132 INFO = -4 133 ELSE IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN 134 INFO = -8 135 END IF 136 IF( INFO.NE.0 ) THEN 137 CALL XERBLA( 'CGEQP3', -INFO ) 138 RETURN 139 ELSE IF( LQUERY ) THEN 140 RETURN 141 END IF 142* 143* Quick return if possible. 144* 145 IF( MINMN.EQ.0 ) THEN 146 WORK( 1 ) = 1 147 RETURN 148 END IF 149* 150* Move initial columns up front. 151* 152 NFXD = 1 153 DO 10 J = 1, N 154 IF( JPVT( J ).NE.0 ) THEN 155 IF( J.NE.NFXD ) THEN 156 CALL CSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 ) 157 JPVT( J ) = JPVT( NFXD ) 158 JPVT( NFXD ) = J 159 ELSE 160 JPVT( J ) = J 161 END IF 162 NFXD = NFXD + 1 163 ELSE 164 JPVT( J ) = J 165 END IF 166 10 CONTINUE 167 NFXD = NFXD - 1 168* 169* Factorize fixed columns 170* ======================= 171* 172* Compute the QR factorization of fixed columns and update 173* remaining columns. 174* 175 IF( NFXD.GT.0 ) THEN 176 NA = MIN( M, NFXD ) 177*CC CALL CGEQR2( M, NA, A, LDA, TAU, WORK, INFO ) 178 CALL CGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO ) 179 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 180 IF( NA.LT.N ) THEN 181*CC CALL CUNM2R( 'Left', 'Conjugate Transpose', M, N-NA, 182*CC $ NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK, 183*CC $ INFO ) 184 CALL CUNMQR( 'Left', 'Conjugate Transpose', M, N-NA, NA, A, 185 $ LDA, TAU, A( 1, NA+1 ), LDA, WORK, LWORK, 186 $ INFO ) 187 IWS = MAX( IWS, INT( WORK( 1 ) ) ) 188 END IF 189 END IF 190* 191* Factorize free columns 192* ====================== 193* 194 IF( NFXD.LT.MINMN ) THEN 195* 196 SM = M - NFXD 197 SN = N - NFXD 198 SMINMN = MINMN - NFXD 199* 200* Determine the block size. 201* 202 NB = ILAENV( INB, 'CGEQRF', ' ', SM, SN, -1, -1 ) 203 NBMIN = 2 204 NX = 0 205* 206 IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN 207* 208* Determine when to cross over from blocked to unblocked code. 209* 210 NX = MAX( 0, ILAENV( IXOVER, 'CGEQRF', ' ', SM, SN, -1, 211 $ -1 ) ) 212* 213* 214 IF( NX.LT.SMINMN ) THEN 215* 216* Determine if workspace is large enough for blocked code. 217* 218 MINWS = ( SN+1 )*NB 219 IWS = MAX( IWS, MINWS ) 220 IF( LWORK.LT.MINWS ) THEN 221* 222* Not enough workspace to use optimal NB: Reduce NB and 223* determine the minimum value of NB. 224* 225 NB = LWORK / ( SN+1 ) 226 NBMIN = MAX( 2, ILAENV( INBMIN, 'CGEQRF', ' ', SM, SN, 227 $ -1, -1 ) ) 228* 229* 230 END IF 231 END IF 232 END IF 233* 234* Initialize partial column norms. The first N elements of work 235* store the exact column norms. 236* 237 DO 20 J = NFXD + 1, N 238 RWORK( J ) = SCNRM2( SM, A( NFXD+1, J ), 1 ) 239 RWORK( N+J ) = RWORK( J ) 240 20 CONTINUE 241* 242 IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND. 243 $ ( NX.LT.SMINMN ) ) THEN 244* 245* Use blocked code initially. 246* 247 J = NFXD + 1 248* 249* Compute factorization: while loop. 250* 251* 252 TOPBMN = MINMN - NX 253 30 CONTINUE 254 IF( J.LE.TOPBMN ) THEN 255 JB = MIN( NB, TOPBMN-J+1 ) 256* 257* Factorize JB columns among columns J:N. 258* 259 CALL CLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA, 260 $ JPVT( J ), TAU( J ), RWORK( J ), 261 $ RWORK( N+J ), WORK( 1 ), WORK( JB+1 ), 262 $ N-J+1 ) 263* 264 J = J + FJB 265 GO TO 30 266 END IF 267 ELSE 268 J = NFXD + 1 269 END IF 270* 271* Use unblocked code to factor the last or only block. 272* 273* 274 IF( J.LE.MINMN ) 275 $ CALL CLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ), 276 $ TAU( J ), RWORK( J ), RWORK( N+J ), WORK( 1 ) ) 277* 278 END IF 279* 280 WORK( 1 ) = IWS 281 RETURN 282* 283* End of CGEQP3 284* 285 END 286