3 /* Table of constant values */
5 static integer c__1 = 1;
6 static doublereal c_b14 = 1.;
7 static doublereal c_b25 = -1.;
9 /* Subroutine */ int dlarfb_(char *side, char *trans, char *direct, char *
10 storev, integer *m, integer *n, integer *k, doublereal *v, integer *
11 ldv, doublereal *t, integer *ldt, doublereal *c__, integer *ldc,
12 doublereal *work, integer *ldwork)
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 /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
21 integer *, doublereal *, doublereal *, integer *, doublereal *,
22 integer *, doublereal *, doublereal *, integer *);
23 extern logical lsame_(char *, char *);
24 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
25 doublereal *, integer *), dtrmm_(char *, char *, char *, char *,
26 integer *, integer *, doublereal *, doublereal *, integer *,
27 doublereal *, integer *);
31 /* -- LAPACK auxiliary routine (version 3.1) -- */
32 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
35 /* .. Scalar Arguments .. */
37 /* .. Array Arguments .. */
43 /* DLARFB applies a real block reflector H or its transpose H' to a */
44 /* real m by n matrix C, from either the left or the right. */
49 /* SIDE (input) CHARACTER*1 */
50 /* = 'L': apply H or H' from the Left */
51 /* = 'R': apply H or H' from the Right */
53 /* TRANS (input) CHARACTER*1 */
54 /* = 'N': apply H (No transpose) */
55 /* = 'T': apply H' (Transpose) */
57 /* DIRECT (input) CHARACTER*1 */
58 /* Indicates how H is formed from a product of elementary */
60 /* = 'F': H = H(1) H(2) . . . H(k) (Forward) */
61 /* = 'B': H = H(k) . . . H(2) H(1) (Backward) */
63 /* STOREV (input) CHARACTER*1 */
64 /* Indicates how the vectors which define the elementary */
65 /* reflectors are stored: */
66 /* = 'C': Columnwise */
69 /* M (input) INTEGER */
70 /* The number of rows of the matrix C. */
72 /* N (input) INTEGER */
73 /* The number of columns of the matrix C. */
75 /* K (input) INTEGER */
76 /* The order of the matrix T (= the number of elementary */
77 /* reflectors whose product defines the block reflector). */
79 /* V (input) DOUBLE PRECISION array, dimension */
80 /* (LDV,K) if STOREV = 'C' */
81 /* (LDV,M) if STOREV = 'R' and SIDE = 'L' */
82 /* (LDV,N) if STOREV = 'R' and SIDE = 'R' */
83 /* The matrix V. See further details. */
85 /* LDV (input) INTEGER */
86 /* The leading dimension of the array V. */
87 /* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); */
88 /* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); */
89 /* if STOREV = 'R', LDV >= K. */
91 /* T (input) DOUBLE PRECISION array, dimension (LDT,K) */
92 /* The triangular k by k matrix T in the representation of the */
93 /* block reflector. */
95 /* LDT (input) INTEGER */
96 /* The leading dimension of the array T. LDT >= K. */
98 /* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
99 /* On entry, the m by n matrix C. */
100 /* On exit, C is overwritten by H*C or H'*C or C*H or C*H'. */
102 /* LDC (input) INTEGER */
103 /* The leading dimension of the array C. LDA >= max(1,M). */
105 /* WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) */
107 /* LDWORK (input) INTEGER */
108 /* The leading dimension of the array WORK. */
109 /* If SIDE = 'L', LDWORK >= max(1,N); */
110 /* if SIDE = 'R', LDWORK >= max(1,M). */
112 /* ===================================================================== */
114 /* .. Parameters .. */
116 /* .. Local Scalars .. */
118 /* .. External Functions .. */
120 /* .. External Subroutines .. */
122 /* .. Executable Statements .. */
124 /* Quick return if possible */
126 /* Parameter adjustments */
128 v_offset = 1 + v_dim1;
131 t_offset = 1 + t_dim1;
134 c_offset = 1 + c_dim1;
137 work_offset = 1 + work_dim1;
141 if (*m <= 0 || *n <= 0) {
145 if (lsame_(trans, "N")) {
146 *(unsigned char *)transt = 'T';
148 *(unsigned char *)transt = 'N';
151 if (lsame_(storev, "C")) {
153 if (lsame_(direct, "F")) {
155 /* Let V = ( V1 ) (first K rows) */
157 /* where V1 is unit lower triangular. */
159 if (lsame_(side, "L")) {
161 /* Form H * C or H' * C where C = ( C1 ) */
164 /* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) */
169 for (j = 1; j <= i__1; ++j) {
170 dcopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
177 dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
178 &v[v_offset], ldv, &work[work_offset], ldwork);
181 /* W := W + C2'*V2 */
184 dgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
185 c__[*k + 1 + c_dim1], ldc, &v[*k + 1 + v_dim1],
186 ldv, &c_b14, &work[work_offset], ldwork);
189 /* W := W * T' or W * T */
191 dtrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
192 t_offset], ldt, &work[work_offset], ldwork);
194 /* C := C - V * W' */
198 /* C2 := C2 - V2 * W' */
201 dgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
202 v[*k + 1 + v_dim1], ldv, &work[work_offset],
203 ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
208 dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
209 v[v_offset], ldv, &work[work_offset], ldwork);
214 for (j = 1; j <= i__1; ++j) {
216 for (i__ = 1; i__ <= i__2; ++i__) {
217 c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
223 } else if (lsame_(side, "R")) {
225 /* Form C * H or C * H' where C = ( C1 C2 ) */
227 /* W := C * V = (C1*V1 + C2*V2) (stored in WORK) */
232 for (j = 1; j <= i__1; ++j) {
233 dcopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
234 work_dim1 + 1], &c__1);
240 dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
241 &v[v_offset], ldv, &work[work_offset], ldwork);
244 /* W := W + C2 * V2 */
247 dgemm_("No transpose", "No transpose", m, k, &i__1, &
248 c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k +
249 1 + v_dim1], ldv, &c_b14, &work[work_offset],
253 /* W := W * T or W * T' */
255 dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
256 t_offset], ldt, &work[work_offset], ldwork);
258 /* C := C - W * V' */
262 /* C2 := C2 - W * V2' */
265 dgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
266 work[work_offset], ldwork, &v[*k + 1 + v_dim1],
267 ldv, &c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc);
272 dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
273 v[v_offset], ldv, &work[work_offset], ldwork);
278 for (j = 1; j <= i__1; ++j) {
280 for (i__ = 1; i__ <= i__2; ++i__) {
281 c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
291 /* ( V2 ) (last K rows) */
292 /* where V2 is unit upper triangular. */
294 if (lsame_(side, "L")) {
296 /* Form H * C or H' * C where C = ( C1 ) */
299 /* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) */
304 for (j = 1; j <= i__1; ++j) {
305 dcopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
306 work_dim1 + 1], &c__1);
312 dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
313 &v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
317 /* W := W + C1'*V1 */
320 dgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
321 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
322 work[work_offset], ldwork);
325 /* W := W * T' or W * T */
327 dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
328 t_offset], ldt, &work[work_offset], ldwork);
330 /* C := C - V * W' */
334 /* C1 := C1 - V1 * W' */
337 dgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
338 v[v_offset], ldv, &work[work_offset], ldwork, &
339 c_b14, &c__[c_offset], ldc)
345 dtrmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
346 v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
352 for (j = 1; j <= i__1; ++j) {
354 for (i__ = 1; i__ <= i__2; ++i__) {
355 c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
362 } else if (lsame_(side, "R")) {
364 /* Form C * H or C * H' where C = ( C1 C2 ) */
366 /* W := C * V = (C1*V1 + C2*V2) (stored in WORK) */
371 for (j = 1; j <= i__1; ++j) {
372 dcopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
373 j * work_dim1 + 1], &c__1);
379 dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
380 &v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
384 /* W := W + C1 * V1 */
387 dgemm_("No transpose", "No transpose", m, k, &i__1, &
388 c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, &
389 c_b14, &work[work_offset], ldwork);
392 /* W := W * T or W * T' */
394 dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
395 t_offset], ldt, &work[work_offset], ldwork);
397 /* C := C - W * V' */
401 /* C1 := C1 - W * V1' */
404 dgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
405 work[work_offset], ldwork, &v[v_offset], ldv, &
406 c_b14, &c__[c_offset], ldc)
412 dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
413 v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
419 for (j = 1; j <= i__1; ++j) {
421 for (i__ = 1; i__ <= i__2; ++i__) {
422 c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
431 } else if (lsame_(storev, "R")) {
433 if (lsame_(direct, "F")) {
435 /* Let V = ( V1 V2 ) (V1: first K columns) */
436 /* where V1 is unit upper triangular. */
438 if (lsame_(side, "L")) {
440 /* Form H * C or H' * C where C = ( C1 ) */
443 /* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
448 for (j = 1; j <= i__1; ++j) {
449 dcopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
456 dtrmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
457 v[v_offset], ldv, &work[work_offset], ldwork);
460 /* W := W + C2'*V2' */
463 dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
464 c__[*k + 1 + c_dim1], ldc, &v[(*k + 1) * v_dim1 +
465 1], ldv, &c_b14, &work[work_offset], ldwork);
468 /* W := W * T' or W * T */
470 dtrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
471 t_offset], ldt, &work[work_offset], ldwork);
473 /* C := C - V' * W' */
477 /* C2 := C2 - V2' * W' */
480 dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[(
481 *k + 1) * v_dim1 + 1], ldv, &work[work_offset],
482 ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
487 dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
488 &v[v_offset], ldv, &work[work_offset], ldwork);
493 for (j = 1; j <= i__1; ++j) {
495 for (i__ = 1; i__ <= i__2; ++i__) {
496 c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
502 } else if (lsame_(side, "R")) {
504 /* Form C * H or C * H' where C = ( C1 C2 ) */
506 /* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
511 for (j = 1; j <= i__1; ++j) {
512 dcopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
513 work_dim1 + 1], &c__1);
519 dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
520 v[v_offset], ldv, &work[work_offset], ldwork);
523 /* W := W + C2 * V2' */
526 dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
527 c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) *
528 v_dim1 + 1], ldv, &c_b14, &work[work_offset],
532 /* W := W * T or W * T' */
534 dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
535 t_offset], ldt, &work[work_offset], ldwork);
541 /* C2 := C2 - W * V2 */
544 dgemm_("No transpose", "No transpose", m, &i__1, k, &
545 c_b25, &work[work_offset], ldwork, &v[(*k + 1) *
546 v_dim1 + 1], ldv, &c_b14, &c__[(*k + 1) * c_dim1
552 dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
553 &v[v_offset], ldv, &work[work_offset], ldwork);
558 for (j = 1; j <= i__1; ++j) {
560 for (i__ = 1; i__ <= i__2; ++i__) {
561 c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
571 /* Let V = ( V1 V2 ) (V2: last K columns) */
572 /* where V2 is unit lower triangular. */
574 if (lsame_(side, "L")) {
576 /* Form H * C or H' * C where C = ( C1 ) */
579 /* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
584 for (j = 1; j <= i__1; ++j) {
585 dcopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
586 work_dim1 + 1], &c__1);
592 dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
593 v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
597 /* W := W + C1'*V1' */
600 dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
601 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
602 work[work_offset], ldwork);
605 /* W := W * T' or W * T */
607 dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
608 t_offset], ldt, &work[work_offset], ldwork);
610 /* C := C - V' * W' */
614 /* C1 := C1 - V1' * W' */
617 dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
618 v_offset], ldv, &work[work_offset], ldwork, &
619 c_b14, &c__[c_offset], ldc);
624 dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
625 &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
626 work_offset], ldwork);
631 for (j = 1; j <= i__1; ++j) {
633 for (i__ = 1; i__ <= i__2; ++i__) {
634 c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
641 } else if (lsame_(side, "R")) {
643 /* Form C * H or C * H' where C = ( C1 C2 ) */
645 /* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
650 for (j = 1; j <= i__1; ++j) {
651 dcopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
652 j * work_dim1 + 1], &c__1);
658 dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
659 v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
663 /* W := W + C1 * V1' */
666 dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
667 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
668 work[work_offset], ldwork);
671 /* W := W * T or W * T' */
673 dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
674 t_offset], ldt, &work[work_offset], ldwork);
680 /* C1 := C1 - W * V1 */
683 dgemm_("No transpose", "No transpose", m, &i__1, k, &
684 c_b25, &work[work_offset], ldwork, &v[v_offset],
685 ldv, &c_b14, &c__[c_offset], ldc);
690 dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
691 &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
692 work_offset], ldwork);
697 for (j = 1; j <= i__1; ++j) {
699 for (i__ = 1; i__ <= i__2; ++i__) {
700 c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *