![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGEQPF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGEQPF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqpf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqpf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqpf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, 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 *> This routine is deprecated and has been replaced by routine SGEQP3. 00038 *> 00039 *> SGEQPF computes a QR factorization with column pivoting of a 00040 *> real M-by-N matrix A: A*P = Q*R. 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 REAL 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 triangular matrix R; the elements 00064 *> below the diagonal, together with the array TAU, 00065 *> represent the orthogonal matrix Q as a product of 00066 *> min(m,n) elementary 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(i) .ne. 0, the i-th column of A is permuted 00079 *> to the front of A*P (a leading column); if JPVT(i) = 0, 00080 *> the i-th column of A is a free column. 00081 *> On exit, if JPVT(i) = k, then the i-th column of A*P 00082 *> was the k-th column of A. 00083 *> \endverbatim 00084 *> 00085 *> \param[out] TAU 00086 *> \verbatim 00087 *> TAU is REAL 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 REAL array, dimension (3*N) 00094 *> \endverbatim 00095 *> 00096 *> \param[out] INFO 00097 *> \verbatim 00098 *> INFO is INTEGER 00099 *> = 0: successful exit 00100 *> < 0: if INFO = -i, the i-th argument had an illegal value 00101 *> \endverbatim 00102 * 00103 * Authors: 00104 * ======== 00105 * 00106 *> \author Univ. of Tennessee 00107 *> \author Univ. of California Berkeley 00108 *> \author Univ. of Colorado Denver 00109 *> \author NAG Ltd. 00110 * 00111 *> \date November 2011 00112 * 00113 *> \ingroup realGEcomputational 00114 * 00115 *> \par Further Details: 00116 * ===================== 00117 *> 00118 *> \verbatim 00119 *> 00120 *> The matrix Q is represented as a product of elementary reflectors 00121 *> 00122 *> Q = H(1) H(2) . . . H(n) 00123 *> 00124 *> Each H(i) has the form 00125 *> 00126 *> H = I - tau * v * v**T 00127 *> 00128 *> where tau is a real scalar, and v is a real vector with 00129 *> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). 00130 *> 00131 *> The matrix P is represented in jpvt as follows: If 00132 *> jpvt(j) = i 00133 *> then the jth column of P is the ith canonical unit vector. 00134 *> 00135 *> Partial column norm updating strategy modified by 00136 *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 00137 *> University of Zagreb, Croatia. 00138 *> -- April 2011 -- 00139 *> For more details see LAPACK Working Note 176. 00140 *> \endverbatim 00141 *> 00142 * ===================================================================== 00143 SUBROUTINE SGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO ) 00144 * 00145 * -- LAPACK computational routine (version 3.4.0) -- 00146 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00148 * November 2011 00149 * 00150 * .. Scalar Arguments .. 00151 INTEGER INFO, LDA, M, N 00152 * .. 00153 * .. Array Arguments .. 00154 INTEGER JPVT( * ) 00155 REAL A( LDA, * ), TAU( * ), WORK( * ) 00156 * .. 00157 * 00158 * ===================================================================== 00159 * 00160 * .. Parameters .. 00161 REAL ZERO, ONE 00162 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00163 * .. 00164 * .. Local Scalars .. 00165 INTEGER I, ITEMP, J, MA, MN, PVT 00166 REAL AII, TEMP, TEMP2, TOL3Z 00167 * .. 00168 * .. External Subroutines .. 00169 EXTERNAL SGEQR2, SLARF, SLARFG, SORM2R, SSWAP, XERBLA 00170 * .. 00171 * .. Intrinsic Functions .. 00172 INTRINSIC ABS, MAX, MIN, SQRT 00173 * .. 00174 * .. External Functions .. 00175 INTEGER ISAMAX 00176 REAL SLAMCH, SNRM2 00177 EXTERNAL ISAMAX, SLAMCH, SNRM2 00178 * .. 00179 * .. Executable Statements .. 00180 * 00181 * Test the input arguments 00182 * 00183 INFO = 0 00184 IF( M.LT.0 ) THEN 00185 INFO = -1 00186 ELSE IF( N.LT.0 ) THEN 00187 INFO = -2 00188 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00189 INFO = -4 00190 END IF 00191 IF( INFO.NE.0 ) THEN 00192 CALL XERBLA( 'SGEQPF', -INFO ) 00193 RETURN 00194 END IF 00195 * 00196 MN = MIN( M, N ) 00197 TOL3Z = SQRT(SLAMCH('Epsilon')) 00198 * 00199 * Move initial columns up front 00200 * 00201 ITEMP = 1 00202 DO 10 I = 1, N 00203 IF( JPVT( I ).NE.0 ) THEN 00204 IF( I.NE.ITEMP ) THEN 00205 CALL SSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 ) 00206 JPVT( I ) = JPVT( ITEMP ) 00207 JPVT( ITEMP ) = I 00208 ELSE 00209 JPVT( I ) = I 00210 END IF 00211 ITEMP = ITEMP + 1 00212 ELSE 00213 JPVT( I ) = I 00214 END IF 00215 10 CONTINUE 00216 ITEMP = ITEMP - 1 00217 * 00218 * Compute the QR factorization and update remaining columns 00219 * 00220 IF( ITEMP.GT.0 ) THEN 00221 MA = MIN( ITEMP, M ) 00222 CALL SGEQR2( M, MA, A, LDA, TAU, WORK, INFO ) 00223 IF( MA.LT.N ) THEN 00224 CALL SORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU, 00225 $ A( 1, MA+1 ), LDA, WORK, INFO ) 00226 END IF 00227 END IF 00228 * 00229 IF( ITEMP.LT.MN ) THEN 00230 * 00231 * Initialize partial column norms. The first n elements of 00232 * work store the exact column norms. 00233 * 00234 DO 20 I = ITEMP + 1, N 00235 WORK( I ) = SNRM2( M-ITEMP, A( ITEMP+1, I ), 1 ) 00236 WORK( N+I ) = WORK( I ) 00237 20 CONTINUE 00238 * 00239 * Compute factorization 00240 * 00241 DO 40 I = ITEMP + 1, MN 00242 * 00243 * Determine ith pivot column and swap if necessary 00244 * 00245 PVT = ( I-1 ) + ISAMAX( N-I+1, WORK( I ), 1 ) 00246 * 00247 IF( PVT.NE.I ) THEN 00248 CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) 00249 ITEMP = JPVT( PVT ) 00250 JPVT( PVT ) = JPVT( I ) 00251 JPVT( I ) = ITEMP 00252 WORK( PVT ) = WORK( I ) 00253 WORK( N+PVT ) = WORK( N+I ) 00254 END IF 00255 * 00256 * Generate elementary reflector H(i) 00257 * 00258 IF( I.LT.M ) THEN 00259 CALL SLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) ) 00260 ELSE 00261 CALL SLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) ) 00262 END IF 00263 * 00264 IF( I.LT.N ) THEN 00265 * 00266 * Apply H(i) to A(i:m,i+1:n) from the left 00267 * 00268 AII = A( I, I ) 00269 A( I, I ) = ONE 00270 CALL SLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ), 00271 $ A( I, I+1 ), LDA, WORK( 2*N+1 ) ) 00272 A( I, I ) = AII 00273 END IF 00274 * 00275 * Update partial column norms 00276 * 00277 DO 30 J = I + 1, N 00278 IF( WORK( J ).NE.ZERO ) THEN 00279 * 00280 * NOTE: The following 4 lines follow from the analysis in 00281 * Lapack Working Note 176. 00282 * 00283 TEMP = ABS( A( I, J ) ) / WORK( J ) 00284 TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 00285 TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2 00286 IF( TEMP2 .LE. TOL3Z ) THEN 00287 IF( M-I.GT.0 ) THEN 00288 WORK( J ) = SNRM2( M-I, A( I+1, J ), 1 ) 00289 WORK( N+J ) = WORK( J ) 00290 ELSE 00291 WORK( J ) = ZERO 00292 WORK( N+J ) = ZERO 00293 END IF 00294 ELSE 00295 WORK( J ) = WORK( J )*SQRT( TEMP ) 00296 END IF 00297 END IF 00298 30 CONTINUE 00299 * 00300 40 CONTINUE 00301 END IF 00302 RETURN 00303 * 00304 * End of SGEQPF 00305 * 00306 END