3 /* Table of constant values */
5 static real c_b4 = -1.f;
6 static real c_b5 = 1.f;
7 static integer c__1 = 1;
8 static real c_b16 = 0.f;
10 /* Subroutine */ int slabrd_(integer *m, integer *n, integer *nb, real *a,
11 integer *lda, real *d__, real *e, real *tauq, real *taup, real *x,
12 integer *ldx, real *y, integer *ldy)
14 /* System generated locals */
15 integer a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset, i__1, i__2,
20 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
21 sgemv_(char *, integer *, integer *, real *, real *, integer *,
22 real *, integer *, real *, real *, integer *), slarfg_(
23 integer *, real *, real *, integer *, real *);
26 /* -- LAPACK auxiliary routine (version 3.1) -- */
27 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
30 /* .. Scalar Arguments .. */
32 /* .. Array Arguments .. */
38 /* SLABRD reduces the first NB rows and columns of a real general */
39 /* m by n matrix A to upper or lower bidiagonal form by an orthogonal */
40 /* transformation Q' * A * P, and returns the matrices X and Y which */
41 /* are needed to apply the transformation to the unreduced part of A. */
43 /* If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower */
44 /* bidiagonal form. */
46 /* This is an auxiliary routine called by SGEBRD */
51 /* M (input) INTEGER */
52 /* The number of rows in the matrix A. */
54 /* N (input) INTEGER */
55 /* The number of columns in the matrix A. */
57 /* NB (input) INTEGER */
58 /* The number of leading rows and columns of A to be reduced. */
60 /* A (input/output) REAL array, dimension (LDA,N) */
61 /* On entry, the m by n general matrix to be reduced. */
62 /* On exit, the first NB rows and columns of the matrix are */
63 /* overwritten; the rest of the array is unchanged. */
64 /* If m >= n, elements on and below the diagonal in the first NB */
65 /* columns, with the array TAUQ, represent the orthogonal */
66 /* matrix Q as a product of elementary reflectors; and */
67 /* elements above the diagonal in the first NB rows, with the */
68 /* array TAUP, represent the orthogonal matrix P as a product */
69 /* of elementary reflectors. */
70 /* If m < n, elements below the diagonal in the first NB */
71 /* columns, with the array TAUQ, represent the orthogonal */
72 /* matrix Q as a product of elementary reflectors, and */
73 /* elements on and above the diagonal in the first NB rows, */
74 /* with the array TAUP, represent the orthogonal matrix P as */
75 /* a product of elementary reflectors. */
76 /* See Further Details. */
78 /* LDA (input) INTEGER */
79 /* The leading dimension of the array A. LDA >= max(1,M). */
81 /* D (output) REAL array, dimension (NB) */
82 /* The diagonal elements of the first NB rows and columns of */
83 /* the reduced matrix. D(i) = A(i,i). */
85 /* E (output) REAL array, dimension (NB) */
86 /* The off-diagonal elements of the first NB rows and columns of */
87 /* the reduced matrix. */
89 /* TAUQ (output) REAL array dimension (NB) */
90 /* The scalar factors of the elementary reflectors which */
91 /* represent the orthogonal matrix Q. See Further Details. */
93 /* TAUP (output) REAL array, dimension (NB) */
94 /* The scalar factors of the elementary reflectors which */
95 /* represent the orthogonal matrix P. See Further Details. */
97 /* X (output) REAL array, dimension (LDX,NB) */
98 /* The m-by-nb matrix X required to update the unreduced part */
101 /* LDX (input) INTEGER */
102 /* The leading dimension of the array X. LDX >= M. */
104 /* Y (output) REAL array, dimension (LDY,NB) */
105 /* The n-by-nb matrix Y required to update the unreduced part */
108 /* LDY (input) INTEGER */
109 /* The leading dimension of the array Y. LDY >= N. */
111 /* Further Details */
112 /* =============== */
114 /* The matrices Q and P are represented as products of elementary */
117 /* Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) */
119 /* Each H(i) and G(i) has the form: */
121 /* H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' */
123 /* where tauq and taup are real scalars, and v and u are real vectors. */
125 /* If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in */
126 /* A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in */
127 /* A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). */
129 /* If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in */
130 /* A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in */
131 /* A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). */
133 /* The elements of the vectors v and u together form the m-by-nb matrix */
134 /* V and the nb-by-n matrix U' which are needed, with X and Y, to apply */
135 /* the transformation to the unreduced part of the matrix, using a block */
136 /* update of the form: A := A - V*Y' - X*U'. */
138 /* The contents of A on exit are illustrated by the following examples */
141 /* m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): */
143 /* ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) */
144 /* ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) */
145 /* ( v1 v2 a a a ) ( v1 1 a a a a ) */
146 /* ( v1 v2 a a a ) ( v1 v2 a a a a ) */
147 /* ( v1 v2 a a a ) ( v1 v2 a a a a ) */
148 /* ( v1 v2 a a a ) */
150 /* where a denotes an element of the original matrix which is unchanged, */
151 /* vi denotes an element of the vector defining H(i), and ui an element */
152 /* of the vector defining G(i). */
154 /* ===================================================================== */
156 /* .. Parameters .. */
158 /* .. Local Scalars .. */
160 /* .. External Subroutines .. */
162 /* .. Intrinsic Functions .. */
164 /* .. Executable Statements .. */
166 /* Quick return if possible */
168 /* Parameter adjustments */
170 a_offset = 1 + a_dim1;
177 x_offset = 1 + x_dim1;
180 y_offset = 1 + y_dim1;
184 if (*m <= 0 || *n <= 0) {
190 /* Reduce to upper bidiagonal form */
193 for (i__ = 1; i__ <= i__1; ++i__) {
195 /* Update A(i:m,i) */
199 sgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + a_dim1], lda,
200 &y[i__ + y_dim1], ldy, &c_b5, &a[i__ + i__ * a_dim1], &
204 sgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + x_dim1], ldx,
205 &a[i__ * a_dim1 + 1], &c__1, &c_b5, &a[i__ + i__ *
208 /* Generate reflection Q(i) to annihilate A(i+1:m,i) */
213 slarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[min(i__3, *m)+ i__ *
214 a_dim1], &c__1, &tauq[i__]);
215 d__[i__] = a[i__ + i__ * a_dim1];
217 a[i__ + i__ * a_dim1] = 1.f;
219 /* Compute Y(i+1:n,i) */
223 sgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + (i__ + 1) *
224 a_dim1], lda, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &
225 y[i__ + 1 + i__ * y_dim1], &c__1);
228 sgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + a_dim1],
229 lda, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &y[i__ *
233 sgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + 1 +
234 y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[
235 i__ + 1 + i__ * y_dim1], &c__1);
238 sgemv_("Transpose", &i__2, &i__3, &c_b5, &x[i__ + x_dim1],
239 ldx, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &y[i__ *
243 sgemv_("Transpose", &i__2, &i__3, &c_b4, &a[(i__ + 1) *
244 a_dim1 + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &c_b5,
245 &y[i__ + 1 + i__ * y_dim1], &c__1);
247 sscal_(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
249 /* Update A(i,i+1:n) */
252 sgemv_("No transpose", &i__2, &i__, &c_b4, &y[i__ + 1 +
253 y_dim1], ldy, &a[i__ + a_dim1], lda, &c_b5, &a[i__ + (
254 i__ + 1) * a_dim1], lda);
257 sgemv_("Transpose", &i__2, &i__3, &c_b4, &a[(i__ + 1) *
258 a_dim1 + 1], lda, &x[i__ + x_dim1], ldx, &c_b5, &a[
259 i__ + (i__ + 1) * a_dim1], lda);
261 /* Generate reflection P(i) to annihilate A(i,i+2:n) */
266 slarfg_(&i__2, &a[i__ + (i__ + 1) * a_dim1], &a[i__ + min(
267 i__3, *n)* a_dim1], lda, &taup[i__]);
268 e[i__] = a[i__ + (i__ + 1) * a_dim1];
269 a[i__ + (i__ + 1) * a_dim1] = 1.f;
271 /* Compute X(i+1:m,i) */
275 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + (i__
276 + 1) * a_dim1], lda, &a[i__ + (i__ + 1) * a_dim1],
277 lda, &c_b16, &x[i__ + 1 + i__ * x_dim1], &c__1);
279 sgemv_("Transpose", &i__2, &i__, &c_b5, &y[i__ + 1 + y_dim1],
280 ldy, &a[i__ + (i__ + 1) * a_dim1], lda, &c_b16, &x[
281 i__ * x_dim1 + 1], &c__1);
283 sgemv_("No transpose", &i__2, &i__, &c_b4, &a[i__ + 1 +
284 a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
285 i__ + 1 + i__ * x_dim1], &c__1);
288 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[(i__ + 1) *
289 a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, &
290 c_b16, &x[i__ * x_dim1 + 1], &c__1);
293 sgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + 1 +
294 x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
295 i__ + 1 + i__ * x_dim1], &c__1);
297 sscal_(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
303 /* Reduce to lower bidiagonal form */
306 for (i__ = 1; i__ <= i__1; ++i__) {
308 /* Update A(i,i:n) */
312 sgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + y_dim1], ldy,
313 &a[i__ + a_dim1], lda, &c_b5, &a[i__ + i__ * a_dim1],
317 sgemv_("Transpose", &i__2, &i__3, &c_b4, &a[i__ * a_dim1 + 1],
318 lda, &x[i__ + x_dim1], ldx, &c_b5, &a[i__ + i__ * a_dim1],
321 /* Generate reflection P(i) to annihilate A(i,i+1:n) */
326 slarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[i__ + min(i__3, *n)*
327 a_dim1], lda, &taup[i__]);
328 d__[i__] = a[i__ + i__ * a_dim1];
330 a[i__ + i__ * a_dim1] = 1.f;
332 /* Compute X(i+1:m,i) */
336 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + i__ *
337 a_dim1], lda, &a[i__ + i__ * a_dim1], lda, &c_b16, &
338 x[i__ + 1 + i__ * x_dim1], &c__1);
341 sgemv_("Transpose", &i__2, &i__3, &c_b5, &y[i__ + y_dim1],
342 ldy, &a[i__ + i__ * a_dim1], lda, &c_b16, &x[i__ *
346 sgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + 1 +
347 a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
348 i__ + 1 + i__ * x_dim1], &c__1);
351 sgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ * a_dim1 +
352 1], lda, &a[i__ + i__ * a_dim1], lda, &c_b16, &x[i__ *
356 sgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + 1 +
357 x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
358 i__ + 1 + i__ * x_dim1], &c__1);
360 sscal_(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
362 /* Update A(i+1:m,i) */
366 sgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + 1 +
367 a_dim1], lda, &y[i__ + y_dim1], ldy, &c_b5, &a[i__ +
368 1 + i__ * a_dim1], &c__1);
370 sgemv_("No transpose", &i__2, &i__, &c_b4, &x[i__ + 1 +
371 x_dim1], ldx, &a[i__ * a_dim1 + 1], &c__1, &c_b5, &a[
372 i__ + 1 + i__ * a_dim1], &c__1);
374 /* Generate reflection Q(i) to annihilate A(i+2:m,i) */
379 slarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *m)+
380 i__ * a_dim1], &c__1, &tauq[i__]);
381 e[i__] = a[i__ + 1 + i__ * a_dim1];
382 a[i__ + 1 + i__ * a_dim1] = 1.f;
384 /* Compute Y(i+1:n,i) */
388 sgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + (i__ +
389 1) * a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1,
390 &c_b16, &y[i__ + 1 + i__ * y_dim1], &c__1);
393 sgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + a_dim1],
394 lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &y[
395 i__ * y_dim1 + 1], &c__1);
398 sgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + 1 +
399 y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[
400 i__ + 1 + i__ * y_dim1], &c__1);
402 sgemv_("Transpose", &i__2, &i__, &c_b5, &x[i__ + 1 + x_dim1],
403 ldx, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &y[
404 i__ * y_dim1 + 1], &c__1);
406 sgemv_("Transpose", &i__, &i__2, &c_b4, &a[(i__ + 1) * a_dim1
407 + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[i__
408 + 1 + i__ * y_dim1], &c__1);
410 sscal_(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);