root / src / lapack / double / dlaqps.f @ 1
Historique | Voir | Annoter | Télécharger (8,08 ko)
| 1 | 1 | equemene | SUBROUTINE DLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, |
|---|---|---|---|
| 2 | 1 | equemene | $ VN2, AUXV, F, LDF ) |
| 3 | 1 | equemene | * |
| 4 | 1 | equemene | * -- LAPACK auxiliary routine (version 3.2.2) -- |
| 5 | 1 | equemene | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 6 | 1 | equemene | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 7 | 1 | equemene | * June 2010 |
| 8 | 1 | equemene | * |
| 9 | 1 | equemene | * .. Scalar Arguments .. |
| 10 | 1 | equemene | INTEGER KB, LDA, LDF, M, N, NB, OFFSET |
| 11 | 1 | equemene | * .. |
| 12 | 1 | equemene | * .. Array Arguments .. |
| 13 | 1 | equemene | INTEGER JPVT( * ) |
| 14 | 1 | equemene | DOUBLE PRECISION A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), |
| 15 | 1 | equemene | $ VN1( * ), VN2( * ) |
| 16 | 1 | equemene | * .. |
| 17 | 1 | equemene | * |
| 18 | 1 | equemene | * Purpose |
| 19 | 1 | equemene | * ======= |
| 20 | 1 | equemene | * |
| 21 | 1 | equemene | * DLAQPS computes a step of QR factorization with column pivoting |
| 22 | 1 | equemene | * of a real M-by-N matrix A by using Blas-3. It tries to factorize |
| 23 | 1 | equemene | * NB columns from A starting from the row OFFSET+1, and updates all |
| 24 | 1 | equemene | * of the matrix with Blas-3 xGEMM. |
| 25 | 1 | equemene | * |
| 26 | 1 | equemene | * In some cases, due to catastrophic cancellations, it cannot |
| 27 | 1 | equemene | * factorize NB columns. Hence, the actual number of factorized |
| 28 | 1 | equemene | * columns is returned in KB. |
| 29 | 1 | equemene | * |
| 30 | 1 | equemene | * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. |
| 31 | 1 | equemene | * |
| 32 | 1 | equemene | * Arguments |
| 33 | 1 | equemene | * ========= |
| 34 | 1 | equemene | * |
| 35 | 1 | equemene | * M (input) INTEGER |
| 36 | 1 | equemene | * The number of rows of the matrix A. M >= 0. |
| 37 | 1 | equemene | * |
| 38 | 1 | equemene | * N (input) INTEGER |
| 39 | 1 | equemene | * The number of columns of the matrix A. N >= 0 |
| 40 | 1 | equemene | * |
| 41 | 1 | equemene | * OFFSET (input) INTEGER |
| 42 | 1 | equemene | * The number of rows of A that have been factorized in |
| 43 | 1 | equemene | * previous steps. |
| 44 | 1 | equemene | * |
| 45 | 1 | equemene | * NB (input) INTEGER |
| 46 | 1 | equemene | * The number of columns to factorize. |
| 47 | 1 | equemene | * |
| 48 | 1 | equemene | * KB (output) INTEGER |
| 49 | 1 | equemene | * The number of columns actually factorized. |
| 50 | 1 | equemene | * |
| 51 | 1 | equemene | * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) |
| 52 | 1 | equemene | * On entry, the M-by-N matrix A. |
| 53 | 1 | equemene | * On exit, block A(OFFSET+1:M,1:KB) is the triangular |
| 54 | 1 | equemene | * factor obtained and block A(1:OFFSET,1:N) has been |
| 55 | 1 | equemene | * accordingly pivoted, but no factorized. |
| 56 | 1 | equemene | * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has |
| 57 | 1 | equemene | * been updated. |
| 58 | 1 | equemene | * |
| 59 | 1 | equemene | * LDA (input) INTEGER |
| 60 | 1 | equemene | * The leading dimension of the array A. LDA >= max(1,M). |
| 61 | 1 | equemene | * |
| 62 | 1 | equemene | * JPVT (input/output) INTEGER array, dimension (N) |
| 63 | 1 | equemene | * JPVT(I) = K <==> Column K of the full matrix A has been |
| 64 | 1 | equemene | * permuted into position I in AP. |
| 65 | 1 | equemene | * |
| 66 | 1 | equemene | * TAU (output) DOUBLE PRECISION array, dimension (KB) |
| 67 | 1 | equemene | * The scalar factors of the elementary reflectors. |
| 68 | 1 | equemene | * |
| 69 | 1 | equemene | * VN1 (input/output) DOUBLE PRECISION array, dimension (N) |
| 70 | 1 | equemene | * The vector with the partial column norms. |
| 71 | 1 | equemene | * |
| 72 | 1 | equemene | * VN2 (input/output) DOUBLE PRECISION array, dimension (N) |
| 73 | 1 | equemene | * The vector with the exact column norms. |
| 74 | 1 | equemene | * |
| 75 | 1 | equemene | * AUXV (input/output) DOUBLE PRECISION array, dimension (NB) |
| 76 | 1 | equemene | * Auxiliar vector. |
| 77 | 1 | equemene | * |
| 78 | 1 | equemene | * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB) |
| 79 | 1 | equemene | * Matrix F' = L*Y'*A. |
| 80 | 1 | equemene | * |
| 81 | 1 | equemene | * LDF (input) INTEGER |
| 82 | 1 | equemene | * The leading dimension of the array F. LDF >= max(1,N). |
| 83 | 1 | equemene | * |
| 84 | 1 | equemene | * Further Details |
| 85 | 1 | equemene | * =============== |
| 86 | 1 | equemene | * |
| 87 | 1 | equemene | * Based on contributions by |
| 88 | 1 | equemene | * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain |
| 89 | 1 | equemene | * X. Sun, Computer Science Dept., Duke University, USA |
| 90 | 1 | equemene | * |
| 91 | 1 | equemene | * Partial column norm updating strategy modified by |
| 92 | 1 | equemene | * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, |
| 93 | 1 | equemene | * University of Zagreb, Croatia. |
| 94 | 1 | equemene | * June 2010 |
| 95 | 1 | equemene | * For more details see LAPACK Working Note 176. |
| 96 | 1 | equemene | * ===================================================================== |
| 97 | 1 | equemene | * |
| 98 | 1 | equemene | * .. Parameters .. |
| 99 | 1 | equemene | DOUBLE PRECISION ZERO, ONE |
| 100 | 1 | equemene | PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
| 101 | 1 | equemene | * .. |
| 102 | 1 | equemene | * .. Local Scalars .. |
| 103 | 1 | equemene | INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK |
| 104 | 1 | equemene | DOUBLE PRECISION AKK, TEMP, TEMP2, TOL3Z |
| 105 | 1 | equemene | * .. |
| 106 | 1 | equemene | * .. External Subroutines .. |
| 107 | 1 | equemene | EXTERNAL DGEMM, DGEMV, DLARFG, DSWAP |
| 108 | 1 | equemene | * .. |
| 109 | 1 | equemene | * .. Intrinsic Functions .. |
| 110 | 1 | equemene | INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT |
| 111 | 1 | equemene | * .. |
| 112 | 1 | equemene | * .. External Functions .. |
| 113 | 1 | equemene | INTEGER IDAMAX |
| 114 | 1 | equemene | DOUBLE PRECISION DLAMCH, DNRM2 |
| 115 | 1 | equemene | EXTERNAL IDAMAX, DLAMCH, DNRM2 |
| 116 | 1 | equemene | * .. |
| 117 | 1 | equemene | * .. Executable Statements .. |
| 118 | 1 | equemene | * |
| 119 | 1 | equemene | LASTRK = MIN( M, N+OFFSET ) |
| 120 | 1 | equemene | LSTICC = 0 |
| 121 | 1 | equemene | K = 0 |
| 122 | 1 | equemene | TOL3Z = SQRT(DLAMCH('Epsilon'))
|
| 123 | 1 | equemene | * |
| 124 | 1 | equemene | * Beginning of while loop. |
| 125 | 1 | equemene | * |
| 126 | 1 | equemene | 10 CONTINUE |
| 127 | 1 | equemene | IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN |
| 128 | 1 | equemene | K = K + 1 |
| 129 | 1 | equemene | RK = OFFSET + K |
| 130 | 1 | equemene | * |
| 131 | 1 | equemene | * Determine ith pivot column and swap if necessary |
| 132 | 1 | equemene | * |
| 133 | 1 | equemene | PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 ) |
| 134 | 1 | equemene | IF( PVT.NE.K ) THEN |
| 135 | 1 | equemene | CALL DSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) |
| 136 | 1 | equemene | CALL DSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) |
| 137 | 1 | equemene | ITEMP = JPVT( PVT ) |
| 138 | 1 | equemene | JPVT( PVT ) = JPVT( K ) |
| 139 | 1 | equemene | JPVT( K ) = ITEMP |
| 140 | 1 | equemene | VN1( PVT ) = VN1( K ) |
| 141 | 1 | equemene | VN2( PVT ) = VN2( K ) |
| 142 | 1 | equemene | END IF |
| 143 | 1 | equemene | * |
| 144 | 1 | equemene | * Apply previous Householder reflectors to column K: |
| 145 | 1 | equemene | * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. |
| 146 | 1 | equemene | * |
| 147 | 1 | equemene | IF( K.GT.1 ) THEN |
| 148 | 1 | equemene | CALL DGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ), |
| 149 | 1 | equemene | $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 ) |
| 150 | 1 | equemene | END IF |
| 151 | 1 | equemene | * |
| 152 | 1 | equemene | * Generate elementary reflector H(k). |
| 153 | 1 | equemene | * |
| 154 | 1 | equemene | IF( RK.LT.M ) THEN |
| 155 | 1 | equemene | CALL DLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) |
| 156 | 1 | equemene | ELSE |
| 157 | 1 | equemene | CALL DLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) |
| 158 | 1 | equemene | END IF |
| 159 | 1 | equemene | * |
| 160 | 1 | equemene | AKK = A( RK, K ) |
| 161 | 1 | equemene | A( RK, K ) = ONE |
| 162 | 1 | equemene | * |
| 163 | 1 | equemene | * Compute Kth column of F: |
| 164 | 1 | equemene | * |
| 165 | 1 | equemene | * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). |
| 166 | 1 | equemene | * |
| 167 | 1 | equemene | IF( K.LT.N ) THEN |
| 168 | 1 | equemene | CALL DGEMV( 'Transpose', M-RK+1, N-K, TAU( K ), |
| 169 | 1 | equemene | $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO, |
| 170 | 1 | equemene | $ F( K+1, K ), 1 ) |
| 171 | 1 | equemene | END IF |
| 172 | 1 | equemene | * |
| 173 | 1 | equemene | * Padding F(1:K,K) with zeros. |
| 174 | 1 | equemene | * |
| 175 | 1 | equemene | DO 20 J = 1, K |
| 176 | 1 | equemene | F( J, K ) = ZERO |
| 177 | 1 | equemene | 20 CONTINUE |
| 178 | 1 | equemene | * |
| 179 | 1 | equemene | * Incremental updating of F: |
| 180 | 1 | equemene | * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' |
| 181 | 1 | equemene | * *A(RK:M,K). |
| 182 | 1 | equemene | * |
| 183 | 1 | equemene | IF( K.GT.1 ) THEN |
| 184 | 1 | equemene | CALL DGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ), |
| 185 | 1 | equemene | $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 ) |
| 186 | 1 | equemene | * |
| 187 | 1 | equemene | CALL DGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF, |
| 188 | 1 | equemene | $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 ) |
| 189 | 1 | equemene | END IF |
| 190 | 1 | equemene | * |
| 191 | 1 | equemene | * Update the current row of A: |
| 192 | 1 | equemene | * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. |
| 193 | 1 | equemene | * |
| 194 | 1 | equemene | IF( K.LT.N ) THEN |
| 195 | 1 | equemene | CALL DGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF, |
| 196 | 1 | equemene | $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA ) |
| 197 | 1 | equemene | END IF |
| 198 | 1 | equemene | * |
| 199 | 1 | equemene | * Update partial column norms. |
| 200 | 1 | equemene | * |
| 201 | 1 | equemene | IF( RK.LT.LASTRK ) THEN |
| 202 | 1 | equemene | DO 30 J = K + 1, N |
| 203 | 1 | equemene | IF( VN1( J ).NE.ZERO ) THEN |
| 204 | 1 | equemene | * |
| 205 | 1 | equemene | * NOTE: The following 4 lines follow from the analysis in |
| 206 | 1 | equemene | * Lapack Working Note 176. |
| 207 | 1 | equemene | * |
| 208 | 1 | equemene | TEMP = ABS( A( RK, J ) ) / VN1( J ) |
| 209 | 1 | equemene | TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) |
| 210 | 1 | equemene | TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 |
| 211 | 1 | equemene | IF( TEMP2 .LE. TOL3Z ) THEN |
| 212 | 1 | equemene | VN2( J ) = DBLE( LSTICC ) |
| 213 | 1 | equemene | LSTICC = J |
| 214 | 1 | equemene | ELSE |
| 215 | 1 | equemene | VN1( J ) = VN1( J )*SQRT( TEMP ) |
| 216 | 1 | equemene | END IF |
| 217 | 1 | equemene | END IF |
| 218 | 1 | equemene | 30 CONTINUE |
| 219 | 1 | equemene | END IF |
| 220 | 1 | equemene | * |
| 221 | 1 | equemene | A( RK, K ) = AKK |
| 222 | 1 | equemene | * |
| 223 | 1 | equemene | * End of while loop. |
| 224 | 1 | equemene | * |
| 225 | 1 | equemene | GO TO 10 |
| 226 | 1 | equemene | END IF |
| 227 | 1 | equemene | KB = K |
| 228 | 1 | equemene | RK = OFFSET + KB |
| 229 | 1 | equemene | * |
| 230 | 1 | equemene | * Apply the block reflector to the rest of the matrix: |
| 231 | 1 | equemene | * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - |
| 232 | 1 | equemene | * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. |
| 233 | 1 | equemene | * |
| 234 | 1 | equemene | IF( KB.LT.MIN( N, M-OFFSET ) ) THEN |
| 235 | 1 | equemene | CALL DGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE, |
| 236 | 1 | equemene | $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE, |
| 237 | 1 | equemene | $ A( RK+1, KB+1 ), LDA ) |
| 238 | 1 | equemene | END IF |
| 239 | 1 | equemene | * |
| 240 | 1 | equemene | * Recomputation of difficult columns. |
| 241 | 1 | equemene | * |
| 242 | 1 | equemene | 40 CONTINUE |
| 243 | 1 | equemene | IF( LSTICC.GT.0 ) THEN |
| 244 | 1 | equemene | ITEMP = NINT( VN2( LSTICC ) ) |
| 245 | 1 | equemene | VN1( LSTICC ) = DNRM2( M-RK, A( RK+1, LSTICC ), 1 ) |
| 246 | 1 | equemene | * |
| 247 | 1 | equemene | * NOTE: The computation of VN1( LSTICC ) relies on the fact that |
| 248 | 1 | equemene | * SNRM2 does not fail on vectors with norm below the value of |
| 249 | 1 | equemene | * SQRT(DLAMCH('S'))
|
| 250 | 1 | equemene | * |
| 251 | 1 | equemene | VN2( LSTICC ) = VN1( LSTICC ) |
| 252 | 1 | equemene | LSTICC = ITEMP |
| 253 | 1 | equemene | GO TO 40 |
| 254 | 1 | equemene | END IF |
| 255 | 1 | equemene | * |
| 256 | 1 | equemene | RETURN |
| 257 | 1 | equemene | * |
| 258 | 1 | equemene | * End of DLAQPS |
| 259 | 1 | equemene | * |
| 260 | 1 | equemene | END |