ergo
template_lapack_larfb.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_LARFB_HEADER
38#define TEMPLATE_LAPACK_LARFB_HEADER
39
40
41template<class Treal>
42int template_lapack_larfb(const char *side, const char *trans, const char *direct, const char *
43 storev, const integer *m, const integer *n, const integer *k, const Treal *v, const integer *
44 ldv, const Treal *t, const integer *ldt, Treal *c__, const integer *ldc,
45 Treal *work, const integer *ldwork)
46{
47/* -- LAPACK auxiliary routine (version 3.0) --
48 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
49 Courant Institute, Argonne National Lab, and Rice University
50 February 29, 1992
51
52
53 Purpose
54 =======
55
56 DLARFB applies a real block reflector H or its transpose H' to a
57 real m by n matrix C, from either the left or the right.
58
59 Arguments
60 =========
61
62 SIDE (input) CHARACTER*1
63 = 'L': apply H or H' from the Left
64 = 'R': apply H or H' from the Right
65
66 TRANS (input) CHARACTER*1
67 = 'N': apply H (No transpose)
68 = 'T': apply H' (Transpose)
69
70 DIRECT (input) CHARACTER*1
71 Indicates how H is formed from a product of elementary
72 reflectors
73 = 'F': H = H(1) H(2) . . . H(k) (Forward)
74 = 'B': H = H(k) . . . H(2) H(1) (Backward)
75
76 STOREV (input) CHARACTER*1
77 Indicates how the vectors which define the elementary
78 reflectors are stored:
79 = 'C': Columnwise
80 = 'R': Rowwise
81
82 M (input) INTEGER
83 The number of rows of the matrix C.
84
85 N (input) INTEGER
86 The number of columns of the matrix C.
87
88 K (input) INTEGER
89 The order of the matrix T (= the number of elementary
90 reflectors whose product defines the block reflector).
91
92 V (input) DOUBLE PRECISION array, dimension
93 (LDV,K) if STOREV = 'C'
94 (LDV,M) if STOREV = 'R' and SIDE = 'L'
95 (LDV,N) if STOREV = 'R' and SIDE = 'R'
96 The matrix V. See further details.
97
98 LDV (input) INTEGER
99 The leading dimension of the array V.
100 If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
101 if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
102 if STOREV = 'R', LDV >= K.
103
104 T (input) DOUBLE PRECISION array, dimension (LDT,K)
105 The triangular k by k matrix T in the representation of the
106 block reflector.
107
108 LDT (input) INTEGER
109 The leading dimension of the array T. LDT >= K.
110
111 C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
112 On entry, the m by n matrix C.
113 On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
114
115 LDC (input) INTEGER
116 The leading dimension of the array C. LDA >= max(1,M).
117
118 WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
119
120 LDWORK (input) INTEGER
121 The leading dimension of the array WORK.
122 If SIDE = 'L', LDWORK >= max(1,N);
123 if SIDE = 'R', LDWORK >= max(1,M).
124
125 =====================================================================
126
127
128 Quick return if possible
129
130 Parameter adjustments */
131 /* Table of constant values */
132 integer c__1 = 1;
133 Treal c_b14 = 1.;
134 Treal c_b25 = -1.;
135
136 /* System generated locals */
137 integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1,
138 work_offset, i__1, i__2;
139 /* Local variables */
140 integer i__, j;
141 char transt[1];
142#define work_ref(a_1,a_2) work[(a_2)*work_dim1 + a_1]
143#define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
144#define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1]
145
146
147 v_dim1 = *ldv;
148 v_offset = 1 + v_dim1 * 1;
149 v -= v_offset;
150 t_dim1 = *ldt;
151 t_offset = 1 + t_dim1 * 1;
152 t -= t_offset;
153 c_dim1 = *ldc;
154 c_offset = 1 + c_dim1 * 1;
155 c__ -= c_offset;
156 work_dim1 = *ldwork;
157 work_offset = 1 + work_dim1 * 1;
158 work -= work_offset;
159
160 /* Function Body */
161 if (*m <= 0 || *n <= 0) {
162 return 0;
163 }
164
165 if (template_blas_lsame(trans, "N")) {
166 *(unsigned char *)transt = 'T';
167 } else {
168 *(unsigned char *)transt = 'N';
169 }
170
171 if (template_blas_lsame(storev, "C")) {
172
173 if (template_blas_lsame(direct, "F")) {
174
175/* Let V = ( V1 ) (first K rows)
176 ( V2 )
177 where V1 is unit lower triangular. */
178
179 if (template_blas_lsame(side, "L")) {
180
181/* Form H * C or H' * C where C = ( C1 )
182 ( C2 )
183
184 W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
185
186 W := C1' */
187
188 i__1 = *k;
189 for (j = 1; j <= i__1; ++j) {
190 template_blas_copy(n, &c___ref(j, 1), ldc, &work_ref(1, j), &c__1);
191/* L10: */
192 }
193
194/* W := W * V1 */
195
196 template_blas_trmm("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
197 &v[v_offset], ldv, &work[work_offset], ldwork);
198 if (*m > *k) {
199
200/* W := W + C2'*V2 */
201
202 i__1 = *m - *k;
203 template_blas_gemm("Transpose", "No transpose", n, k, &i__1, &c_b14, &
204 c___ref(*k + 1, 1), ldc, &v_ref(*k + 1, 1), ldv, &
205 c_b14, &work[work_offset], ldwork);
206 }
207
208/* W := W * T' or W * T */
209
210 template_blas_trmm("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
211 t_offset], ldt, &work[work_offset], ldwork);
212
213/* C := C - V * W' */
214
215 if (*m > *k) {
216
217/* C2 := C2 - V2 * W' */
218
219 i__1 = *m - *k;
220 template_blas_gemm("No transpose", "Transpose", &i__1, n, k, &c_b25, &
221 v_ref(*k + 1, 1), ldv, &work[work_offset], ldwork,
222 &c_b14, &c___ref(*k + 1, 1), ldc);
223 }
224
225/* W := W * V1' */
226
227 template_blas_trmm("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
228 v[v_offset], ldv, &work[work_offset], ldwork);
229
230/* C1 := C1 - W' */
231
232 i__1 = *k;
233 for (j = 1; j <= i__1; ++j) {
234 i__2 = *n;
235 for (i__ = 1; i__ <= i__2; ++i__) {
236 c___ref(j, i__) = c___ref(j, i__) - work_ref(i__, j);
237/* L20: */
238 }
239/* L30: */
240 }
241
242 } else if (template_blas_lsame(side, "R")) {
243
244/* Form C * H or C * H' where C = ( C1 C2 )
245
246 W := C * V = (C1*V1 + C2*V2) (stored in WORK)
247
248 W := C1 */
249
250 i__1 = *k;
251 for (j = 1; j <= i__1; ++j) {
252 template_blas_copy(m, &c___ref(1, j), &c__1, &work_ref(1, j), &c__1);
253/* L40: */
254 }
255
256/* W := W * V1 */
257
258 template_blas_trmm("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
259 &v[v_offset], ldv, &work[work_offset], ldwork);
260 if (*n > *k) {
261
262/* W := W + C2 * V2 */
263
264 i__1 = *n - *k;
265 template_blas_gemm("No transpose", "No transpose", m, k, &i__1, &
266 c_b14, &c___ref(1, *k + 1), ldc, &v_ref(*k + 1, 1)
267 , ldv, &c_b14, &work[work_offset], ldwork);
268 }
269
270/* W := W * T or W * T' */
271
272 template_blas_trmm("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
273 t_offset], ldt, &work[work_offset], ldwork);
274
275/* C := C - W * V' */
276
277 if (*n > *k) {
278
279/* C2 := C2 - W * V2' */
280
281 i__1 = *n - *k;
282 template_blas_gemm("No transpose", "Transpose", m, &i__1, k, &c_b25, &
283 work[work_offset], ldwork, &v_ref(*k + 1, 1), ldv,
284 &c_b14, &c___ref(1, *k + 1), ldc);
285 }
286
287/* W := W * V1' */
288
289 template_blas_trmm("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
290 v[v_offset], ldv, &work[work_offset], ldwork);
291
292/* C1 := C1 - W */
293
294 i__1 = *k;
295 for (j = 1; j <= i__1; ++j) {
296 i__2 = *m;
297 for (i__ = 1; i__ <= i__2; ++i__) {
298 c___ref(i__, j) = c___ref(i__, j) - work_ref(i__, j);
299/* L50: */
300 }
301/* L60: */
302 }
303 }
304
305 } else {
306
307/* Let V = ( V1 )
308 ( V2 ) (last K rows)
309 where V2 is unit upper triangular. */
310
311 if (template_blas_lsame(side, "L")) {
312
313/* Form H * C or H' * C where C = ( C1 )
314 ( C2 )
315
316 W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
317
318 W := C2' */
319
320 i__1 = *k;
321 for (j = 1; j <= i__1; ++j) {
322 template_blas_copy(n, &c___ref(*m - *k + j, 1), ldc, &work_ref(1, j),
323 &c__1);
324/* L70: */
325 }
326
327/* W := W * V2 */
328
329 template_blas_trmm("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
330 &v_ref(*m - *k + 1, 1), ldv, &work[work_offset],
331 ldwork);
332 if (*m > *k) {
333
334/* W := W + C1'*V1 */
335
336 i__1 = *m - *k;
337 template_blas_gemm("Transpose", "No transpose", n, k, &i__1, &c_b14, &
338 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
339 work[work_offset], ldwork);
340 }
341
342/* W := W * T' or W * T */
343
344 template_blas_trmm("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
345 t_offset], ldt, &work[work_offset], ldwork);
346
347/* C := C - V * W' */
348
349 if (*m > *k) {
350
351/* C1 := C1 - V1 * W' */
352
353 i__1 = *m - *k;
354 template_blas_gemm("No transpose", "Transpose", &i__1, n, k, &c_b25, &
355 v[v_offset], ldv, &work[work_offset], ldwork, &
356 c_b14, &c__[c_offset], ldc)
357 ;
358 }
359
360/* W := W * V2' */
361
362 template_blas_trmm("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
363 v_ref(*m - *k + 1, 1), ldv, &work[work_offset],
364 ldwork);
365
366/* C2 := C2 - W' */
367
368 i__1 = *k;
369 for (j = 1; j <= i__1; ++j) {
370 i__2 = *n;
371 for (i__ = 1; i__ <= i__2; ++i__) {
372 c___ref(*m - *k + j, i__) = c___ref(*m - *k + j, i__)
373 - work_ref(i__, j);
374/* L80: */
375 }
376/* L90: */
377 }
378
379 } else if (template_blas_lsame(side, "R")) {
380
381/* Form C * H or C * H' where C = ( C1 C2 )
382
383 W := C * V = (C1*V1 + C2*V2) (stored in WORK)
384
385 W := C2 */
386
387 i__1 = *k;
388 for (j = 1; j <= i__1; ++j) {
389 template_blas_copy(m, &c___ref(1, *n - *k + j), &c__1, &work_ref(1, j)
390 , &c__1);
391/* L100: */
392 }
393
394/* W := W * V2 */
395
396 template_blas_trmm("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
397 &v_ref(*n - *k + 1, 1), ldv, &work[work_offset],
398 ldwork);
399 if (*n > *k) {
400
401/* W := W + C1 * V1 */
402
403 i__1 = *n - *k;
404 template_blas_gemm("No transpose", "No transpose", m, k, &i__1, &
405 c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, &
406 c_b14, &work[work_offset], ldwork);
407 }
408
409/* W := W * T or W * T' */
410
411 template_blas_trmm("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
412 t_offset], ldt, &work[work_offset], ldwork);
413
414/* C := C - W * V' */
415
416 if (*n > *k) {
417
418/* C1 := C1 - W * V1' */
419
420 i__1 = *n - *k;
421 template_blas_gemm("No transpose", "Transpose", m, &i__1, k, &c_b25, &
422 work[work_offset], ldwork, &v[v_offset], ldv, &
423 c_b14, &c__[c_offset], ldc)
424 ;
425 }
426
427/* W := W * V2' */
428
429 template_blas_trmm("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
430 v_ref(*n - *k + 1, 1), ldv, &work[work_offset],
431 ldwork);
432
433/* C2 := C2 - W */
434
435 i__1 = *k;
436 for (j = 1; j <= i__1; ++j) {
437 i__2 = *m;
438 for (i__ = 1; i__ <= i__2; ++i__) {
439 c___ref(i__, *n - *k + j) = c___ref(i__, *n - *k + j)
440 - work_ref(i__, j);
441/* L110: */
442 }
443/* L120: */
444 }
445 }
446 }
447
448 } else if (template_blas_lsame(storev, "R")) {
449
450 if (template_blas_lsame(direct, "F")) {
451
452/* Let V = ( V1 V2 ) (V1: first K columns)
453 where V1 is unit upper triangular. */
454
455 if (template_blas_lsame(side, "L")) {
456
457/* Form H * C or H' * C where C = ( C1 )
458 ( C2 )
459
460 W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
461
462 W := C1' */
463
464 i__1 = *k;
465 for (j = 1; j <= i__1; ++j) {
466 template_blas_copy(n, &c___ref(j, 1), ldc, &work_ref(1, j), &c__1);
467/* L130: */
468 }
469
470/* W := W * V1' */
471
472 template_blas_trmm("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
473 v[v_offset], ldv, &work[work_offset], ldwork);
474 if (*m > *k) {
475
476/* W := W + C2'*V2' */
477
478 i__1 = *m - *k;
479 template_blas_gemm("Transpose", "Transpose", n, k, &i__1, &c_b14, &
480 c___ref(*k + 1, 1), ldc, &v_ref(1, *k + 1), ldv, &
481 c_b14, &work[work_offset], ldwork);
482 }
483
484/* W := W * T' or W * T */
485
486 template_blas_trmm("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
487 t_offset], ldt, &work[work_offset], ldwork);
488
489/* C := C - V' * W' */
490
491 if (*m > *k) {
492
493/* C2 := C2 - V2' * W' */
494
495 i__1 = *m - *k;
496 template_blas_gemm("Transpose", "Transpose", &i__1, n, k, &c_b25, &
497 v_ref(1, *k + 1), ldv, &work[work_offset], ldwork,
498 &c_b14, &c___ref(*k + 1, 1), ldc);
499 }
500
501/* W := W * V1 */
502
503 template_blas_trmm("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
504 &v[v_offset], ldv, &work[work_offset], ldwork);
505
506/* C1 := C1 - W' */
507
508 i__1 = *k;
509 for (j = 1; j <= i__1; ++j) {
510 i__2 = *n;
511 for (i__ = 1; i__ <= i__2; ++i__) {
512 c___ref(j, i__) = c___ref(j, i__) - work_ref(i__, j);
513/* L140: */
514 }
515/* L150: */
516 }
517
518 } else if (template_blas_lsame(side, "R")) {
519
520/* Form C * H or C * H' where C = ( C1 C2 )
521
522 W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
523
524 W := C1 */
525
526 i__1 = *k;
527 for (j = 1; j <= i__1; ++j) {
528 template_blas_copy(m, &c___ref(1, j), &c__1, &work_ref(1, j), &c__1);
529/* L160: */
530 }
531
532/* W := W * V1' */
533
534 template_blas_trmm("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
535 v[v_offset], ldv, &work[work_offset], ldwork);
536 if (*n > *k) {
537
538/* W := W + C2 * V2' */
539
540 i__1 = *n - *k;
541 template_blas_gemm("No transpose", "Transpose", m, k, &i__1, &c_b14, &
542 c___ref(1, *k + 1), ldc, &v_ref(1, *k + 1), ldv, &
543 c_b14, &work[work_offset], ldwork);
544 }
545
546/* W := W * T or W * T' */
547
548 template_blas_trmm("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
549 t_offset], ldt, &work[work_offset], ldwork);
550
551/* C := C - W * V */
552
553 if (*n > *k) {
554
555/* C2 := C2 - W * V2 */
556
557 i__1 = *n - *k;
558 template_blas_gemm("No transpose", "No transpose", m, &i__1, k, &
559 c_b25, &work[work_offset], ldwork, &v_ref(1, *k +
560 1), ldv, &c_b14, &c___ref(1, *k + 1), ldc);
561 }
562
563/* W := W * V1 */
564
565 template_blas_trmm("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
566 &v[v_offset], ldv, &work[work_offset], ldwork);
567
568/* C1 := C1 - W */
569
570 i__1 = *k;
571 for (j = 1; j <= i__1; ++j) {
572 i__2 = *m;
573 for (i__ = 1; i__ <= i__2; ++i__) {
574 c___ref(i__, j) = c___ref(i__, j) - work_ref(i__, j);
575/* L170: */
576 }
577/* L180: */
578 }
579
580 }
581
582 } else {
583
584/* Let V = ( V1 V2 ) (V2: last K columns)
585 where V2 is unit lower triangular. */
586
587 if (template_blas_lsame(side, "L")) {
588
589/* Form H * C or H' * C where C = ( C1 )
590 ( C2 )
591
592 W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
593
594 W := C2' */
595
596 i__1 = *k;
597 for (j = 1; j <= i__1; ++j) {
598 template_blas_copy(n, &c___ref(*m - *k + j, 1), ldc, &work_ref(1, j),
599 &c__1);
600/* L190: */
601 }
602
603/* W := W * V2' */
604
605 template_blas_trmm("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
606 v_ref(1, *m - *k + 1), ldv, &work[work_offset],
607 ldwork);
608 if (*m > *k) {
609
610/* W := W + C1'*V1' */
611
612 i__1 = *m - *k;
613 template_blas_gemm("Transpose", "Transpose", n, k, &i__1, &c_b14, &
614 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
615 work[work_offset], ldwork);
616 }
617
618/* W := W * T' or W * T */
619
620 template_blas_trmm("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
621 t_offset], ldt, &work[work_offset], ldwork);
622
623/* C := C - V' * W' */
624
625 if (*m > *k) {
626
627/* C1 := C1 - V1' * W' */
628
629 i__1 = *m - *k;
630 template_blas_gemm("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
631 v_offset], ldv, &work[work_offset], ldwork, &
632 c_b14, &c__[c_offset], ldc);
633 }
634
635/* W := W * V2 */
636
637 template_blas_trmm("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
638 &v_ref(1, *m - *k + 1), ldv, &work[work_offset],
639 ldwork);
640
641/* C2 := C2 - W' */
642
643 i__1 = *k;
644 for (j = 1; j <= i__1; ++j) {
645 i__2 = *n;
646 for (i__ = 1; i__ <= i__2; ++i__) {
647 c___ref(*m - *k + j, i__) = c___ref(*m - *k + j, i__)
648 - work_ref(i__, j);
649/* L200: */
650 }
651/* L210: */
652 }
653
654 } else if (template_blas_lsame(side, "R")) {
655
656/* Form C * H or C * H' where C = ( C1 C2 )
657
658 W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
659
660 W := C2 */
661
662 i__1 = *k;
663 for (j = 1; j <= i__1; ++j) {
664 template_blas_copy(m, &c___ref(1, *n - *k + j), &c__1, &work_ref(1, j)
665 , &c__1);
666/* L220: */
667 }
668
669/* W := W * V2' */
670
671 template_blas_trmm("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
672 v_ref(1, *n - *k + 1), ldv, &work[work_offset],
673 ldwork);
674 if (*n > *k) {
675
676/* W := W + C1 * V1' */
677
678 i__1 = *n - *k;
679 template_blas_gemm("No transpose", "Transpose", m, k, &i__1, &c_b14, &
680 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
681 work[work_offset], ldwork);
682 }
683
684/* W := W * T or W * T' */
685
686 template_blas_trmm("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
687 t_offset], ldt, &work[work_offset], ldwork);
688
689/* C := C - W * V */
690
691 if (*n > *k) {
692
693/* C1 := C1 - W * V1 */
694
695 i__1 = *n - *k;
696 template_blas_gemm("No transpose", "No transpose", m, &i__1, k, &
697 c_b25, &work[work_offset], ldwork, &v[v_offset],
698 ldv, &c_b14, &c__[c_offset], ldc);
699 }
700
701/* W := W * V2 */
702
703 template_blas_trmm("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
704 &v_ref(1, *n - *k + 1), ldv, &work[work_offset],
705 ldwork);
706
707/* C1 := C1 - W */
708
709 i__1 = *k;
710 for (j = 1; j <= i__1; ++j) {
711 i__2 = *m;
712 for (i__ = 1; i__ <= i__2; ++i__) {
713 c___ref(i__, *n - *k + j) = c___ref(i__, *n - *k + j)
714 - work_ref(i__, j);
715/* L230: */
716 }
717/* L240: */
718 }
719
720 }
721
722 }
723 }
724
725 return 0;
726
727/* End of DLARFB */
728
729} /* dlarfb_ */
730
731#undef v_ref
732#undef c___ref
733#undef work_ref
734
735
736#endif
side
Definition: Matrix.h:75
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:42
int template_blas_gemm(const char *transa, const char *transb, const integer *m, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition: template_blas_gemm.h:43
int template_blas_trmm(const char *side, const char *uplo, const char *transa, const char *diag, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, Treal *b, const integer *ldb)
Definition: template_blas_trmm.h:43
#define c___ref(a_1, a_2)
#define work_ref(a_1, a_2)
#define v_ref(a_1, a_2)
int template_lapack_larfb(const char *side, const char *trans, const char *direct, const char *storev, const integer *m, const integer *n, const integer *k, const Treal *v, const integer *ldv, const Treal *t, const integer *ldt, Treal *c__, const integer *ldc, Treal *work, const integer *ldwork)
Definition: template_lapack_larfb.h:42