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