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