|
2329
|
1 SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 INTEGER K, LDA, LDT, LDY, N, NB |
|
|
9 * .. |
|
|
10 * .. Array Arguments .. |
|
|
11 DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), |
|
|
12 $ Y( LDY, NB ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * DLAHRD reduces the first NB columns of a real general n-by-(n-k+1) |
|
|
19 * matrix A so that elements below the k-th subdiagonal are zero. The |
|
|
20 * reduction is performed by an orthogonal similarity transformation |
|
|
21 * Q' * A * Q. The routine returns the matrices V and T which determine |
|
|
22 * Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. |
|
|
23 * |
|
7034
|
24 * This is an OBSOLETE auxiliary routine. |
|
|
25 * This routine will be 'deprecated' in a future release. |
|
|
26 * Please use the new routine DLAHR2 instead. |
|
2329
|
27 * |
|
|
28 * Arguments |
|
|
29 * ========= |
|
|
30 * |
|
|
31 * N (input) INTEGER |
|
|
32 * The order of the matrix A. |
|
|
33 * |
|
|
34 * K (input) INTEGER |
|
|
35 * The offset for the reduction. Elements below the k-th |
|
|
36 * subdiagonal in the first NB columns are reduced to zero. |
|
|
37 * |
|
|
38 * NB (input) INTEGER |
|
|
39 * The number of columns to be reduced. |
|
|
40 * |
|
|
41 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1) |
|
|
42 * On entry, the n-by-(n-k+1) general matrix A. |
|
|
43 * On exit, the elements on and above the k-th subdiagonal in |
|
|
44 * the first NB columns are overwritten with the corresponding |
|
|
45 * elements of the reduced matrix; the elements below the k-th |
|
|
46 * subdiagonal, with the array TAU, represent the matrix Q as a |
|
|
47 * product of elementary reflectors. The other columns of A are |
|
|
48 * unchanged. See Further Details. |
|
|
49 * |
|
|
50 * LDA (input) INTEGER |
|
|
51 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
52 * |
|
|
53 * TAU (output) DOUBLE PRECISION array, dimension (NB) |
|
|
54 * The scalar factors of the elementary reflectors. See Further |
|
|
55 * Details. |
|
|
56 * |
|
3333
|
57 * T (output) DOUBLE PRECISION array, dimension (LDT,NB) |
|
2329
|
58 * The upper triangular matrix T. |
|
|
59 * |
|
|
60 * LDT (input) INTEGER |
|
|
61 * The leading dimension of the array T. LDT >= NB. |
|
|
62 * |
|
|
63 * Y (output) DOUBLE PRECISION array, dimension (LDY,NB) |
|
|
64 * The n-by-nb matrix Y. |
|
|
65 * |
|
|
66 * LDY (input) INTEGER |
|
|
67 * The leading dimension of the array Y. LDY >= N. |
|
|
68 * |
|
|
69 * Further Details |
|
|
70 * =============== |
|
|
71 * |
|
|
72 * The matrix Q is represented as a product of nb elementary reflectors |
|
|
73 * |
|
|
74 * Q = H(1) H(2) . . . H(nb). |
|
|
75 * |
|
|
76 * Each H(i) has the form |
|
|
77 * |
|
|
78 * H(i) = I - tau * v * v' |
|
|
79 * |
|
|
80 * where tau is a real scalar, and v is a real vector with |
|
|
81 * v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in |
|
|
82 * A(i+k+1:n,i), and tau in TAU(i). |
|
|
83 * |
|
|
84 * The elements of the vectors v together form the (n-k+1)-by-nb matrix |
|
|
85 * V which is needed, with T and Y, to apply the transformation to the |
|
|
86 * unreduced part of the matrix, using an update of the form: |
|
|
87 * A := (I - V*T*V') * (A - Y*V'). |
|
|
88 * |
|
|
89 * The contents of A on exit are illustrated by the following example |
|
|
90 * with n = 7, k = 3 and nb = 2: |
|
|
91 * |
|
|
92 * ( a h a a a ) |
|
|
93 * ( a h a a a ) |
|
|
94 * ( a h a a a ) |
|
|
95 * ( h h a a a ) |
|
|
96 * ( v1 h a a a ) |
|
|
97 * ( v1 v2 a a a ) |
|
|
98 * ( v1 v2 a a a ) |
|
|
99 * |
|
|
100 * where a denotes an element of the original matrix A, h denotes a |
|
|
101 * modified element of the upper Hessenberg matrix H, and vi denotes an |
|
|
102 * element of the vector defining H(i). |
|
|
103 * |
|
|
104 * ===================================================================== |
|
|
105 * |
|
|
106 * .. Parameters .. |
|
|
107 DOUBLE PRECISION ZERO, ONE |
|
|
108 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
|
|
109 * .. |
|
|
110 * .. Local Scalars .. |
|
|
111 INTEGER I |
|
|
112 DOUBLE PRECISION EI |
|
|
113 * .. |
|
|
114 * .. External Subroutines .. |
|
|
115 EXTERNAL DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV |
|
|
116 * .. |
|
|
117 * .. Intrinsic Functions .. |
|
|
118 INTRINSIC MIN |
|
|
119 * .. |
|
|
120 * .. Executable Statements .. |
|
|
121 * |
|
|
122 * Quick return if possible |
|
|
123 * |
|
|
124 IF( N.LE.1 ) |
|
|
125 $ RETURN |
|
|
126 * |
|
|
127 DO 10 I = 1, NB |
|
|
128 IF( I.GT.1 ) THEN |
|
|
129 * |
|
|
130 * Update A(1:n,i) |
|
|
131 * |
|
|
132 * Compute i-th column of A - Y * V' |
|
|
133 * |
|
|
134 CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, |
|
|
135 $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 ) |
|
|
136 * |
|
|
137 * Apply I - V * T' * V' to this column (call it b) from the |
|
|
138 * left, using the last column of T as workspace |
|
|
139 * |
|
|
140 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) |
|
|
141 * ( V2 ) ( b2 ) |
|
|
142 * |
|
|
143 * where V1 is unit lower triangular |
|
|
144 * |
|
|
145 * w := V1' * b1 |
|
|
146 * |
|
|
147 CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) |
|
|
148 CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ), |
|
|
149 $ LDA, T( 1, NB ), 1 ) |
|
|
150 * |
|
|
151 * w := w + V2'*b2 |
|
|
152 * |
|
|
153 CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), |
|
|
154 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) |
|
|
155 * |
|
|
156 * w := T'*w |
|
|
157 * |
|
|
158 CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT, |
|
|
159 $ T( 1, NB ), 1 ) |
|
|
160 * |
|
|
161 * b2 := b2 - V2*w |
|
|
162 * |
|
|
163 CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ), |
|
|
164 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) |
|
|
165 * |
|
|
166 * b1 := b1 - V1*w |
|
|
167 * |
|
|
168 CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1, |
|
|
169 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) |
|
|
170 CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) |
|
|
171 * |
|
|
172 A( K+I-1, I-1 ) = EI |
|
|
173 END IF |
|
|
174 * |
|
|
175 * Generate the elementary reflector H(i) to annihilate |
|
|
176 * A(k+i+1:n,i) |
|
|
177 * |
|
|
178 CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, |
|
|
179 $ TAU( I ) ) |
|
|
180 EI = A( K+I, I ) |
|
|
181 A( K+I, I ) = ONE |
|
|
182 * |
|
|
183 * Compute Y(1:n,i) |
|
|
184 * |
|
|
185 CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA, |
|
|
186 $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 ) |
|
|
187 CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA, |
|
|
188 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) |
|
|
189 CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1, |
|
|
190 $ ONE, Y( 1, I ), 1 ) |
|
|
191 CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 ) |
|
|
192 * |
|
|
193 * Compute T(1:i,i) |
|
|
194 * |
|
|
195 CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) |
|
|
196 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT, |
|
|
197 $ T( 1, I ), 1 ) |
|
|
198 T( I, I ) = TAU( I ) |
|
|
199 * |
|
|
200 10 CONTINUE |
|
|
201 A( K+NB, NB ) = EI |
|
|
202 * |
|
|
203 RETURN |
|
|
204 * |
|
|
205 * End of DLAHRD |
|
|
206 * |
|
|
207 END |
