|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 00002 $ VN2, AUXV, F, LDF ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER KB, LDA, LDF, M, N, NB, OFFSET 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER JPVT( * ) 00014 DOUBLE PRECISION A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), 00015 $ VN1( * ), VN2( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLAQPS computes a step of QR factorization with column pivoting 00022 * of a real M-by-N matrix A by using Blas-3. It tries to factorize 00023 * NB columns from A starting from the row OFFSET+1, and updates all 00024 * of the matrix with Blas-3 xGEMM. 00025 * 00026 * In some cases, due to catastrophic cancellations, it cannot 00027 * factorize NB columns. Hence, the actual number of factorized 00028 * columns is returned in KB. 00029 * 00030 * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * M (input) INTEGER 00036 * The number of rows of the matrix A. M >= 0. 00037 * 00038 * N (input) INTEGER 00039 * The number of columns of the matrix A. N >= 0 00040 * 00041 * OFFSET (input) INTEGER 00042 * The number of rows of A that have been factorized in 00043 * previous steps. 00044 * 00045 * NB (input) INTEGER 00046 * The number of columns to factorize. 00047 * 00048 * KB (output) INTEGER 00049 * The number of columns actually factorized. 00050 * 00051 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00052 * On entry, the M-by-N matrix A. 00053 * On exit, block A(OFFSET+1:M,1:KB) is the triangular 00054 * factor obtained and block A(1:OFFSET,1:N) has been 00055 * accordingly pivoted, but no factorized. 00056 * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has 00057 * been updated. 00058 * 00059 * LDA (input) INTEGER 00060 * The leading dimension of the array A. LDA >= max(1,M). 00061 * 00062 * JPVT (input/output) INTEGER array, dimension (N) 00063 * JPVT(I) = K <==> Column K of the full matrix A has been 00064 * permuted into position I in AP. 00065 * 00066 * TAU (output) DOUBLE PRECISION array, dimension (KB) 00067 * The scalar factors of the elementary reflectors. 00068 * 00069 * VN1 (input/output) DOUBLE PRECISION array, dimension (N) 00070 * The vector with the partial column norms. 00071 * 00072 * VN2 (input/output) DOUBLE PRECISION array, dimension (N) 00073 * The vector with the exact column norms. 00074 * 00075 * AUXV (input/output) DOUBLE PRECISION array, dimension (NB) 00076 * Auxiliar vector. 00077 * 00078 * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB) 00079 * Matrix F' = L*Y'*A. 00080 * 00081 * LDF (input) INTEGER 00082 * The leading dimension of the array F. LDF >= max(1,N). 00083 * 00084 * Further Details 00085 * =============== 00086 * 00087 * Based on contributions by 00088 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 00089 * X. Sun, Computer Science Dept., Duke University, USA 00090 * 00091 * Partial column norm updating strategy modified by 00092 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 00093 * University of Zagreb, Croatia. 00094 * June 2010 00095 * For more details see LAPACK Working Note 176. 00096 * ===================================================================== 00097 * 00098 * .. Parameters .. 00099 DOUBLE PRECISION ZERO, ONE 00100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00101 * .. 00102 * .. Local Scalars .. 00103 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK 00104 DOUBLE PRECISION AKK, TEMP, TEMP2, TOL3Z 00105 * .. 00106 * .. External Subroutines .. 00107 EXTERNAL DGEMM, DGEMV, DLARFG, DSWAP 00108 * .. 00109 * .. Intrinsic Functions .. 00110 INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT 00111 * .. 00112 * .. External Functions .. 00113 INTEGER IDAMAX 00114 DOUBLE PRECISION DLAMCH, DNRM2 00115 EXTERNAL IDAMAX, DLAMCH, DNRM2 00116 * .. 00117 * .. Executable Statements .. 00118 * 00119 LASTRK = MIN( M, N+OFFSET ) 00120 LSTICC = 0 00121 K = 0 00122 TOL3Z = SQRT(DLAMCH('Epsilon')) 00123 * 00124 * Beginning of while loop. 00125 * 00126 10 CONTINUE 00127 IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN 00128 K = K + 1 00129 RK = OFFSET + K 00130 * 00131 * Determine ith pivot column and swap if necessary 00132 * 00133 PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 ) 00134 IF( PVT.NE.K ) THEN 00135 CALL DSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) 00136 CALL DSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) 00137 ITEMP = JPVT( PVT ) 00138 JPVT( PVT ) = JPVT( K ) 00139 JPVT( K ) = ITEMP 00140 VN1( PVT ) = VN1( K ) 00141 VN2( PVT ) = VN2( K ) 00142 END IF 00143 * 00144 * Apply previous Householder reflectors to column K: 00145 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. 00146 * 00147 IF( K.GT.1 ) THEN 00148 CALL DGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ), 00149 $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 ) 00150 END IF 00151 * 00152 * Generate elementary reflector H(k). 00153 * 00154 IF( RK.LT.M ) THEN 00155 CALL DLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) 00156 ELSE 00157 CALL DLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) 00158 END IF 00159 * 00160 AKK = A( RK, K ) 00161 A( RK, K ) = ONE 00162 * 00163 * Compute Kth column of F: 00164 * 00165 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). 00166 * 00167 IF( K.LT.N ) THEN 00168 CALL DGEMV( 'Transpose', M-RK+1, N-K, TAU( K ), 00169 $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO, 00170 $ F( K+1, K ), 1 ) 00171 END IF 00172 * 00173 * Padding F(1:K,K) with zeros. 00174 * 00175 DO 20 J = 1, K 00176 F( J, K ) = ZERO 00177 20 CONTINUE 00178 * 00179 * Incremental updating of F: 00180 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' 00181 * *A(RK:M,K). 00182 * 00183 IF( K.GT.1 ) THEN 00184 CALL DGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ), 00185 $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 ) 00186 * 00187 CALL DGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF, 00188 $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 ) 00189 END IF 00190 * 00191 * Update the current row of A: 00192 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. 00193 * 00194 IF( K.LT.N ) THEN 00195 CALL DGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF, 00196 $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA ) 00197 END IF 00198 * 00199 * Update partial column norms. 00200 * 00201 IF( RK.LT.LASTRK ) THEN 00202 DO 30 J = K + 1, N 00203 IF( VN1( J ).NE.ZERO ) THEN 00204 * 00205 * NOTE: The following 4 lines follow from the analysis in 00206 * Lapack Working Note 176. 00207 * 00208 TEMP = ABS( A( RK, J ) ) / VN1( J ) 00209 TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 00210 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 00211 IF( TEMP2 .LE. TOL3Z ) THEN 00212 VN2( J ) = DBLE( LSTICC ) 00213 LSTICC = J 00214 ELSE 00215 VN1( J ) = VN1( J )*SQRT( TEMP ) 00216 END IF 00217 END IF 00218 30 CONTINUE 00219 END IF 00220 * 00221 A( RK, K ) = AKK 00222 * 00223 * End of while loop. 00224 * 00225 GO TO 10 00226 END IF 00227 KB = K 00228 RK = OFFSET + KB 00229 * 00230 * Apply the block reflector to the rest of the matrix: 00231 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - 00232 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. 00233 * 00234 IF( KB.LT.MIN( N, M-OFFSET ) ) THEN 00235 CALL DGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE, 00236 $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE, 00237 $ A( RK+1, KB+1 ), LDA ) 00238 END IF 00239 * 00240 * Recomputation of difficult columns. 00241 * 00242 40 CONTINUE 00243 IF( LSTICC.GT.0 ) THEN 00244 ITEMP = NINT( VN2( LSTICC ) ) 00245 VN1( LSTICC ) = DNRM2( M-RK, A( RK+1, LSTICC ), 1 ) 00246 * 00247 * NOTE: The computation of VN1( LSTICC ) relies on the fact that 00248 * SNRM2 does not fail on vectors with norm below the value of 00249 * SQRT(DLAMCH('S')) 00250 * 00251 VN2( LSTICC ) = VN1( LSTICC ) 00252 LSTICC = ITEMP 00253 GO TO 40 00254 END IF 00255 * 00256 RETURN 00257 * 00258 * End of DLAQPS 00259 * 00260 END
1.7.3 
