Mercurial > octave-nkf

comparison libcruft/lapack/clahrd.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE CLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER K, LDA, LDT, LDY, N, NB
9 * ..
10 * .. Array Arguments ..
11 COMPLEX A( LDA, * ), T( LDT, NB ), TAU( NB ),
12 $ Y( LDY, NB )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * CLAHRD reduces the first NB columns of a complex 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 a unitary 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 *
24 * This is an OBSOLETE auxiliary routine.
25 * This routine will be 'deprecated' in a future release.
26 * Please use the new routine CLAHR2 instead.
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) COMPLEX 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) COMPLEX array, dimension (NB)
54 * The scalar factors of the elementary reflectors. See Further
55 * Details.
56 *
57 * T (output) COMPLEX array, dimension (LDT,NB)
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) COMPLEX 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 >= max(1,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 complex scalar, and v is a complex 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 COMPLEX ZERO, ONE
108 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
109 $ ONE = ( 1.0E+0, 0.0E+0 ) )
110 * ..
111 * .. Local Scalars ..
112 INTEGER I
113 COMPLEX EI
114 * ..
115 * .. External Subroutines ..
116 EXTERNAL CAXPY, CCOPY, CGEMV, CLACGV, CLARFG, CSCAL,
117 $ CTRMV
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC MIN
121 * ..
122 * .. Executable Statements ..
123 *
124 * Quick return if possible
125 *
126 IF( N.LE.1 )
127 $ RETURN
128 *
129 DO 10 I = 1, NB
130 IF( I.GT.1 ) THEN
131 *
132 * Update A(1:n,i)
133 *
134 * Compute i-th column of A - Y * V'
135 *
136 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
137 CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
138 $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
139 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
140 *
141 * Apply I - V * T' * V' to this column (call it b) from the
142 * left, using the last column of T as workspace
143 *
144 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
145 * ( V2 ) ( b2 )
146 *
147 * where V1 is unit lower triangular
148 *
149 * w := V1' * b1
150 *
151 CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
152 CALL CTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
153 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
154 *
155 * w := w + V2'*b2
156 *
157 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
158 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
159 $ T( 1, NB ), 1 )
160 *
161 * w := T'*w
162 *
163 CALL CTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
164 $ T, LDT, T( 1, NB ), 1 )
165 *
166 * b2 := b2 - V2*w
167 *
168 CALL CGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
169 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
170 *
171 * b1 := b1 - V1*w
172 *
173 CALL CTRMV( 'Lower', 'No transpose', 'Unit', I-1,
174 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
175 CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
176 *
177 A( K+I-1, I-1 ) = EI
178 END IF
179 *
180 * Generate the elementary reflector H(i) to annihilate
181 * A(k+i+1:n,i)
182 *
183 EI = A( K+I, I )
184 CALL CLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
185 $ TAU( I ) )
186 A( K+I, I ) = ONE
187 *
188 * Compute Y(1:n,i)
189 *
190 CALL CGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
191 $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
192 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
193 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
194 $ 1 )
195 CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
196 $ ONE, Y( 1, I ), 1 )
197 CALL CSCAL( N, TAU( I ), Y( 1, I ), 1 )
198 *
199 * Compute T(1:i,i)
200 *
201 CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
202 CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
203 $ T( 1, I ), 1 )
204 T( I, I ) = TAU( I )
205 *
206 10 CONTINUE
207 A( K+NB, NB ) = EI
208 *
209 RETURN
210 *
211 * End of CLAHRD
212 *
213 END