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