ergo
template_lapack_larfb.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
41 template<class Treal>
42 int 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
int integer
Definition: template_blas_common.h:40
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 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
side
Definition: Matrix.h:75
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:42
#define c___ref(a_1, a_2)
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
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46