[22] | 1 | *> \brief \b SLAQPS |
---|
| 2 | * |
---|
| 3 | * =========== DOCUMENTATION =========== |
---|
| 4 | * |
---|
| 5 | * Online html documentation available at |
---|
| 6 | * http://www.netlib.org/lapack/explore-html/ |
---|
| 7 | * |
---|
| 8 | *> \htmlonly |
---|
| 9 | *> Download SLAQPS + dependencies |
---|
| 10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqps.f"> |
---|
| 11 | *> [TGZ]</a> |
---|
| 12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqps.f"> |
---|
| 13 | *> [ZIP]</a> |
---|
| 14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqps.f"> |
---|
| 15 | *> [TXT]</a> |
---|
| 16 | *> \endhtmlonly |
---|
| 17 | * |
---|
| 18 | * Definition: |
---|
| 19 | * =========== |
---|
| 20 | * |
---|
| 21 | * SUBROUTINE SLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, |
---|
| 22 | * VN2, AUXV, F, LDF ) |
---|
| 23 | * |
---|
| 24 | * .. Scalar Arguments .. |
---|
| 25 | * INTEGER KB, LDA, LDF, M, N, NB, OFFSET |
---|
| 26 | * .. |
---|
| 27 | * .. Array Arguments .. |
---|
| 28 | * INTEGER JPVT( * ) |
---|
| 29 | * REAL A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), |
---|
| 30 | * $ VN1( * ), VN2( * ) |
---|
| 31 | * .. |
---|
| 32 | * |
---|
| 33 | * |
---|
| 34 | *> \par Purpose: |
---|
| 35 | * ============= |
---|
| 36 | *> |
---|
| 37 | *> \verbatim |
---|
| 38 | *> |
---|
| 39 | *> SLAQPS computes a step of QR factorization with column pivoting |
---|
| 40 | *> of a real M-by-N matrix A by using Blas-3. It tries to factorize |
---|
| 41 | *> NB columns from A starting from the row OFFSET+1, and updates all |
---|
| 42 | *> of the matrix with Blas-3 xGEMM. |
---|
| 43 | *> |
---|
| 44 | *> In some cases, due to catastrophic cancellations, it cannot |
---|
| 45 | *> factorize NB columns. Hence, the actual number of factorized |
---|
| 46 | *> columns is returned in KB. |
---|
| 47 | *> |
---|
| 48 | *> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. |
---|
| 49 | *> \endverbatim |
---|
| 50 | * |
---|
| 51 | * Arguments: |
---|
| 52 | * ========== |
---|
| 53 | * |
---|
| 54 | *> \param[in] M |
---|
| 55 | *> \verbatim |
---|
| 56 | *> M is INTEGER |
---|
| 57 | *> The number of rows of the matrix A. M >= 0. |
---|
| 58 | *> \endverbatim |
---|
| 59 | *> |
---|
| 60 | *> \param[in] N |
---|
| 61 | *> \verbatim |
---|
| 62 | *> N is INTEGER |
---|
| 63 | *> The number of columns of the matrix A. N >= 0 |
---|
| 64 | *> \endverbatim |
---|
| 65 | *> |
---|
| 66 | *> \param[in] OFFSET |
---|
| 67 | *> \verbatim |
---|
| 68 | *> OFFSET is INTEGER |
---|
| 69 | *> The number of rows of A that have been factorized in |
---|
| 70 | *> previous steps. |
---|
| 71 | *> \endverbatim |
---|
| 72 | *> |
---|
| 73 | *> \param[in] NB |
---|
| 74 | *> \verbatim |
---|
| 75 | *> NB is INTEGER |
---|
| 76 | *> The number of columns to factorize. |
---|
| 77 | *> \endverbatim |
---|
| 78 | *> |
---|
| 79 | *> \param[out] KB |
---|
| 80 | *> \verbatim |
---|
| 81 | *> KB is INTEGER |
---|
| 82 | *> The number of columns actually factorized. |
---|
| 83 | *> \endverbatim |
---|
| 84 | *> |
---|
| 85 | *> \param[in,out] A |
---|
| 86 | *> \verbatim |
---|
| 87 | *> A is REAL array, dimension (LDA,N) |
---|
| 88 | *> On entry, the M-by-N matrix A. |
---|
| 89 | *> On exit, block A(OFFSET+1:M,1:KB) is the triangular |
---|
| 90 | *> factor obtained and block A(1:OFFSET,1:N) has been |
---|
| 91 | *> accordingly pivoted, but no factorized. |
---|
| 92 | *> The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has |
---|
| 93 | *> been updated. |
---|
| 94 | *> \endverbatim |
---|
| 95 | *> |
---|
| 96 | *> \param[in] LDA |
---|
| 97 | *> \verbatim |
---|
| 98 | *> LDA is INTEGER |
---|
| 99 | *> The leading dimension of the array A. LDA >= max(1,M). |
---|
| 100 | *> \endverbatim |
---|
| 101 | *> |
---|
| 102 | *> \param[in,out] JPVT |
---|
| 103 | *> \verbatim |
---|
| 104 | *> JPVT is INTEGER array, dimension (N) |
---|
| 105 | *> JPVT(I) = K <==> Column K of the full matrix A has been |
---|
| 106 | *> permuted into position I in AP. |
---|
| 107 | *> \endverbatim |
---|
| 108 | *> |
---|
| 109 | *> \param[out] TAU |
---|
| 110 | *> \verbatim |
---|
| 111 | *> TAU is REAL array, dimension (KB) |
---|
| 112 | *> The scalar factors of the elementary reflectors. |
---|
| 113 | *> \endverbatim |
---|
| 114 | *> |
---|
| 115 | *> \param[in,out] VN1 |
---|
| 116 | *> \verbatim |
---|
| 117 | *> VN1 is REAL array, dimension (N) |
---|
| 118 | *> The vector with the partial column norms. |
---|
| 119 | *> \endverbatim |
---|
| 120 | *> |
---|
| 121 | *> \param[in,out] VN2 |
---|
| 122 | *> \verbatim |
---|
| 123 | *> VN2 is REAL array, dimension (N) |
---|
| 124 | *> The vector with the exact column norms. |
---|
| 125 | *> \endverbatim |
---|
| 126 | *> |
---|
| 127 | *> \param[in,out] AUXV |
---|
| 128 | *> \verbatim |
---|
| 129 | *> AUXV is REAL array, dimension (NB) |
---|
| 130 | *> Auxiliar vector. |
---|
| 131 | *> \endverbatim |
---|
| 132 | *> |
---|
| 133 | *> \param[in,out] F |
---|
| 134 | *> \verbatim |
---|
| 135 | *> F is REAL array, dimension (LDF,NB) |
---|
| 136 | *> Matrix F**T = L*Y**T*A. |
---|
| 137 | *> \endverbatim |
---|
| 138 | *> |
---|
| 139 | *> \param[in] LDF |
---|
| 140 | *> \verbatim |
---|
| 141 | *> LDF is INTEGER |
---|
| 142 | *> The leading dimension of the array F. LDF >= max(1,N). |
---|
| 143 | *> \endverbatim |
---|
| 144 | * |
---|
| 145 | * Authors: |
---|
| 146 | * ======== |
---|
| 147 | * |
---|
| 148 | *> \author Univ. of Tennessee |
---|
| 149 | *> \author Univ. of California Berkeley |
---|
| 150 | *> \author Univ. of Colorado Denver |
---|
| 151 | *> \author NAG Ltd. |
---|
| 152 | * |
---|
| 153 | *> \date November 2011 |
---|
| 154 | * |
---|
| 155 | *> \ingroup realOTHERauxiliary |
---|
| 156 | * |
---|
| 157 | *> \par Contributors: |
---|
| 158 | * ================== |
---|
| 159 | *> |
---|
| 160 | *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain |
---|
| 161 | *> X. Sun, Computer Science Dept., Duke University, USA |
---|
| 162 | *> |
---|
| 163 | *> \n |
---|
| 164 | *> Partial column norm updating strategy modified on April 2011 |
---|
| 165 | *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, |
---|
| 166 | *> University of Zagreb, Croatia. |
---|
| 167 | * |
---|
| 168 | *> \par References: |
---|
| 169 | * ================ |
---|
| 170 | *> |
---|
| 171 | *> LAPACK Working Note 176 |
---|
| 172 | * |
---|
| 173 | *> \htmlonly |
---|
| 174 | *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> |
---|
| 175 | *> \endhtmlonly |
---|
| 176 | * |
---|
| 177 | * ===================================================================== |
---|
| 178 | SUBROUTINE SLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, |
---|
| 179 | $ VN2, AUXV, F, LDF ) |
---|
| 180 | * |
---|
| 181 | * -- LAPACK auxiliary routine (version 3.4.0) -- |
---|
| 182 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
| 183 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
| 184 | * November 2011 |
---|
| 185 | * |
---|
| 186 | * .. Scalar Arguments .. |
---|
| 187 | INTEGER KB, LDA, LDF, M, N, NB, OFFSET |
---|
| 188 | * .. |
---|
| 189 | * .. Array Arguments .. |
---|
| 190 | INTEGER JPVT( * ) |
---|
| 191 | REAL A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), |
---|
| 192 | $ VN1( * ), VN2( * ) |
---|
| 193 | * .. |
---|
| 194 | * |
---|
| 195 | * ===================================================================== |
---|
| 196 | * |
---|
| 197 | * .. Parameters .. |
---|
| 198 | REAL ZERO, ONE |
---|
| 199 | PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) |
---|
| 200 | * .. |
---|
| 201 | * .. Local Scalars .. |
---|
| 202 | INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK |
---|
| 203 | REAL AKK, TEMP, TEMP2, TOL3Z |
---|
| 204 | * .. |
---|
| 205 | * .. External Subroutines .. |
---|
| 206 | EXTERNAL SGEMM, SGEMV, SLARFG, SSWAP |
---|
| 207 | * .. |
---|
| 208 | * .. Intrinsic Functions .. |
---|
| 209 | INTRINSIC ABS, MAX, MIN, NINT, REAL, SQRT |
---|
| 210 | * .. |
---|
| 211 | * .. External Functions .. |
---|
| 212 | INTEGER ISAMAX |
---|
| 213 | REAL SLAMCH, SNRM2 |
---|
| 214 | EXTERNAL ISAMAX, SLAMCH, SNRM2 |
---|
| 215 | * .. |
---|
| 216 | * .. Executable Statements .. |
---|
| 217 | * |
---|
| 218 | LASTRK = MIN( M, N+OFFSET ) |
---|
| 219 | LSTICC = 0 |
---|
| 220 | K = 0 |
---|
| 221 | TOL3Z = SQRT(SLAMCH('Epsilon')) |
---|
| 222 | * |
---|
| 223 | * Beginning of while loop. |
---|
| 224 | * |
---|
| 225 | 10 CONTINUE |
---|
| 226 | IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN |
---|
| 227 | K = K + 1 |
---|
| 228 | RK = OFFSET + K |
---|
| 229 | * |
---|
| 230 | * Determine ith pivot column and swap if necessary |
---|
| 231 | * |
---|
| 232 | PVT = ( K-1 ) + ISAMAX( N-K+1, VN1( K ), 1 ) |
---|
| 233 | IF( PVT.NE.K ) THEN |
---|
| 234 | CALL SSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) |
---|
| 235 | CALL SSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) |
---|
| 236 | ITEMP = JPVT( PVT ) |
---|
| 237 | JPVT( PVT ) = JPVT( K ) |
---|
| 238 | JPVT( K ) = ITEMP |
---|
| 239 | VN1( PVT ) = VN1( K ) |
---|
| 240 | VN2( PVT ) = VN2( K ) |
---|
| 241 | END IF |
---|
| 242 | * |
---|
| 243 | * Apply previous Householder reflectors to column K: |
---|
| 244 | * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**T. |
---|
| 245 | * |
---|
| 246 | IF( K.GT.1 ) THEN |
---|
| 247 | CALL SGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ), |
---|
| 248 | $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 ) |
---|
| 249 | END IF |
---|
| 250 | * |
---|
| 251 | * Generate elementary reflector H(k). |
---|
| 252 | * |
---|
| 253 | IF( RK.LT.M ) THEN |
---|
| 254 | CALL SLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) |
---|
| 255 | ELSE |
---|
| 256 | CALL SLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) |
---|
| 257 | END IF |
---|
| 258 | * |
---|
| 259 | AKK = A( RK, K ) |
---|
| 260 | A( RK, K ) = ONE |
---|
| 261 | * |
---|
| 262 | * Compute Kth column of F: |
---|
| 263 | * |
---|
| 264 | * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**T*A(RK:M,K). |
---|
| 265 | * |
---|
| 266 | IF( K.LT.N ) THEN |
---|
| 267 | CALL SGEMV( 'Transpose', M-RK+1, N-K, TAU( K ), |
---|
| 268 | $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO, |
---|
| 269 | $ F( K+1, K ), 1 ) |
---|
| 270 | END IF |
---|
| 271 | * |
---|
| 272 | * Padding F(1:K,K) with zeros. |
---|
| 273 | * |
---|
| 274 | DO 20 J = 1, K |
---|
| 275 | F( J, K ) = ZERO |
---|
| 276 | 20 CONTINUE |
---|
| 277 | * |
---|
| 278 | * Incremental updating of F: |
---|
| 279 | * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**T |
---|
| 280 | * *A(RK:M,K). |
---|
| 281 | * |
---|
| 282 | IF( K.GT.1 ) THEN |
---|
| 283 | CALL SGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ), |
---|
| 284 | $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 ) |
---|
| 285 | * |
---|
| 286 | CALL SGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF, |
---|
| 287 | $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 ) |
---|
| 288 | END IF |
---|
| 289 | * |
---|
| 290 | * Update the current row of A: |
---|
| 291 | * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**T. |
---|
| 292 | * |
---|
| 293 | IF( K.LT.N ) THEN |
---|
| 294 | CALL SGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF, |
---|
| 295 | $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA ) |
---|
| 296 | END IF |
---|
| 297 | * |
---|
| 298 | * Update partial column norms. |
---|
| 299 | * |
---|
| 300 | IF( RK.LT.LASTRK ) THEN |
---|
| 301 | DO 30 J = K + 1, N |
---|
| 302 | IF( VN1( J ).NE.ZERO ) THEN |
---|
| 303 | * |
---|
| 304 | * NOTE: The following 4 lines follow from the analysis in |
---|
| 305 | * Lapack Working Note 176. |
---|
| 306 | * |
---|
| 307 | TEMP = ABS( A( RK, J ) ) / VN1( J ) |
---|
| 308 | TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) |
---|
| 309 | TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 |
---|
| 310 | IF( TEMP2 .LE. TOL3Z ) THEN |
---|
| 311 | VN2( J ) = REAL( LSTICC ) |
---|
| 312 | LSTICC = J |
---|
| 313 | ELSE |
---|
| 314 | VN1( J ) = VN1( J )*SQRT( TEMP ) |
---|
| 315 | END IF |
---|
| 316 | END IF |
---|
| 317 | 30 CONTINUE |
---|
| 318 | END IF |
---|
| 319 | * |
---|
| 320 | A( RK, K ) = AKK |
---|
| 321 | * |
---|
| 322 | * End of while loop. |
---|
| 323 | * |
---|
| 324 | GO TO 10 |
---|
| 325 | END IF |
---|
| 326 | KB = K |
---|
| 327 | RK = OFFSET + KB |
---|
| 328 | * |
---|
| 329 | * Apply the block reflector to the rest of the matrix: |
---|
| 330 | * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - |
---|
| 331 | * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**T. |
---|
| 332 | * |
---|
| 333 | IF( KB.LT.MIN( N, M-OFFSET ) ) THEN |
---|
| 334 | CALL SGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE, |
---|
| 335 | $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE, |
---|
| 336 | $ A( RK+1, KB+1 ), LDA ) |
---|
| 337 | END IF |
---|
| 338 | * |
---|
| 339 | * Recomputation of difficult columns. |
---|
| 340 | * |
---|
| 341 | 40 CONTINUE |
---|
| 342 | IF( LSTICC.GT.0 ) THEN |
---|
| 343 | ITEMP = NINT( VN2( LSTICC ) ) |
---|
| 344 | VN1( LSTICC ) = SNRM2( M-RK, A( RK+1, LSTICC ), 1 ) |
---|
| 345 | * |
---|
| 346 | * NOTE: The computation of VN1( LSTICC ) relies on the fact that |
---|
| 347 | * SNRM2 does not fail on vectors with norm below the value of |
---|
| 348 | * SQRT(DLAMCH('S')) |
---|
| 349 | * |
---|
| 350 | VN2( LSTICC ) = VN1( LSTICC ) |
---|
| 351 | LSTICC = ITEMP |
---|
| 352 | GO TO 40 |
---|
| 353 | END IF |
---|
| 354 | * |
---|
| 355 | RETURN |
---|
| 356 | * |
---|
| 357 | * End of SLAQPS |
---|
| 358 | * |
---|
| 359 | END |
---|