1 SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
2 $ WORK )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, M, N, OFFSET
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQP2 computes a QR factorization with column pivoting of
22 * the block A(OFFSET+1:M,1:N).
23 * The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
24 *
25 * Arguments
26 * =========
27 *
28 * M (input) INTEGER
29 * The number of rows of the matrix A. M >= 0.
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A. N >= 0.
33 *
34 * OFFSET (input) INTEGER
35 * The number of rows of the matrix A that must be pivoted
36 * but no factorized. OFFSET >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the M-by-N matrix A.
40 * On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
41 * the triangular factor obtained; the elements in block
42 * A(OFFSET+1:M,1:N) below the diagonal, together with the
43 * array TAU, represent the orthogonal matrix Q as a product of
44 * elementary reflectors. Block A(1:OFFSET,1:N) has been
45 * accordingly pivoted, but no factorized.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,M).
49 *
50 * JPVT (input/output) INTEGER array, dimension (N)
51 * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
52 * to the front of A*P (a leading column); if JPVT(i) = 0,
53 * the i-th column of A is a free column.
54 * On exit, if JPVT(i) = k, then the i-th column of A*P
55 * was the k-th column of A.
56 *
57 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
58 * The scalar factors of the elementary reflectors.
59 *
60 * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
61 * The vector with the partial column norms.
62 *
63 * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
64 * The vector with the exact column norms.
65 *
66 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
67 *
68 * Further Details
69 * ===============
70 *
71 * Based on contributions by
72 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
73 * X. Sun, Computer Science Dept., Duke University, USA
74 *
75 * Partial column norm updating strategy modified by
76 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
77 * University of Zagreb, Croatia.
78 * -- April 2011 --
79 * For more details see LAPACK Working Note 176.
80 * =====================================================================
81 *
82 * .. Parameters ..
83 DOUBLE PRECISION ZERO, ONE
84 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, ITEMP, J, MN, OFFPI, PVT
88 DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL DLARF, DLARFG, DSWAP
92 * ..
93 * .. Intrinsic Functions ..
94 INTRINSIC ABS, MAX, MIN, SQRT
95 * ..
96 * .. External Functions ..
97 INTEGER IDAMAX
98 DOUBLE PRECISION DLAMCH, DNRM2
99 EXTERNAL IDAMAX, DLAMCH, DNRM2
100 * ..
101 * .. Executable Statements ..
102 *
103 MN = MIN( M-OFFSET, N )
104 TOL3Z = SQRT(DLAMCH('Epsilon'))
105 *
106 * Compute factorization.
107 *
108 DO 20 I = 1, MN
109 *
110 OFFPI = OFFSET + I
111 *
112 * Determine ith pivot column and swap if necessary.
113 *
114 PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
115 *
116 IF( PVT.NE.I ) THEN
117 CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
118 ITEMP = JPVT( PVT )
119 JPVT( PVT ) = JPVT( I )
120 JPVT( I ) = ITEMP
121 VN1( PVT ) = VN1( I )
122 VN2( PVT ) = VN2( I )
123 END IF
124 *
125 * Generate elementary reflector H(i).
126 *
127 IF( OFFPI.LT.M ) THEN
128 CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
129 $ TAU( I ) )
130 ELSE
131 CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
132 END IF
133 *
134 IF( I.LE.N ) THEN
135 *
136 * Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
137 *
138 AII = A( OFFPI, I )
139 A( OFFPI, I ) = ONE
140 CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
141 $ TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
142 A( OFFPI, I ) = AII
143 END IF
144 *
145 * Update partial column norms.
146 *
147 DO 10 J = I + 1, N
148 IF( VN1( J ).NE.ZERO ) THEN
149 *
150 * NOTE: The following 4 lines follow from the analysis in
151 * Lapack Working Note 176.
152 *
153 TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
154 TEMP = MAX( TEMP, ZERO )
155 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
156 IF( TEMP2 .LE. TOL3Z ) THEN
157 IF( OFFPI.LT.M ) THEN
158 VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
159 VN2( J ) = VN1( J )
160 ELSE
161 VN1( J ) = ZERO
162 VN2( J ) = ZERO
163 END IF
164 ELSE
165 VN1( J ) = VN1( J )*SQRT( TEMP )
166 END IF
167 END IF
168 10 CONTINUE
169 *
170 20 CONTINUE
171 *
172 RETURN
173 *
174 * End of DLAQP2
175 *
176 END
2 $ WORK )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, M, N, OFFSET
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQP2 computes a QR factorization with column pivoting of
22 * the block A(OFFSET+1:M,1:N).
23 * The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
24 *
25 * Arguments
26 * =========
27 *
28 * M (input) INTEGER
29 * The number of rows of the matrix A. M >= 0.
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A. N >= 0.
33 *
34 * OFFSET (input) INTEGER
35 * The number of rows of the matrix A that must be pivoted
36 * but no factorized. OFFSET >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the M-by-N matrix A.
40 * On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
41 * the triangular factor obtained; the elements in block
42 * A(OFFSET+1:M,1:N) below the diagonal, together with the
43 * array TAU, represent the orthogonal matrix Q as a product of
44 * elementary reflectors. Block A(1:OFFSET,1:N) has been
45 * accordingly pivoted, but no factorized.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,M).
49 *
50 * JPVT (input/output) INTEGER array, dimension (N)
51 * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
52 * to the front of A*P (a leading column); if JPVT(i) = 0,
53 * the i-th column of A is a free column.
54 * On exit, if JPVT(i) = k, then the i-th column of A*P
55 * was the k-th column of A.
56 *
57 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
58 * The scalar factors of the elementary reflectors.
59 *
60 * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
61 * The vector with the partial column norms.
62 *
63 * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
64 * The vector with the exact column norms.
65 *
66 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
67 *
68 * Further Details
69 * ===============
70 *
71 * Based on contributions by
72 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
73 * X. Sun, Computer Science Dept., Duke University, USA
74 *
75 * Partial column norm updating strategy modified by
76 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
77 * University of Zagreb, Croatia.
78 * -- April 2011 --
79 * For more details see LAPACK Working Note 176.
80 * =====================================================================
81 *
82 * .. Parameters ..
83 DOUBLE PRECISION ZERO, ONE
84 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, ITEMP, J, MN, OFFPI, PVT
88 DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL DLARF, DLARFG, DSWAP
92 * ..
93 * .. Intrinsic Functions ..
94 INTRINSIC ABS, MAX, MIN, SQRT
95 * ..
96 * .. External Functions ..
97 INTEGER IDAMAX
98 DOUBLE PRECISION DLAMCH, DNRM2
99 EXTERNAL IDAMAX, DLAMCH, DNRM2
100 * ..
101 * .. Executable Statements ..
102 *
103 MN = MIN( M-OFFSET, N )
104 TOL3Z = SQRT(DLAMCH('Epsilon'))
105 *
106 * Compute factorization.
107 *
108 DO 20 I = 1, MN
109 *
110 OFFPI = OFFSET + I
111 *
112 * Determine ith pivot column and swap if necessary.
113 *
114 PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
115 *
116 IF( PVT.NE.I ) THEN
117 CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
118 ITEMP = JPVT( PVT )
119 JPVT( PVT ) = JPVT( I )
120 JPVT( I ) = ITEMP
121 VN1( PVT ) = VN1( I )
122 VN2( PVT ) = VN2( I )
123 END IF
124 *
125 * Generate elementary reflector H(i).
126 *
127 IF( OFFPI.LT.M ) THEN
128 CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
129 $ TAU( I ) )
130 ELSE
131 CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
132 END IF
133 *
134 IF( I.LE.N ) THEN
135 *
136 * Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
137 *
138 AII = A( OFFPI, I )
139 A( OFFPI, I ) = ONE
140 CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
141 $ TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
142 A( OFFPI, I ) = AII
143 END IF
144 *
145 * Update partial column norms.
146 *
147 DO 10 J = I + 1, N
148 IF( VN1( J ).NE.ZERO ) THEN
149 *
150 * NOTE: The following 4 lines follow from the analysis in
151 * Lapack Working Note 176.
152 *
153 TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
154 TEMP = MAX( TEMP, ZERO )
155 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
156 IF( TEMP2 .LE. TOL3Z ) THEN
157 IF( OFFPI.LT.M ) THEN
158 VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
159 VN2( J ) = VN1( J )
160 ELSE
161 VN1( J ) = ZERO
162 VN2( J ) = ZERO
163 END IF
164 ELSE
165 VN1( J ) = VN1( J )*SQRT( TEMP )
166 END IF
167 END IF
168 10 CONTINUE
169 *
170 20 CONTINUE
171 *
172 RETURN
173 *
174 * End of DLAQP2
175 *
176 END