001: SUBROUTINE DLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1,
002: $ VN2, AUXV, F, LDF )
003: *
004: * -- LAPACK auxiliary routine (version 3.2) --
005: * -- LAPACK is a software package provided by Univ. of Tennessee, --
006: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: * November 2006
008: *
009: * .. Scalar Arguments ..
010: INTEGER KB, LDA, LDF, M, N, NB, OFFSET
011: * ..
012: * .. Array Arguments ..
013: INTEGER JPVT( * )
014: DOUBLE PRECISION A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ),
015: $ VN1( * ), VN2( * )
016: * ..
017: *
018: * Purpose
019: * =======
020: *
021: * DLAQPS computes a step of QR factorization with column pivoting
022: * of a real M-by-N matrix A by using Blas-3. It tries to factorize
023: * NB columns from A starting from the row OFFSET+1, and updates all
024: * of the matrix with Blas-3 xGEMM.
025: *
026: * In some cases, due to catastrophic cancellations, it cannot
027: * factorize NB columns. Hence, the actual number of factorized
028: * columns is returned in KB.
029: *
030: * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
031: *
032: * Arguments
033: * =========
034: *
035: * M (input) INTEGER
036: * The number of rows of the matrix A. M >= 0.
037: *
038: * N (input) INTEGER
039: * The number of columns of the matrix A. N >= 0
040: *
041: * OFFSET (input) INTEGER
042: * The number of rows of A that have been factorized in
043: * previous steps.
044: *
045: * NB (input) INTEGER
046: * The number of columns to factorize.
047: *
048: * KB (output) INTEGER
049: * The number of columns actually factorized.
050: *
051: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
052: * On entry, the M-by-N matrix A.
053: * On exit, block A(OFFSET+1:M,1:KB) is the triangular
054: * factor obtained and block A(1:OFFSET,1:N) has been
055: * accordingly pivoted, but no factorized.
056: * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
057: * been updated.
058: *
059: * LDA (input) INTEGER
060: * The leading dimension of the array A. LDA >= max(1,M).
061: *
062: * JPVT (input/output) INTEGER array, dimension (N)
063: * JPVT(I) = K <==> Column K of the full matrix A has been
064: * permuted into position I in AP.
065: *
066: * TAU (output) DOUBLE PRECISION array, dimension (KB)
067: * The scalar factors of the elementary reflectors.
068: *
069: * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
070: * The vector with the partial column norms.
071: *
072: * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
073: * The vector with the exact column norms.
074: *
075: * AUXV (input/output) DOUBLE PRECISION array, dimension (NB)
076: * Auxiliar vector.
077: *
078: * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB)
079: * Matrix F' = L*Y'*A.
080: *
081: * LDF (input) INTEGER
082: * The leading dimension of the array F. LDF >= max(1,N).
083: *
084: * Further Details
085: * ===============
086: *
087: * Based on contributions by
088: * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
089: * X. Sun, Computer Science Dept., Duke University, USA
090: *
091: * Partial column norm updating strategy modified by
092: * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
093: * University of Zagreb, Croatia.
094: * June 2006.
095: * For more details see LAPACK Working Note 176.
096: * =====================================================================
097: *
098: * .. Parameters ..
099: 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 ..
107: EXTERNAL DGEMM, DGEMV, DLARFP, DSWAP
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
155: CALL DLARFP( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) )
156: ELSE
157: CALL DLARFP( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) )
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
261: