LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaqp2.f
Go to the documentation of this file.
1*> \brief \b DLAQP2 computes a QR factorization with column pivoting of the matrix block.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAQP2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqp2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqp2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqp2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
22* WORK )
23*
24* .. Scalar Arguments ..
25* INTEGER LDA, M, N, OFFSET
26* ..
27* .. Array Arguments ..
28* INTEGER JPVT( * )
29* DOUBLE PRECISION A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
30* $ WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLAQP2 computes a QR factorization with column pivoting of
40*> the block A(OFFSET+1:M,1:N).
41*> The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] M
48*> \verbatim
49*> M is INTEGER
50*> The number of rows of the matrix A. M >= 0.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The number of columns of the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in] OFFSET
60*> \verbatim
61*> OFFSET is INTEGER
62*> The number of rows of the matrix A that must be pivoted
63*> but no factorized. OFFSET >= 0.
64*> \endverbatim
65*>
66*> \param[in,out] A
67*> \verbatim
68*> A is DOUBLE PRECISION array, dimension (LDA,N)
69*> On entry, the M-by-N matrix A.
70*> On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
71*> the triangular factor obtained; the elements in block
72*> A(OFFSET+1:M,1:N) below the diagonal, together with the
73*> array TAU, represent the orthogonal matrix Q as a product of
74*> elementary reflectors. Block A(1:OFFSET,1:N) has been
75*> accordingly pivoted, but no factorized.
76*> \endverbatim
77*>
78*> \param[in] LDA
79*> \verbatim
80*> LDA is INTEGER
81*> The leading dimension of the array A. LDA >= max(1,M).
82*> \endverbatim
83*>
84*> \param[in,out] JPVT
85*> \verbatim
86*> JPVT is INTEGER array, dimension (N)
87*> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
88*> to the front of A*P (a leading column); if JPVT(i) = 0,
89*> the i-th column of A is a free column.
90*> On exit, if JPVT(i) = k, then the i-th column of A*P
91*> was the k-th column of A.
92*> \endverbatim
93*>
94*> \param[out] TAU
95*> \verbatim
96*> TAU is DOUBLE PRECISION array, dimension (min(M,N))
97*> The scalar factors of the elementary reflectors.
98*> \endverbatim
99*>
100*> \param[in,out] VN1
101*> \verbatim
102*> VN1 is DOUBLE PRECISION array, dimension (N)
103*> The vector with the partial column norms.
104*> \endverbatim
105*>
106*> \param[in,out] VN2
107*> \verbatim
108*> VN2 is DOUBLE PRECISION array, dimension (N)
109*> The vector with the exact column norms.
110*> \endverbatim
111*>
112*> \param[out] WORK
113*> \verbatim
114*> WORK is DOUBLE PRECISION array, dimension (N)
115*> \endverbatim
116*
117* Authors:
118* ========
119*
120*> \author Univ. of Tennessee
121*> \author Univ. of California Berkeley
122*> \author Univ. of Colorado Denver
123*> \author NAG Ltd.
124*
125*> \ingroup laqp2
126*
127*> \par Contributors:
128* ==================
129*>
130*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
131*> X. Sun, Computer Science Dept., Duke University, USA
132*> \n
133*> Partial column norm updating strategy modified on April 2011
134*> Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
135*> University of Zagreb, Croatia.
136*
137*> \par References:
138* ================
139*>
140*> LAPACK Working Note 176
141*
142*> \htmlonly
143*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
144*> \endhtmlonly
145*
146* =====================================================================
147 SUBROUTINE dlaqp2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
148 $ WORK )
149*
150* -- LAPACK auxiliary routine --
151* -- LAPACK is a software package provided by Univ. of Tennessee, --
152* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153*
154* .. Scalar Arguments ..
155 INTEGER LDA, M, N, OFFSET
156* ..
157* .. Array Arguments ..
158 INTEGER JPVT( * )
159 DOUBLE PRECISION A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
160 $ work( * )
161* ..
162*
163* =====================================================================
164*
165* .. Parameters ..
166 DOUBLE PRECISION ZERO, ONE
167 parameter( zero = 0.0d+0, one = 1.0d+0 )
168* ..
169* .. Local Scalars ..
170 INTEGER I, ITEMP, J, MN, OFFPI, PVT
171 DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z
172* ..
173* .. External Subroutines ..
174 EXTERNAL dlarf, dlarfg, dswap
175* ..
176* .. Intrinsic Functions ..
177 INTRINSIC abs, max, min, sqrt
178* ..
179* .. External Functions ..
180 INTEGER IDAMAX
181 DOUBLE PRECISION DLAMCH, DNRM2
182 EXTERNAL idamax, dlamch, dnrm2
183* ..
184* .. Executable Statements ..
185*
186 mn = min( m-offset, n )
187 tol3z = sqrt(dlamch('Epsilon'))
188*
189* Compute factorization.
190*
191 DO 20 i = 1, mn
192*
193 offpi = offset + i
194*
195* Determine ith pivot column and swap if necessary.
196*
197 pvt = ( i-1 ) + idamax( n-i+1, vn1( i ), 1 )
198*
199 IF( pvt.NE.i ) THEN
200 CALL dswap( m, a( 1, pvt ), 1, a( 1, i ), 1 )
201 itemp = jpvt( pvt )
202 jpvt( pvt ) = jpvt( i )
203 jpvt( i ) = itemp
204 vn1( pvt ) = vn1( i )
205 vn2( pvt ) = vn2( i )
206 END IF
207*
208* Generate elementary reflector H(i).
209*
210 IF( offpi.LT.m ) THEN
211 CALL dlarfg( m-offpi+1, a( offpi, i ), a( offpi+1, i ), 1,
212 $ tau( i ) )
213 ELSE
214 CALL dlarfg( 1, a( m, i ), a( m, i ), 1, tau( i ) )
215 END IF
216*
217 IF( i.LT.n ) THEN
218*
219* Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
220*
221 aii = a( offpi, i )
222 a( offpi, i ) = one
223 CALL dlarf( 'Left', m-offpi+1, n-i, a( offpi, i ), 1,
224 $ tau( i ), a( offpi, i+1 ), lda, work( 1 ) )
225 a( offpi, i ) = aii
226 END IF
227*
228* Update partial column norms.
229*
230 DO 10 j = i + 1, n
231 IF( vn1( j ).NE.zero ) THEN
232*
233* NOTE: The following 4 lines follow from the analysis in
234* Lapack Working Note 176.
235*
236 temp = one - ( abs( a( offpi, j ) ) / vn1( j ) )**2
237 temp = max( temp, zero )
238 temp2 = temp*( vn1( j ) / vn2( j ) )**2
239 IF( temp2 .LE. tol3z ) THEN
240 IF( offpi.LT.m ) THEN
241 vn1( j ) = dnrm2( m-offpi, a( offpi+1, j ), 1 )
242 vn2( j ) = vn1( j )
243 ELSE
244 vn1( j ) = zero
245 vn2( j ) = zero
246 END IF
247 ELSE
248 vn1( j ) = vn1( j )*sqrt( temp )
249 END IF
250 END IF
251 10 CONTINUE
252*
253 20 CONTINUE
254*
255 RETURN
256*
257* End of DLAQP2
258*
259 END
subroutine dlaqp2(m, n, offset, a, lda, jpvt, tau, vn1, vn2, work)
DLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition dlaqp2.f:149
subroutine dlarf(side, m, n, v, incv, tau, c, ldc, work)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition dlarf.f:124
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82