3 /* Table of constant values */
5 static integer c__1 = 1;
6 static real c_b14 = 1.f;
7 static real c_b25 = -1.f;
9 /* Subroutine */ int slarfb_(char *side, char *trans, char *direct, char *
10 storev, integer *m, integer *n, integer *k, real *v, integer *ldv,
11 real *t, integer *ldt, real *c__, integer *ldc, real *work, integer *
14 /* System generated locals */
15 integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1,
16 work_offset, i__1, i__2;
20 extern logical lsame_(char *, char *);
21 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
22 integer *, real *, real *, integer *, real *, integer *, real *,
23 real *, integer *), scopy_(integer *, real *,
24 integer *, real *, integer *), strmm_(char *, char *, char *,
25 char *, integer *, integer *, real *, 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 /* SLARFB applies a real block reflector H or its transpose H' to a */
43 /* real m by n matrix C, from either the left or the right. */
48 /* SIDE (input) CHARACTER*1 */
49 /* = 'L': apply H or H' from the Left */
50 /* = 'R': apply H or H' from the Right */
52 /* TRANS (input) CHARACTER*1 */
53 /* = 'N': apply H (No transpose) */
54 /* = 'T': apply H' (Transpose) */
56 /* DIRECT (input) CHARACTER*1 */
57 /* Indicates how H is formed from a product of elementary */
59 /* = 'F': H = H(1) H(2) . . . H(k) (Forward) */
60 /* = 'B': H = H(k) . . . H(2) H(1) (Backward) */
62 /* STOREV (input) CHARACTER*1 */
63 /* Indicates how the vectors which define the elementary */
64 /* reflectors are stored: */
65 /* = 'C': Columnwise */
68 /* M (input) INTEGER */
69 /* The number of rows of the matrix C. */
71 /* N (input) INTEGER */
72 /* The number of columns of the matrix C. */
74 /* K (input) INTEGER */
75 /* The order of the matrix T (= the number of elementary */
76 /* reflectors whose product defines the block reflector). */
78 /* V (input) REAL array, dimension */
79 /* (LDV,K) if STOREV = 'C' */
80 /* (LDV,M) if STOREV = 'R' and SIDE = 'L' */
81 /* (LDV,N) if STOREV = 'R' and SIDE = 'R' */
82 /* The matrix V. See further details. */
84 /* LDV (input) INTEGER */
85 /* The leading dimension of the array V. */
86 /* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); */
87 /* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); */
88 /* if STOREV = 'R', LDV >= K. */
90 /* T (input) REAL array, dimension (LDT,K) */
91 /* The triangular k by k matrix T in the representation of the */
92 /* block reflector. */
94 /* LDT (input) INTEGER */
95 /* The leading dimension of the array T. LDT >= K. */
97 /* C (input/output) REAL array, dimension (LDC,N) */
98 /* On entry, the m by n matrix C. */
99 /* On exit, C is overwritten by H*C or H'*C or C*H or C*H'. */
101 /* LDC (input) INTEGER */
102 /* The leading dimension of the array C. LDA >= max(1,M). */
104 /* WORK (workspace) REAL array, dimension (LDWORK,K) */
106 /* LDWORK (input) INTEGER */
107 /* The leading dimension of the array WORK. */
108 /* If SIDE = 'L', LDWORK >= max(1,N); */
109 /* if SIDE = 'R', LDWORK >= max(1,M). */
111 /* ===================================================================== */
113 /* .. Parameters .. */
115 /* .. Local Scalars .. */
117 /* .. External Functions .. */
119 /* .. External Subroutines .. */
121 /* .. Executable Statements .. */
123 /* Quick return if possible */
125 /* Parameter adjustments */
127 v_offset = 1 + v_dim1;
130 t_offset = 1 + t_dim1;
133 c_offset = 1 + c_dim1;
136 work_offset = 1 + work_dim1;
140 if (*m <= 0 || *n <= 0) {
144 if (lsame_(trans, "N")) {
145 *(unsigned char *)transt = 'T';
147 *(unsigned char *)transt = 'N';
150 if (lsame_(storev, "C")) {
152 if (lsame_(direct, "F")) {
154 /* Let V = ( V1 ) (first K rows) */
156 /* where V1 is unit lower triangular. */
158 if (lsame_(side, "L")) {
160 /* Form H * C or H' * C where C = ( C1 ) */
163 /* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) */
168 for (j = 1; j <= i__1; ++j) {
169 scopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
176 strmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
177 &v[v_offset], ldv, &work[work_offset], ldwork);
180 /* W := W + C2'*V2 */
183 sgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
184 c__[*k + 1 + c_dim1], ldc, &v[*k + 1 + v_dim1],
185 ldv, &c_b14, &work[work_offset], ldwork);
188 /* W := W * T' or W * T */
190 strmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
191 t_offset], ldt, &work[work_offset], ldwork);
193 /* C := C - V * W' */
197 /* C2 := C2 - V2 * W' */
200 sgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
201 v[*k + 1 + v_dim1], ldv, &work[work_offset],
202 ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
207 strmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
208 v[v_offset], ldv, &work[work_offset], ldwork);
213 for (j = 1; j <= i__1; ++j) {
215 for (i__ = 1; i__ <= i__2; ++i__) {
216 c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
222 } else if (lsame_(side, "R")) {
224 /* Form C * H or C * H' where C = ( C1 C2 ) */
226 /* W := C * V = (C1*V1 + C2*V2) (stored in WORK) */
231 for (j = 1; j <= i__1; ++j) {
232 scopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
233 work_dim1 + 1], &c__1);
239 strmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
240 &v[v_offset], ldv, &work[work_offset], ldwork);
243 /* W := W + C2 * V2 */
246 sgemm_("No transpose", "No transpose", m, k, &i__1, &
247 c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k +
248 1 + v_dim1], ldv, &c_b14, &work[work_offset],
252 /* W := W * T or W * T' */
254 strmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
255 t_offset], ldt, &work[work_offset], ldwork);
257 /* C := C - W * V' */
261 /* C2 := C2 - W * V2' */
264 sgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
265 work[work_offset], ldwork, &v[*k + 1 + v_dim1],
266 ldv, &c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc);
271 strmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
272 v[v_offset], ldv, &work[work_offset], ldwork);
277 for (j = 1; j <= i__1; ++j) {
279 for (i__ = 1; i__ <= i__2; ++i__) {
280 c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
290 /* ( V2 ) (last K rows) */
291 /* where V2 is unit upper triangular. */
293 if (lsame_(side, "L")) {
295 /* Form H * C or H' * C where C = ( C1 ) */
298 /* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) */
303 for (j = 1; j <= i__1; ++j) {
304 scopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
305 work_dim1 + 1], &c__1);
311 strmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
312 &v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
316 /* W := W + C1'*V1 */
319 sgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
320 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
321 work[work_offset], ldwork);
324 /* W := W * T' or W * T */
326 strmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
327 t_offset], ldt, &work[work_offset], ldwork);
329 /* C := C - V * W' */
333 /* C1 := C1 - V1 * W' */
336 sgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
337 v[v_offset], ldv, &work[work_offset], ldwork, &
338 c_b14, &c__[c_offset], ldc)
344 strmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
345 v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
351 for (j = 1; j <= i__1; ++j) {
353 for (i__ = 1; i__ <= i__2; ++i__) {
354 c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
361 } else if (lsame_(side, "R")) {
363 /* Form C * H or C * H' where C = ( C1 C2 ) */
365 /* W := C * V = (C1*V1 + C2*V2) (stored in WORK) */
370 for (j = 1; j <= i__1; ++j) {
371 scopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
372 j * work_dim1 + 1], &c__1);
378 strmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
379 &v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
383 /* W := W + C1 * V1 */
386 sgemm_("No transpose", "No transpose", m, k, &i__1, &
387 c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, &
388 c_b14, &work[work_offset], ldwork);
391 /* W := W * T or W * T' */
393 strmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
394 t_offset], ldt, &work[work_offset], ldwork);
396 /* C := C - W * V' */
400 /* C1 := C1 - W * V1' */
403 sgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
404 work[work_offset], ldwork, &v[v_offset], ldv, &
405 c_b14, &c__[c_offset], ldc)
411 strmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
412 v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
418 for (j = 1; j <= i__1; ++j) {
420 for (i__ = 1; i__ <= i__2; ++i__) {
421 c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
430 } else if (lsame_(storev, "R")) {
432 if (lsame_(direct, "F")) {
434 /* Let V = ( V1 V2 ) (V1: first K columns) */
435 /* where V1 is unit upper triangular. */
437 if (lsame_(side, "L")) {
439 /* Form H * C or H' * C where C = ( C1 ) */
442 /* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
447 for (j = 1; j <= i__1; ++j) {
448 scopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
455 strmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
456 v[v_offset], ldv, &work[work_offset], ldwork);
459 /* W := W + C2'*V2' */
462 sgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
463 c__[*k + 1 + c_dim1], ldc, &v[(*k + 1) * v_dim1 +
464 1], ldv, &c_b14, &work[work_offset], ldwork);
467 /* W := W * T' or W * T */
469 strmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
470 t_offset], ldt, &work[work_offset], ldwork);
472 /* C := C - V' * W' */
476 /* C2 := C2 - V2' * W' */
479 sgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[(
480 *k + 1) * v_dim1 + 1], ldv, &work[work_offset],
481 ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
486 strmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
487 &v[v_offset], ldv, &work[work_offset], ldwork);
492 for (j = 1; j <= i__1; ++j) {
494 for (i__ = 1; i__ <= i__2; ++i__) {
495 c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
501 } else if (lsame_(side, "R")) {
503 /* Form C * H or C * H' where C = ( C1 C2 ) */
505 /* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
510 for (j = 1; j <= i__1; ++j) {
511 scopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
512 work_dim1 + 1], &c__1);
518 strmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
519 v[v_offset], ldv, &work[work_offset], ldwork);
522 /* W := W + C2 * V2' */
525 sgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
526 c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) *
527 v_dim1 + 1], ldv, &c_b14, &work[work_offset],
531 /* W := W * T or W * T' */
533 strmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
534 t_offset], ldt, &work[work_offset], ldwork);
540 /* C2 := C2 - W * V2 */
543 sgemm_("No transpose", "No transpose", m, &i__1, k, &
544 c_b25, &work[work_offset], ldwork, &v[(*k + 1) *
545 v_dim1 + 1], ldv, &c_b14, &c__[(*k + 1) * c_dim1
551 strmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
552 &v[v_offset], ldv, &work[work_offset], ldwork);
557 for (j = 1; j <= i__1; ++j) {
559 for (i__ = 1; i__ <= i__2; ++i__) {
560 c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
570 /* Let V = ( V1 V2 ) (V2: last K columns) */
571 /* where V2 is unit lower triangular. */
573 if (lsame_(side, "L")) {
575 /* Form H * C or H' * C where C = ( C1 ) */
578 /* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
583 for (j = 1; j <= i__1; ++j) {
584 scopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
585 work_dim1 + 1], &c__1);
591 strmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
592 v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
596 /* W := W + C1'*V1' */
599 sgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
600 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
601 work[work_offset], ldwork);
604 /* W := W * T' or W * T */
606 strmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
607 t_offset], ldt, &work[work_offset], ldwork);
609 /* C := C - V' * W' */
613 /* C1 := C1 - V1' * W' */
616 sgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
617 v_offset], ldv, &work[work_offset], ldwork, &
618 c_b14, &c__[c_offset], ldc);
623 strmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
624 &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
625 work_offset], ldwork);
630 for (j = 1; j <= i__1; ++j) {
632 for (i__ = 1; i__ <= i__2; ++i__) {
633 c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
640 } else if (lsame_(side, "R")) {
642 /* Form C * H or C * H' where C = ( C1 C2 ) */
644 /* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
649 for (j = 1; j <= i__1; ++j) {
650 scopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
651 j * work_dim1 + 1], &c__1);
657 strmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
658 v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
662 /* W := W + C1 * V1' */
665 sgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
666 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
667 work[work_offset], ldwork);
670 /* W := W * T or W * T' */
672 strmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
673 t_offset], ldt, &work[work_offset], ldwork);
679 /* C1 := C1 - W * V1 */
682 sgemm_("No transpose", "No transpose", m, &i__1, k, &
683 c_b25, &work[work_offset], ldwork, &v[v_offset],
684 ldv, &c_b14, &c__[c_offset], ldc);
689 strmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
690 &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
691 work_offset], ldwork);
696 for (j = 1; j <= i__1; ++j) {
698 for (i__ = 1; i__ <= i__2; ++i__) {
699 c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *