1 SUBROUTINE ZLAQP2( 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 VN1( * ), VN2( * )
15 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZLAQP2 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 COMPLEX*16 CONE
85 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
86 $ CONE = ( 1.0D+0, 0.0D+0 ) )
87 * ..
88 * .. Local Scalars ..
89 INTEGER I, ITEMP, J, MN, OFFPI, PVT
90 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
91 COMPLEX*16 AII
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZLARF, ZLARFG, ZSWAP
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC ABS, DCONJG, MAX, MIN, SQRT
98 * ..
99 * .. External Functions ..
100 INTEGER IDAMAX
101 DOUBLE PRECISION DLAMCH, DZNRM2
102 EXTERNAL IDAMAX, DLAMCH, DZNRM2
103 * ..
104 * .. Executable Statements ..
105 *
106 MN = MIN( M-OFFSET, N )
107 TOL3Z = SQRT(DLAMCH('Epsilon'))
108 *
109 * Compute factorization.
110 *
111 DO 20 I = 1, MN
112 *
113 OFFPI = OFFSET + I
114 *
115 * Determine ith pivot column and swap if necessary.
116 *
117 PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
118 *
119 IF( PVT.NE.I ) THEN
120 CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
121 ITEMP = JPVT( PVT )
122 JPVT( PVT ) = JPVT( I )
123 JPVT( I ) = ITEMP
124 VN1( PVT ) = VN1( I )
125 VN2( PVT ) = VN2( I )
126 END IF
127 *
128 * Generate elementary reflector H(i).
129 *
130 IF( OFFPI.LT.M ) THEN
131 CALL ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
132 $ TAU( I ) )
133 ELSE
134 CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
135 END IF
136 *
137 IF( I.LT.N ) THEN
138 *
139 * Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
140 *
141 AII = A( OFFPI, I )
142 A( OFFPI, I ) = CONE
143 CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
144 $ DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
145 $ WORK( 1 ) )
146 A( OFFPI, I ) = AII
147 END IF
148 *
149 * Update partial column norms.
150 *
151 DO 10 J = I + 1, N
152 IF( VN1( J ).NE.ZERO ) THEN
153 *
154 * NOTE: The following 4 lines follow from the analysis in
155 * Lapack Working Note 176.
156 *
157 TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
158 TEMP = MAX( TEMP, ZERO )
159 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
160 IF( TEMP2 .LE. TOL3Z ) THEN
161 IF( OFFPI.LT.M ) THEN
162 VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
163 VN2( J ) = VN1( J )
164 ELSE
165 VN1( J ) = ZERO
166 VN2( J ) = ZERO
167 END IF
168 ELSE
169 VN1( J ) = VN1( J )*SQRT( TEMP )
170 END IF
171 END IF
172 10 CONTINUE
173 *
174 20 CONTINUE
175 *
176 RETURN
177 *
178 * End of ZLAQP2
179 *
180 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 VN1( * ), VN2( * )
15 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZLAQP2 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 COMPLEX*16 CONE
85 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
86 $ CONE = ( 1.0D+0, 0.0D+0 ) )
87 * ..
88 * .. Local Scalars ..
89 INTEGER I, ITEMP, J, MN, OFFPI, PVT
90 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
91 COMPLEX*16 AII
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZLARF, ZLARFG, ZSWAP
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC ABS, DCONJG, MAX, MIN, SQRT
98 * ..
99 * .. External Functions ..
100 INTEGER IDAMAX
101 DOUBLE PRECISION DLAMCH, DZNRM2
102 EXTERNAL IDAMAX, DLAMCH, DZNRM2
103 * ..
104 * .. Executable Statements ..
105 *
106 MN = MIN( M-OFFSET, N )
107 TOL3Z = SQRT(DLAMCH('Epsilon'))
108 *
109 * Compute factorization.
110 *
111 DO 20 I = 1, MN
112 *
113 OFFPI = OFFSET + I
114 *
115 * Determine ith pivot column and swap if necessary.
116 *
117 PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
118 *
119 IF( PVT.NE.I ) THEN
120 CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
121 ITEMP = JPVT( PVT )
122 JPVT( PVT ) = JPVT( I )
123 JPVT( I ) = ITEMP
124 VN1( PVT ) = VN1( I )
125 VN2( PVT ) = VN2( I )
126 END IF
127 *
128 * Generate elementary reflector H(i).
129 *
130 IF( OFFPI.LT.M ) THEN
131 CALL ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
132 $ TAU( I ) )
133 ELSE
134 CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
135 END IF
136 *
137 IF( I.LT.N ) THEN
138 *
139 * Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
140 *
141 AII = A( OFFPI, I )
142 A( OFFPI, I ) = CONE
143 CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
144 $ DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
145 $ WORK( 1 ) )
146 A( OFFPI, I ) = AII
147 END IF
148 *
149 * Update partial column norms.
150 *
151 DO 10 J = I + 1, N
152 IF( VN1( J ).NE.ZERO ) THEN
153 *
154 * NOTE: The following 4 lines follow from the analysis in
155 * Lapack Working Note 176.
156 *
157 TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
158 TEMP = MAX( TEMP, ZERO )
159 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
160 IF( TEMP2 .LE. TOL3Z ) THEN
161 IF( OFFPI.LT.M ) THEN
162 VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
163 VN2( J ) = VN1( J )
164 ELSE
165 VN1( J ) = ZERO
166 VN2( J ) = ZERO
167 END IF
168 ELSE
169 VN1( J ) = VN1( J )*SQRT( TEMP )
170 END IF
171 END IF
172 10 CONTINUE
173 *
174 20 CONTINUE
175 *
176 RETURN
177 *
178 * End of ZLAQP2
179 *
180 END