1: *> \brief \b DLAHRD reduces the first nb columns of a general rectangular matrix A so that elements below the k-th subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAHRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahrd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER K, LDA, LDT, LDY, N, NB
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ),
28: * $ Y( LDY, NB )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> This routine is deprecated and has been replaced by routine DLAHR2.
38: *>
39: *> DLAHRD reduces the first NB columns of a real general n-by-(n-k+1)
40: *> matrix A so that elements below the k-th subdiagonal are zero. The
41: *> reduction is performed by an orthogonal similarity transformation
42: *> Q**T * A * Q. The routine returns the matrices V and T which determine
43: *> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] N
50: *> \verbatim
51: *> N is INTEGER
52: *> The order of the matrix A.
53: *> \endverbatim
54: *>
55: *> \param[in] K
56: *> \verbatim
57: *> K is INTEGER
58: *> The offset for the reduction. Elements below the k-th
59: *> subdiagonal in the first NB columns are reduced to zero.
60: *> \endverbatim
61: *>
62: *> \param[in] NB
63: *> \verbatim
64: *> NB is INTEGER
65: *> The number of columns to be reduced.
66: *> \endverbatim
67: *>
68: *> \param[in,out] A
69: *> \verbatim
70: *> A is DOUBLE PRECISION array, dimension (LDA,N-K+1)
71: *> On entry, the n-by-(n-k+1) general matrix A.
72: *> On exit, the elements on and above the k-th subdiagonal in
73: *> the first NB columns are overwritten with the corresponding
74: *> elements of the reduced matrix; the elements below the k-th
75: *> subdiagonal, with the array TAU, represent the matrix Q as a
76: *> product of elementary reflectors. The other columns of A are
77: *> unchanged. See Further Details.
78: *> \endverbatim
79: *>
80: *> \param[in] LDA
81: *> \verbatim
82: *> LDA is INTEGER
83: *> The leading dimension of the array A. LDA >= max(1,N).
84: *> \endverbatim
85: *>
86: *> \param[out] TAU
87: *> \verbatim
88: *> TAU is DOUBLE PRECISION array, dimension (NB)
89: *> The scalar factors of the elementary reflectors. See Further
90: *> Details.
91: *> \endverbatim
92: *>
93: *> \param[out] T
94: *> \verbatim
95: *> T is DOUBLE PRECISION array, dimension (LDT,NB)
96: *> The upper triangular matrix T.
97: *> \endverbatim
98: *>
99: *> \param[in] LDT
100: *> \verbatim
101: *> LDT is INTEGER
102: *> The leading dimension of the array T. LDT >= NB.
103: *> \endverbatim
104: *>
105: *> \param[out] Y
106: *> \verbatim
107: *> Y is DOUBLE PRECISION array, dimension (LDY,NB)
108: *> The n-by-nb matrix Y.
109: *> \endverbatim
110: *>
111: *> \param[in] LDY
112: *> \verbatim
113: *> LDY is INTEGER
114: *> The leading dimension of the array Y. LDY >= 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 doubleOTHERauxiliary
126: *
127: *> \par Further Details:
128: * =====================
129: *>
130: *> \verbatim
131: *>
132: *> The matrix Q is represented as a product of nb elementary reflectors
133: *>
134: *> Q = H(1) H(2) . . . H(nb).
135: *>
136: *> Each H(i) has the form
137: *>
138: *> H(i) = I - tau * v * v**T
139: *>
140: *> where tau is a real scalar, and v is a real vector with
141: *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
142: *> A(i+k+1:n,i), and tau in TAU(i).
143: *>
144: *> The elements of the vectors v together form the (n-k+1)-by-nb matrix
145: *> V which is needed, with T and Y, to apply the transformation to the
146: *> unreduced part of the matrix, using an update of the form:
147: *> A := (I - V*T*V**T) * (A - Y*V**T).
148: *>
149: *> The contents of A on exit are illustrated by the following example
150: *> with n = 7, k = 3 and nb = 2:
151: *>
152: *> ( a h a a a )
153: *> ( a h a a a )
154: *> ( a h a a a )
155: *> ( h h a a a )
156: *> ( v1 h a a a )
157: *> ( v1 v2 a a a )
158: *> ( v1 v2 a a a )
159: *>
160: *> where a denotes an element of the original matrix A, h denotes a
161: *> modified element of the upper Hessenberg matrix H, and vi denotes an
162: *> element of the vector defining H(i).
163: *> \endverbatim
164: *>
165: * =====================================================================
166: SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
167: *
168: * -- LAPACK auxiliary routine --
169: * -- LAPACK is a software package provided by Univ. of Tennessee, --
170: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171: *
172: * .. Scalar Arguments ..
173: INTEGER K, LDA, LDT, LDY, N, NB
174: * ..
175: * .. Array Arguments ..
176: DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ),
177: $ Y( LDY, NB )
178: * ..
179: *
180: * =====================================================================
181: *
182: * .. Parameters ..
183: DOUBLE PRECISION ZERO, ONE
184: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
185: * ..
186: * .. Local Scalars ..
187: INTEGER I
188: DOUBLE PRECISION EI
189: * ..
190: * .. External Subroutines ..
191: EXTERNAL DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV
192: * ..
193: * .. Intrinsic Functions ..
194: INTRINSIC MIN
195: * ..
196: * .. Executable Statements ..
197: *
198: * Quick return if possible
199: *
200: IF( N.LE.1 )
201: $ RETURN
202: *
203: DO 10 I = 1, NB
204: IF( I.GT.1 ) THEN
205: *
206: * Update A(1:n,i)
207: *
208: * Compute i-th column of A - Y * V**T
209: *
210: CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
211: $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
212: *
213: * Apply I - V * T**T * V**T to this column (call it b) from the
214: * left, using the last column of T as workspace
215: *
216: * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
217: * ( V2 ) ( b2 )
218: *
219: * where V1 is unit lower triangular
220: *
221: * w := V1**T * b1
222: *
223: CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
224: CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
225: $ LDA, T( 1, NB ), 1 )
226: *
227: * w := w + V2**T *b2
228: *
229: CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
230: $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
231: *
232: * w := T**T *w
233: *
234: CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
235: $ T( 1, NB ), 1 )
236: *
237: * b2 := b2 - V2*w
238: *
239: CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
240: $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
241: *
242: * b1 := b1 - V1*w
243: *
244: CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1,
245: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
246: CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
247: *
248: A( K+I-1, I-1 ) = EI
249: END IF
250: *
251: * Generate the elementary reflector H(i) to annihilate
252: * A(k+i+1:n,i)
253: *
254: CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
255: $ TAU( I ) )
256: EI = A( K+I, I )
257: A( K+I, I ) = ONE
258: *
259: * Compute Y(1:n,i)
260: *
261: CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
262: $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
263: CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
264: $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
265: CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
266: $ ONE, Y( 1, I ), 1 )
267: CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 )
268: *
269: * Compute T(1:i,i)
270: *
271: CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
272: CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
273: $ T( 1, I ), 1 )
274: T( I, I ) = TAU( I )
275: *
276: 10 CONTINUE
277: A( K+NB, NB ) = EI
278: *
279: RETURN
280: *
281: * End of DLAHRD
282: *
283: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlahrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack