![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLAQPS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLAQPS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqps.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqps.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqps.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 00022 * VN2, AUXV, F, LDF ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER KB, LDA, LDF, M, N, NB, OFFSET 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER JPVT( * ) 00029 * DOUBLE PRECISION VN1( * ), VN2( * ) 00030 * COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZLAQPS computes a step of QR factorization with column pivoting 00040 *> of a complex M-by-N matrix A by using Blas-3. It tries to factorize 00041 *> NB columns from A starting from the row OFFSET+1, and updates all 00042 *> of the matrix with Blas-3 xGEMM. 00043 *> 00044 *> In some cases, due to catastrophic cancellations, it cannot 00045 *> factorize NB columns. Hence, the actual number of factorized 00046 *> columns is returned in KB. 00047 *> 00048 *> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. 00049 *> \endverbatim 00050 * 00051 * Arguments: 00052 * ========== 00053 * 00054 *> \param[in] M 00055 *> \verbatim 00056 *> M is INTEGER 00057 *> The number of rows of the matrix A. M >= 0. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] N 00061 *> \verbatim 00062 *> N is INTEGER 00063 *> The number of columns of the matrix A. N >= 0 00064 *> \endverbatim 00065 *> 00066 *> \param[in] OFFSET 00067 *> \verbatim 00068 *> OFFSET is INTEGER 00069 *> The number of rows of A that have been factorized in 00070 *> previous steps. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] NB 00074 *> \verbatim 00075 *> NB is INTEGER 00076 *> The number of columns to factorize. 00077 *> \endverbatim 00078 *> 00079 *> \param[out] KB 00080 *> \verbatim 00081 *> KB is INTEGER 00082 *> The number of columns actually factorized. 00083 *> \endverbatim 00084 *> 00085 *> \param[in,out] A 00086 *> \verbatim 00087 *> A is COMPLEX*16 array, dimension (LDA,N) 00088 *> On entry, the M-by-N matrix A. 00089 *> On exit, block A(OFFSET+1:M,1:KB) is the triangular 00090 *> factor obtained and block A(1:OFFSET,1:N) has been 00091 *> accordingly pivoted, but no factorized. 00092 *> The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has 00093 *> been updated. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDA 00097 *> \verbatim 00098 *> LDA is INTEGER 00099 *> The leading dimension of the array A. LDA >= max(1,M). 00100 *> \endverbatim 00101 *> 00102 *> \param[in,out] JPVT 00103 *> \verbatim 00104 *> JPVT is INTEGER array, dimension (N) 00105 *> JPVT(I) = K <==> Column K of the full matrix A has been 00106 *> permuted into position I in AP. 00107 *> \endverbatim 00108 *> 00109 *> \param[out] TAU 00110 *> \verbatim 00111 *> TAU is COMPLEX*16 array, dimension (KB) 00112 *> The scalar factors of the elementary reflectors. 00113 *> \endverbatim 00114 *> 00115 *> \param[in,out] VN1 00116 *> \verbatim 00117 *> VN1 is DOUBLE PRECISION array, dimension (N) 00118 *> The vector with the partial column norms. 00119 *> \endverbatim 00120 *> 00121 *> \param[in,out] VN2 00122 *> \verbatim 00123 *> VN2 is DOUBLE PRECISION array, dimension (N) 00124 *> The vector with the exact column norms. 00125 *> \endverbatim 00126 *> 00127 *> \param[in,out] AUXV 00128 *> \verbatim 00129 *> AUXV is COMPLEX*16 array, dimension (NB) 00130 *> Auxiliar vector. 00131 *> \endverbatim 00132 *> 00133 *> \param[in,out] F 00134 *> \verbatim 00135 *> F is COMPLEX*16 array, dimension (LDF,NB) 00136 *> Matrix F**H = L * Y**H * A. 00137 *> \endverbatim 00138 *> 00139 *> \param[in] LDF 00140 *> \verbatim 00141 *> LDF is INTEGER 00142 *> The leading dimension of the array F. LDF >= max(1,N). 00143 *> \endverbatim 00144 * 00145 * Authors: 00146 * ======== 00147 * 00148 *> \author Univ. of Tennessee 00149 *> \author Univ. of California Berkeley 00150 *> \author Univ. of Colorado Denver 00151 *> \author NAG Ltd. 00152 * 00153 *> \date November 2011 00154 * 00155 *> \ingroup complex16OTHERauxiliary 00156 * 00157 *> \par Contributors: 00158 * ================== 00159 *> 00160 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 00161 *> X. Sun, Computer Science Dept., Duke University, USA 00162 *> \n 00163 *> Partial column norm updating strategy modified on April 2011 00164 *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 00165 *> University of Zagreb, Croatia. 00166 * 00167 *> \par References: 00168 * ================ 00169 *> 00170 *> LAPACK Working Note 176 00171 * 00172 *> \htmlonly 00173 *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> 00174 *> \endhtmlonly 00175 * 00176 * ===================================================================== 00177 SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 00178 $ VN2, AUXV, F, LDF ) 00179 * 00180 * -- LAPACK auxiliary routine (version 3.4.0) -- 00181 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00183 * November 2011 00184 * 00185 * .. Scalar Arguments .. 00186 INTEGER KB, LDA, LDF, M, N, NB, OFFSET 00187 * .. 00188 * .. Array Arguments .. 00189 INTEGER JPVT( * ) 00190 DOUBLE PRECISION VN1( * ), VN2( * ) 00191 COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) 00192 * .. 00193 * 00194 * ===================================================================== 00195 * 00196 * .. Parameters .. 00197 DOUBLE PRECISION ZERO, ONE 00198 COMPLEX*16 CZERO, CONE 00199 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 00200 $ CZERO = ( 0.0D+0, 0.0D+0 ), 00201 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00202 * .. 00203 * .. Local Scalars .. 00204 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK 00205 DOUBLE PRECISION TEMP, TEMP2, TOL3Z 00206 COMPLEX*16 AKK 00207 * .. 00208 * .. External Subroutines .. 00209 EXTERNAL ZGEMM, ZGEMV, ZLARFG, ZSWAP 00210 * .. 00211 * .. Intrinsic Functions .. 00212 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, NINT, SQRT 00213 * .. 00214 * .. External Functions .. 00215 INTEGER IDAMAX 00216 DOUBLE PRECISION DLAMCH, DZNRM2 00217 EXTERNAL IDAMAX, DLAMCH, DZNRM2 00218 * .. 00219 * .. Executable Statements .. 00220 * 00221 LASTRK = MIN( M, N+OFFSET ) 00222 LSTICC = 0 00223 K = 0 00224 TOL3Z = SQRT(DLAMCH('Epsilon')) 00225 * 00226 * Beginning of while loop. 00227 * 00228 10 CONTINUE 00229 IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN 00230 K = K + 1 00231 RK = OFFSET + K 00232 * 00233 * Determine ith pivot column and swap if necessary 00234 * 00235 PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 ) 00236 IF( PVT.NE.K ) THEN 00237 CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) 00238 CALL ZSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) 00239 ITEMP = JPVT( PVT ) 00240 JPVT( PVT ) = JPVT( K ) 00241 JPVT( K ) = ITEMP 00242 VN1( PVT ) = VN1( K ) 00243 VN2( PVT ) = VN2( K ) 00244 END IF 00245 * 00246 * Apply previous Householder reflectors to column K: 00247 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**H. 00248 * 00249 IF( K.GT.1 ) THEN 00250 DO 20 J = 1, K - 1 00251 F( K, J ) = DCONJG( F( K, J ) ) 00252 20 CONTINUE 00253 CALL ZGEMV( 'No transpose', M-RK+1, K-1, -CONE, A( RK, 1 ), 00254 $ LDA, F( K, 1 ), LDF, CONE, A( RK, K ), 1 ) 00255 DO 30 J = 1, K - 1 00256 F( K, J ) = DCONJG( F( K, J ) ) 00257 30 CONTINUE 00258 END IF 00259 * 00260 * Generate elementary reflector H(k). 00261 * 00262 IF( RK.LT.M ) THEN 00263 CALL ZLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) 00264 ELSE 00265 CALL ZLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) 00266 END IF 00267 * 00268 AKK = A( RK, K ) 00269 A( RK, K ) = CONE 00270 * 00271 * Compute Kth column of F: 00272 * 00273 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K). 00274 * 00275 IF( K.LT.N ) THEN 00276 CALL ZGEMV( 'Conjugate transpose', M-RK+1, N-K, TAU( K ), 00277 $ A( RK, K+1 ), LDA, A( RK, K ), 1, CZERO, 00278 $ F( K+1, K ), 1 ) 00279 END IF 00280 * 00281 * Padding F(1:K,K) with zeros. 00282 * 00283 DO 40 J = 1, K 00284 F( J, K ) = CZERO 00285 40 CONTINUE 00286 * 00287 * Incremental updating of F: 00288 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H 00289 * *A(RK:M,K). 00290 * 00291 IF( K.GT.1 ) THEN 00292 CALL ZGEMV( 'Conjugate transpose', M-RK+1, K-1, -TAU( K ), 00293 $ A( RK, 1 ), LDA, A( RK, K ), 1, CZERO, 00294 $ AUXV( 1 ), 1 ) 00295 * 00296 CALL ZGEMV( 'No transpose', N, K-1, CONE, F( 1, 1 ), LDF, 00297 $ AUXV( 1 ), 1, CONE, F( 1, K ), 1 ) 00298 END IF 00299 * 00300 * Update the current row of A: 00301 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H. 00302 * 00303 IF( K.LT.N ) THEN 00304 CALL ZGEMM( 'No transpose', 'Conjugate transpose', 1, N-K, 00305 $ K, -CONE, A( RK, 1 ), LDA, F( K+1, 1 ), LDF, 00306 $ CONE, A( RK, K+1 ), LDA ) 00307 END IF 00308 * 00309 * Update partial column norms. 00310 * 00311 IF( RK.LT.LASTRK ) THEN 00312 DO 50 J = K + 1, N 00313 IF( VN1( J ).NE.ZERO ) THEN 00314 * 00315 * NOTE: The following 4 lines follow from the analysis in 00316 * Lapack Working Note 176. 00317 * 00318 TEMP = ABS( A( RK, J ) ) / VN1( J ) 00319 TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 00320 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 00321 IF( TEMP2 .LE. TOL3Z ) THEN 00322 VN2( J ) = DBLE( LSTICC ) 00323 LSTICC = J 00324 ELSE 00325 VN1( J ) = VN1( J )*SQRT( TEMP ) 00326 END IF 00327 END IF 00328 50 CONTINUE 00329 END IF 00330 * 00331 A( RK, K ) = AKK 00332 * 00333 * End of while loop. 00334 * 00335 GO TO 10 00336 END IF 00337 KB = K 00338 RK = OFFSET + KB 00339 * 00340 * Apply the block reflector to the rest of the matrix: 00341 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - 00342 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H. 00343 * 00344 IF( KB.LT.MIN( N, M-OFFSET ) ) THEN 00345 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-RK, N-KB, 00346 $ KB, -CONE, A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, 00347 $ CONE, A( RK+1, KB+1 ), LDA ) 00348 END IF 00349 * 00350 * Recomputation of difficult columns. 00351 * 00352 60 CONTINUE 00353 IF( LSTICC.GT.0 ) THEN 00354 ITEMP = NINT( VN2( LSTICC ) ) 00355 VN1( LSTICC ) = DZNRM2( M-RK, A( RK+1, LSTICC ), 1 ) 00356 * 00357 * NOTE: The computation of VN1( LSTICC ) relies on the fact that 00358 * SNRM2 does not fail on vectors with norm below the value of 00359 * SQRT(DLAMCH('S')) 00360 * 00361 VN2( LSTICC ) = VN1( LSTICC ) 00362 LSTICC = ITEMP 00363 GO TO 60 00364 END IF 00365 * 00366 RETURN 00367 * 00368 * End of ZLAQPS 00369 * 00370 END