Annotation of rpl/lapack/lapack/dlahrd.f, revision 1.15
1.12 bertrand 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.
1.9 bertrand 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: *>
1.15 ! bertrand 37: *> This routine is deprecated and has been replaced by routine DLAHR2.
! 38: *>
1.9 bertrand 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: *
1.15 ! bertrand 125: *> \date November 2015
1.9 bertrand 126: *
127: *> \ingroup doubleOTHERauxiliary
128: *
129: *> \par Further Details:
130: * =====================
131: *>
132: *> \verbatim
133: *>
134: *> The matrix Q is represented as a product of nb elementary reflectors
135: *>
136: *> Q = H(1) H(2) . . . H(nb).
137: *>
138: *> Each H(i) has the form
139: *>
140: *> H(i) = I - tau * v * v**T
141: *>
142: *> where tau is a real scalar, and v is a real vector with
143: *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
144: *> A(i+k+1:n,i), and tau in TAU(i).
145: *>
146: *> The elements of the vectors v together form the (n-k+1)-by-nb matrix
147: *> V which is needed, with T and Y, to apply the transformation to the
148: *> unreduced part of the matrix, using an update of the form:
149: *> A := (I - V*T*V**T) * (A - Y*V**T).
150: *>
151: *> The contents of A on exit are illustrated by the following example
152: *> with n = 7, k = 3 and nb = 2:
153: *>
154: *> ( a h a a a )
155: *> ( a h a a a )
156: *> ( a h a a a )
157: *> ( h h a a a )
158: *> ( v1 h a a a )
159: *> ( v1 v2 a a a )
160: *> ( v1 v2 a a a )
161: *>
162: *> where a denotes an element of the original matrix A, h denotes a
163: *> modified element of the upper Hessenberg matrix H, and vi denotes an
164: *> element of the vector defining H(i).
165: *> \endverbatim
166: *>
167: * =====================================================================
1.1 bertrand 168: SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
169: *
1.15 ! bertrand 170: * -- LAPACK auxiliary routine (version 3.6.0) --
1.1 bertrand 171: * -- LAPACK is a software package provided by Univ. of Tennessee, --
172: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 ! bertrand 173: * November 2015
1.1 bertrand 174: *
175: * .. Scalar Arguments ..
176: INTEGER K, LDA, LDT, LDY, N, NB
177: * ..
178: * .. Array Arguments ..
179: DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ),
180: $ Y( LDY, NB )
181: * ..
182: *
183: * =====================================================================
184: *
185: * .. Parameters ..
186: DOUBLE PRECISION ZERO, ONE
187: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
188: * ..
189: * .. Local Scalars ..
190: INTEGER I
191: DOUBLE PRECISION EI
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV
195: * ..
196: * .. Intrinsic Functions ..
197: INTRINSIC MIN
198: * ..
199: * .. Executable Statements ..
200: *
201: * Quick return if possible
202: *
203: IF( N.LE.1 )
204: $ RETURN
205: *
206: DO 10 I = 1, NB
207: IF( I.GT.1 ) THEN
208: *
209: * Update A(1:n,i)
210: *
1.8 bertrand 211: * Compute i-th column of A - Y * V**T
1.1 bertrand 212: *
213: CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
214: $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
215: *
1.8 bertrand 216: * Apply I - V * T**T * V**T to this column (call it b) from the
1.1 bertrand 217: * left, using the last column of T as workspace
218: *
219: * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
220: * ( V2 ) ( b2 )
221: *
222: * where V1 is unit lower triangular
223: *
1.8 bertrand 224: * w := V1**T * b1
1.1 bertrand 225: *
226: CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
227: CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
228: $ LDA, T( 1, NB ), 1 )
229: *
1.8 bertrand 230: * w := w + V2**T *b2
1.1 bertrand 231: *
232: CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
233: $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
234: *
1.8 bertrand 235: * w := T**T *w
1.1 bertrand 236: *
237: CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
238: $ T( 1, NB ), 1 )
239: *
240: * b2 := b2 - V2*w
241: *
242: CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
243: $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
244: *
245: * b1 := b1 - V1*w
246: *
247: CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1,
248: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
249: CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
250: *
251: A( K+I-1, I-1 ) = EI
252: END IF
253: *
254: * Generate the elementary reflector H(i) to annihilate
255: * A(k+i+1:n,i)
256: *
257: CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
258: $ TAU( I ) )
259: EI = A( K+I, I )
260: A( K+I, I ) = ONE
261: *
262: * Compute Y(1:n,i)
263: *
264: CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
265: $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
266: CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
267: $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
268: CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
269: $ ONE, Y( 1, I ), 1 )
270: CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 )
271: *
272: * Compute T(1:i,i)
273: *
274: CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
275: CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
276: $ T( 1, I ), 1 )
277: T( I, I ) = TAU( I )
278: *
279: 10 CONTINUE
280: A( K+NB, NB ) = EI
281: *
282: RETURN
283: *
284: * End of DLAHRD
285: *
286: 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