ergo
template_lapack_larft.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_LARFT_HEADER
38 #define TEMPLATE_LAPACK_LARFT_HEADER
39 
40 
41 template<class Treal>
42 int template_lapack_larft(const char *direct, const char *storev, const integer *n, const integer *
43  k, Treal *v, const integer *ldv, const Treal *tau, Treal *t,
44  const integer *ldt)
45 {
46 /* -- LAPACK auxiliary routine (version 3.0) --
47  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48  Courant Institute, Argonne National Lab, and Rice University
49  February 29, 1992
50 
51 
52  Purpose
53  =======
54 
55  DLARFT forms the triangular factor T of a real block reflector H
56  of order n, which is defined as a product of k elementary reflectors.
57 
58  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
59 
60  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
61 
62  If STOREV = 'C', the vector which defines the elementary reflector
63  H(i) is stored in the i-th column of the array V, and
64 
65  H = I - V * T * V'
66 
67  If STOREV = 'R', the vector which defines the elementary reflector
68  H(i) is stored in the i-th row of the array V, and
69 
70  H = I - V' * T * V
71 
72  Arguments
73  =========
74 
75  DIRECT (input) CHARACTER*1
76  Specifies the order in which the elementary reflectors are
77  multiplied to form the block reflector:
78  = 'F': H = H(1) H(2) . . . H(k) (Forward)
79  = 'B': H = H(k) . . . H(2) H(1) (Backward)
80 
81  STOREV (input) CHARACTER*1
82  Specifies how the vectors which define the elementary
83  reflectors are stored (see also Further Details):
84  = 'C': columnwise
85  = 'R': rowwise
86 
87  N (input) INTEGER
88  The order of the block reflector H. N >= 0.
89 
90  K (input) INTEGER
91  The order of the triangular factor T (= the number of
92  elementary reflectors). K >= 1.
93 
94  V (input/output) DOUBLE PRECISION array, dimension
95  (LDV,K) if STOREV = 'C'
96  (LDV,N) if STOREV = 'R'
97  The matrix V. See further details.
98 
99  LDV (input) INTEGER
100  The leading dimension of the array V.
101  If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
102 
103  TAU (input) DOUBLE PRECISION array, dimension (K)
104  TAU(i) must contain the scalar factor of the elementary
105  reflector H(i).
106 
107  T (output) DOUBLE PRECISION array, dimension (LDT,K)
108  The k by k triangular factor T of the block reflector.
109  If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
110  lower triangular. The rest of the array is not used.
111 
112  LDT (input) INTEGER
113  The leading dimension of the array T. LDT >= K.
114 
115  Further Details
116  ===============
117 
118  The shape of the matrix V and the storage of the vectors which define
119  the H(i) is best illustrated by the following example with n = 5 and
120  k = 3. The elements equal to 1 are not stored; the corresponding
121  array elements are modified but restored on exit. The rest of the
122  array is not used.
123 
124  DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
125 
126  V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
127  ( v1 1 ) ( 1 v2 v2 v2 )
128  ( v1 v2 1 ) ( 1 v3 v3 )
129  ( v1 v2 v3 )
130  ( v1 v2 v3 )
131 
132  DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
133 
134  V = ( v1 v2 v3 ) V = ( v1 v1 1 )
135  ( v1 v2 v3 ) ( v2 v2 v2 1 )
136  ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
137  ( 1 v3 )
138  ( 1 )
139 
140  =====================================================================
141 
142 
143  Quick return if possible
144 
145  Parameter adjustments */
146  /* Table of constant values */
147  integer c__1 = 1;
148  Treal c_b8 = 0.;
149 
150  /* System generated locals */
151  integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3;
152  Treal d__1;
153  /* Local variables */
154  integer i__, j;
155  Treal vii;
156 #define t_ref(a_1,a_2) t[(a_2)*t_dim1 + a_1]
157 #define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1]
158 
159 
160  v_dim1 = *ldv;
161  v_offset = 1 + v_dim1 * 1;
162  v -= v_offset;
163  --tau;
164  t_dim1 = *ldt;
165  t_offset = 1 + t_dim1 * 1;
166  t -= t_offset;
167 
168  /* Function Body */
169  if (*n == 0) {
170  return 0;
171  }
172 
173  if (template_blas_lsame(direct, "F")) {
174  i__1 = *k;
175  for (i__ = 1; i__ <= i__1; ++i__) {
176  if (tau[i__] == 0.) {
177 
178 /* H(i) = I */
179 
180  i__2 = i__;
181  for (j = 1; j <= i__2; ++j) {
182  t_ref(j, i__) = 0.;
183 /* L10: */
184  }
185  } else {
186 
187 /* general case */
188 
189  vii = v_ref(i__, i__);
190  v_ref(i__, i__) = 1.;
191  if (template_blas_lsame(storev, "C")) {
192 
193 /* T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) */
194 
195  i__2 = *n - i__ + 1;
196  i__3 = i__ - 1;
197  d__1 = -tau[i__];
198  template_blas_gemv("Transpose", &i__2, &i__3, &d__1, &v_ref(i__, 1),
199  ldv, &v_ref(i__, i__), &c__1, &c_b8, &t_ref(1,
200  i__), &c__1);
201  } else {
202 
203 /* T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' */
204 
205  i__2 = i__ - 1;
206  i__3 = *n - i__ + 1;
207  d__1 = -tau[i__];
208  template_blas_gemv("No transpose", &i__2, &i__3, &d__1, &v_ref(1, i__)
209  , ldv, &v_ref(i__, i__), ldv, &c_b8, &t_ref(1,
210  i__), &c__1);
211  }
212  v_ref(i__, i__) = vii;
213 
214 /* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */
215 
216  i__2 = i__ - 1;
217  template_blas_trmv("Upper", "No transpose", "Non-unit", &i__2, &t[
218  t_offset], ldt, &t_ref(1, i__), &c__1);
219  t_ref(i__, i__) = tau[i__];
220  }
221 /* L20: */
222  }
223  } else {
224  for (i__ = *k; i__ >= 1; --i__) {
225  if (tau[i__] == 0.) {
226 
227 /* H(i) = I */
228 
229  i__1 = *k;
230  for (j = i__; j <= i__1; ++j) {
231  t_ref(j, i__) = 0.;
232 /* L30: */
233  }
234  } else {
235 
236 /* general case */
237 
238  if (i__ < *k) {
239  if (template_blas_lsame(storev, "C")) {
240  vii = v_ref(*n - *k + i__, i__);
241  v_ref(*n - *k + i__, i__) = 1.;
242 
243 /* T(i+1:k,i) :=
244  - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) */
245 
246  i__1 = *n - *k + i__;
247  i__2 = *k - i__;
248  d__1 = -tau[i__];
249  template_blas_gemv("Transpose", &i__1, &i__2, &d__1, &v_ref(1,
250  i__ + 1), ldv, &v_ref(1, i__), &c__1, &c_b8, &
251  t_ref(i__ + 1, i__), &c__1);
252  v_ref(*n - *k + i__, i__) = vii;
253  } else {
254  vii = v_ref(i__, *n - *k + i__);
255  v_ref(i__, *n - *k + i__) = 1.;
256 
257 /* T(i+1:k,i) :=
258  - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' */
259 
260  i__1 = *k - i__;
261  i__2 = *n - *k + i__;
262  d__1 = -tau[i__];
263  template_blas_gemv("No transpose", &i__1, &i__2, &d__1, &v_ref(
264  i__ + 1, 1), ldv, &v_ref(i__, 1), ldv, &c_b8,
265  &t_ref(i__ + 1, i__), &c__1);
266  v_ref(i__, *n - *k + i__) = vii;
267  }
268 
269 /* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */
270 
271  i__1 = *k - i__;
272  template_blas_trmv("Lower", "No transpose", "Non-unit", &i__1, &t_ref(
273  i__ + 1, i__ + 1), ldt, &t_ref(i__ + 1, i__), &
274  c__1);
275  }
276  t_ref(i__, i__) = tau[i__];
277  }
278 /* L40: */
279  }
280  }
281  return 0;
282 
283 /* End of DLARFT */
284 
285 } /* dlarft_ */
286 
287 #undef v_ref
288 #undef t_ref
289 
290 
291 #endif
int integer
Definition: template_blas_common.h:40
int template_blas_trmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition: template_blas_trmv.h:42
#define t_ref(a_1, a_2)
#define v_ref(a_1, a_2)
int template_blas_gemv(const char *trans, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal *beta, Treal *y, const integer *incy)
Definition: template_blas_gemv.h:43
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int template_lapack_larft(const char *direct, const char *storev, const integer *n, const integer *k, Treal *v, const integer *ldv, const Treal *tau, Treal *t, const integer *ldt)
Definition: template_lapack_larft.h:42