3 /* Table of constant values */
5 static real c_b5 = -1.f;
6 static real c_b6 = 1.f;
7 static integer c__1 = 1;
8 static real c_b16 = 0.f;
10 /* Subroutine */ int slatrd_(char *uplo, integer *n, integer *nb, real *a,
11 integer *lda, real *e, real *tau, real *w, integer *ldw)
13 /* System generated locals */
14 integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;
18 extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
20 extern logical lsame_(char *, char *);
21 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
22 sgemv_(char *, integer *, integer *, real *, real *, integer *,
23 real *, integer *, real *, real *, integer *), saxpy_(
24 integer *, real *, real *, integer *, real *, integer *), ssymv_(
25 char *, integer *, real *, real *, integer *, real *, integer *,
26 real *, real *, integer *), slarfg_(integer *, real *,
27 real *, integer *, real *);
30 /* -- LAPACK auxiliary routine (version 3.1) -- */
31 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
34 /* .. Scalar Arguments .. */
36 /* .. Array Arguments .. */
42 /* SLATRD reduces NB rows and columns of a real symmetric matrix A to */
43 /* symmetric tridiagonal form by an orthogonal similarity */
44 /* transformation Q' * A * Q, and returns the matrices V and W which are */
45 /* needed to apply the transformation to the unreduced part of A. */
47 /* If UPLO = 'U', SLATRD reduces the last NB rows and columns of a */
48 /* matrix, of which the upper triangle is supplied; */
49 /* if UPLO = 'L', SLATRD reduces the first NB rows and columns of a */
50 /* matrix, of which the lower triangle is supplied. */
52 /* This is an auxiliary routine called by SSYTRD. */
57 /* UPLO (input) CHARACTER*1 */
58 /* Specifies whether the upper or lower triangular part of the */
59 /* symmetric matrix A is stored: */
60 /* = 'U': Upper triangular */
61 /* = 'L': Lower triangular */
63 /* N (input) INTEGER */
64 /* The order of the matrix A. */
66 /* NB (input) INTEGER */
67 /* The number of rows and columns to be reduced. */
69 /* A (input/output) REAL array, dimension (LDA,N) */
70 /* On entry, the symmetric matrix A. If UPLO = 'U', the leading */
71 /* n-by-n upper triangular part of A contains the upper */
72 /* triangular part of the matrix A, and the strictly lower */
73 /* triangular part of A is not referenced. If UPLO = 'L', the */
74 /* leading n-by-n lower triangular part of A contains the lower */
75 /* triangular part of the matrix A, and the strictly upper */
76 /* triangular part of A is not referenced. */
78 /* if UPLO = 'U', the last NB columns have been reduced to */
79 /* tridiagonal form, with the diagonal elements overwriting */
80 /* the diagonal elements of A; the elements above the diagonal */
81 /* with the array TAU, represent the orthogonal matrix Q as a */
82 /* product of elementary reflectors; */
83 /* if UPLO = 'L', the first NB columns have been reduced to */
84 /* tridiagonal form, with the diagonal elements overwriting */
85 /* the diagonal elements of A; the elements below the diagonal */
86 /* with the array TAU, represent the orthogonal matrix Q as a */
87 /* product of elementary reflectors. */
88 /* See Further Details. */
90 /* LDA (input) INTEGER */
91 /* The leading dimension of the array A. LDA >= (1,N). */
93 /* E (output) REAL array, dimension (N-1) */
94 /* If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal */
95 /* elements of the last NB columns of the reduced matrix; */
96 /* if UPLO = 'L', E(1:nb) contains the subdiagonal elements of */
97 /* the first NB columns of the reduced matrix. */
99 /* TAU (output) REAL array, dimension (N-1) */
100 /* The scalar factors of the elementary reflectors, stored in */
101 /* TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. */
102 /* See Further Details. */
104 /* W (output) REAL array, dimension (LDW,NB) */
105 /* The n-by-nb matrix W required to update the unreduced part */
108 /* LDW (input) INTEGER */
109 /* The leading dimension of the array W. LDW >= max(1,N). */
111 /* Further Details */
112 /* =============== */
114 /* If UPLO = 'U', the matrix Q is represented as a product of elementary */
117 /* Q = H(n) H(n-1) . . . H(n-nb+1). */
119 /* Each H(i) has the form */
121 /* H(i) = I - tau * v * v' */
123 /* where tau is a real scalar, and v is a real vector with */
124 /* v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), */
125 /* and tau in TAU(i-1). */
127 /* If UPLO = 'L', the matrix Q is represented as a product of elementary */
130 /* Q = H(1) H(2) . . . H(nb). */
132 /* Each H(i) has the form */
134 /* H(i) = I - tau * v * v' */
136 /* where tau is a real scalar, and v is a real vector with */
137 /* v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), */
138 /* and tau in TAU(i). */
140 /* The elements of the vectors v together form the n-by-nb matrix V */
141 /* which is needed, with W, to apply the transformation to the unreduced */
142 /* part of the matrix, using a symmetric rank-2k update of the form: */
143 /* A := A - V*W' - W*V'. */
145 /* The contents of A on exit are illustrated by the following examples */
146 /* with n = 5 and nb = 2: */
148 /* if UPLO = 'U': if UPLO = 'L': */
150 /* ( a a a v4 v5 ) ( d ) */
151 /* ( a a v4 v5 ) ( 1 d ) */
152 /* ( a 1 v5 ) ( v1 1 a ) */
153 /* ( d 1 ) ( v1 v2 a a ) */
154 /* ( d ) ( v1 v2 a a a ) */
156 /* where d denotes a diagonal element of the reduced matrix, a denotes */
157 /* an element of the original matrix that is unchanged, and vi denotes */
158 /* an element of the vector defining H(i). */
160 /* ===================================================================== */
162 /* .. Parameters .. */
164 /* .. Local Scalars .. */
166 /* .. External Subroutines .. */
168 /* .. External Functions .. */
170 /* .. Intrinsic Functions .. */
172 /* .. Executable Statements .. */
174 /* Quick return if possible */
176 /* Parameter adjustments */
178 a_offset = 1 + a_dim1;
183 w_offset = 1 + w_dim1;
191 if (lsame_(uplo, "U")) {
193 /* Reduce last NB columns of upper triangle */
196 for (i__ = *n; i__ >= i__1; --i__) {
200 /* Update A(1:i,i) */
203 sgemv_("No transpose", &i__, &i__2, &c_b5, &a[(i__ + 1) *
204 a_dim1 + 1], lda, &w[i__ + (iw + 1) * w_dim1], ldw, &
205 c_b6, &a[i__ * a_dim1 + 1], &c__1);
207 sgemv_("No transpose", &i__, &i__2, &c_b5, &w[(iw + 1) *
208 w_dim1 + 1], ldw, &a[i__ + (i__ + 1) * a_dim1], lda, &
209 c_b6, &a[i__ * a_dim1 + 1], &c__1);
213 /* Generate elementary reflector H(i) to annihilate */
217 slarfg_(&i__2, &a[i__ - 1 + i__ * a_dim1], &a[i__ * a_dim1 +
218 1], &c__1, &tau[i__ - 1]);
219 e[i__ - 1] = a[i__ - 1 + i__ * a_dim1];
220 a[i__ - 1 + i__ * a_dim1] = 1.f;
222 /* Compute W(1:i-1,i) */
225 ssymv_("Upper", &i__2, &c_b6, &a[a_offset], lda, &a[i__ *
226 a_dim1 + 1], &c__1, &c_b16, &w[iw * w_dim1 + 1], &
231 sgemv_("Transpose", &i__2, &i__3, &c_b6, &w[(iw + 1) *
232 w_dim1 + 1], ldw, &a[i__ * a_dim1 + 1], &c__1, &
233 c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);
236 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[(i__ + 1) *
237 a_dim1 + 1], lda, &w[i__ + 1 + iw * w_dim1], &
238 c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);
241 sgemv_("Transpose", &i__2, &i__3, &c_b6, &a[(i__ + 1) *
242 a_dim1 + 1], lda, &a[i__ * a_dim1 + 1], &c__1, &
243 c_b16, &w[i__ + 1 + iw * w_dim1], &c__1);
246 sgemv_("No transpose", &i__2, &i__3, &c_b5, &w[(iw + 1) *
247 w_dim1 + 1], ldw, &w[i__ + 1 + iw * w_dim1], &
248 c__1, &c_b6, &w[iw * w_dim1 + 1], &c__1);
251 sscal_(&i__2, &tau[i__ - 1], &w[iw * w_dim1 + 1], &c__1);
253 alpha = tau[i__ - 1] * -.5f * sdot_(&i__2, &w[iw * w_dim1 + 1]
254 , &c__1, &a[i__ * a_dim1 + 1], &c__1);
256 saxpy_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &w[iw *
264 /* Reduce first NB columns of lower triangle */
267 for (i__ = 1; i__ <= i__1; ++i__) {
269 /* Update A(i:n,i) */
273 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + a_dim1], lda,
274 &w[i__ + w_dim1], ldw, &c_b6, &a[i__ + i__ * a_dim1], &
278 sgemv_("No transpose", &i__2, &i__3, &c_b5, &w[i__ + w_dim1], ldw,
279 &a[i__ + a_dim1], lda, &c_b6, &a[i__ + i__ * a_dim1], &
283 /* Generate elementary reflector H(i) to annihilate */
289 slarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *n)+
290 i__ * a_dim1], &c__1, &tau[i__]);
291 e[i__] = a[i__ + 1 + i__ * a_dim1];
292 a[i__ + 1 + i__ * a_dim1] = 1.f;
294 /* Compute W(i+1:n,i) */
297 ssymv_("Lower", &i__2, &c_b6, &a[i__ + 1 + (i__ + 1) * a_dim1]
298 , lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
299 i__ + 1 + i__ * w_dim1], &c__1);
302 sgemv_("Transpose", &i__2, &i__3, &c_b6, &w[i__ + 1 + w_dim1],
303 ldw, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
304 i__ * w_dim1 + 1], &c__1);
307 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 +
308 a_dim1], lda, &w[i__ * w_dim1 + 1], &c__1, &c_b6, &w[
309 i__ + 1 + i__ * w_dim1], &c__1);
312 sgemv_("Transpose", &i__2, &i__3, &c_b6, &a[i__ + 1 + a_dim1],
313 lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &w[
314 i__ * w_dim1 + 1], &c__1);
317 sgemv_("No transpose", &i__2, &i__3, &c_b5, &w[i__ + 1 +
318 w_dim1], ldw, &w[i__ * w_dim1 + 1], &c__1, &c_b6, &w[
319 i__ + 1 + i__ * w_dim1], &c__1);
321 sscal_(&i__2, &tau[i__], &w[i__ + 1 + i__ * w_dim1], &c__1);
323 alpha = tau[i__] * -.5f * sdot_(&i__2, &w[i__ + 1 + i__ *
324 w_dim1], &c__1, &a[i__ + 1 + i__ * a_dim1], &c__1);
326 saxpy_(&i__2, &alpha, &a[i__ + 1 + i__ * a_dim1], &c__1, &w[
327 i__ + 1 + i__ * w_dim1], &c__1);