ergo
template_lapack_sygst.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_SYGST_HEADER
38 #define TEMPLATE_LAPACK_SYGST_HEADER
39 
40 #include "template_lapack_sygs2.h"
41 
42 template<class Treal>
43 int template_lapack_sygst(const integer *itype, const char *uplo, const integer *n,
44  Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *
45  info)
46 {
47 /* -- LAPACK routine (version 3.0) --
48  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
49  Courant Institute, Argonne National Lab, and Rice University
50  September 30, 1994
51 
52 
53  Purpose
54  =======
55 
56  DSYGST reduces a real symmetric-definite generalized eigenproblem
57  to standard form.
58 
59  If ITYPE = 1, the problem is A*x = lambda*B*x,
60  and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
61 
62  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
63  B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
64 
65  B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
66 
67  Arguments
68  =========
69 
70  ITYPE (input) INTEGER
71  = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
72  = 2 or 3: compute U*A*U**T or L**T*A*L.
73 
74  UPLO (input) CHARACTER
75  = 'U': Upper triangle of A is stored and B is factored as
76  U**T*U;
77  = 'L': Lower triangle of A is stored and B is factored as
78  L*L**T.
79 
80  N (input) INTEGER
81  The order of the matrices A and B. N >= 0.
82 
83  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
84  On entry, the symmetric matrix A. If UPLO = 'U', the leading
85  N-by-N upper triangular part of A contains the upper
86  triangular part of the matrix A, and the strictly lower
87  triangular part of A is not referenced. If UPLO = 'L', the
88  leading N-by-N lower triangular part of A contains the lower
89  triangular part of the matrix A, and the strictly upper
90  triangular part of A is not referenced.
91 
92  On exit, if INFO = 0, the transformed matrix, stored in the
93  same format as A.
94 
95  LDA (input) INTEGER
96  The leading dimension of the array A. LDA >= max(1,N).
97 
98  B (input) DOUBLE PRECISION array, dimension (LDB,N)
99  The triangular factor from the Cholesky factorization of B,
100  as returned by DPOTRF.
101 
102  LDB (input) INTEGER
103  The leading dimension of the array B. LDB >= max(1,N).
104 
105  INFO (output) INTEGER
106  = 0: successful exit
107  < 0: if INFO = -i, the i-th argument had an illegal value
108 
109  =====================================================================
110 
111 
112  Test the input parameters.
113 
114  Parameter adjustments */
115  /* Table of constant values */
116  integer c__1 = 1;
117  integer c_n1 = -1;
118  Treal c_b14 = 1.;
119  Treal c_b16 = -.5;
120  Treal c_b19 = -1.;
121  Treal c_b52 = .5;
122 
123  /* System generated locals */
124  integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
125  /* Local variables */
126  integer k;
127  logical upper;
128  integer kb;
129  integer nb;
130 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
131 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
132 
133 
134  a_dim1 = *lda;
135  a_offset = 1 + a_dim1 * 1;
136  a -= a_offset;
137  b_dim1 = *ldb;
138  b_offset = 1 + b_dim1 * 1;
139  b -= b_offset;
140 
141  /* Function Body */
142  *info = 0;
143  upper = template_blas_lsame(uplo, "U");
144  if (*itype < 1 || *itype > 3) {
145  *info = -1;
146  } else if (! upper && ! template_blas_lsame(uplo, "L")) {
147  *info = -2;
148  } else if (*n < 0) {
149  *info = -3;
150  } else if (*lda < maxMACRO(1,*n)) {
151  *info = -5;
152  } else if (*ldb < maxMACRO(1,*n)) {
153  *info = -7;
154  }
155  if (*info != 0) {
156  i__1 = -(*info);
157  template_blas_erbla("SYGST ", &i__1);
158  return 0;
159  }
160 
161 /* Quick return if possible */
162 
163  if (*n == 0) {
164  return 0;
165  }
166 
167 /* Determine the block size for this environment. */
168 
169  nb = template_lapack_ilaenv(&c__1, "DSYGST", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
170  ftnlen)1);
171 
172  if (nb <= 1 || nb >= *n) {
173 
174 /* Use unblocked code */
175 
176  template_lapack_sygs2(itype, uplo, n, &a[a_offset], lda, &b[b_offset], ldb, info);
177  } else {
178 
179 /* Use blocked code */
180 
181  if (*itype == 1) {
182  if (upper) {
183 
184 /* Compute inv(U')*A*inv(U) */
185 
186  i__1 = *n;
187  i__2 = nb;
188  for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
189 /* Computing MIN */
190  i__3 = *n - k + 1;
191  kb = minMACRO(i__3,nb);
192 
193 /* Update the upper triangle of A(k:n,k:n) */
194 
195  template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
196  ldb, info);
197  if (k + kb <= *n) {
198  i__3 = *n - k - kb + 1;
199  template_blas_trsm("Left", uplo, "Transpose", "Non-unit", &kb, &
200  i__3, &c_b14, &b_ref(k, k), ldb, &a_ref(k, k
201  + kb), lda);
202  i__3 = *n - k - kb + 1;
203  template_blas_symm("Left", uplo, &kb, &i__3, &c_b16, &a_ref(k, k),
204  lda, &b_ref(k, k + kb), ldb, &c_b14, &a_ref(
205  k, k + kb), lda);
206  i__3 = *n - k - kb + 1;
207  template_blas_syr2k(uplo, "Transpose", &i__3, &kb, &c_b19, &a_ref(
208  k, k + kb), lda, &b_ref(k, k + kb), ldb, &
209  c_b14, &a_ref(k + kb, k + kb), lda);
210  i__3 = *n - k - kb + 1;
211  template_blas_symm("Left", uplo, &kb, &i__3, &c_b16, &a_ref(k, k),
212  lda, &b_ref(k, k + kb), ldb, &c_b14, &a_ref(
213  k, k + kb), lda);
214  i__3 = *n - k - kb + 1;
215  template_blas_trsm("Right", uplo, "No transpose", "Non-unit", &kb,
216  &i__3, &c_b14, &b_ref(k + kb, k + kb), ldb, &
217  a_ref(k, k + kb), lda);
218  }
219 /* L10: */
220  }
221  } else {
222 
223 /* Compute inv(L)*A*inv(L') */
224 
225  i__2 = *n;
226  i__1 = nb;
227  for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
228 /* Computing MIN */
229  i__3 = *n - k + 1;
230  kb = minMACRO(i__3,nb);
231 
232 /* Update the lower triangle of A(k:n,k:n) */
233 
234  template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
235  ldb, info);
236  if (k + kb <= *n) {
237  i__3 = *n - k - kb + 1;
238  template_blas_trsm("Right", uplo, "Transpose", "Non-unit", &i__3,
239  &kb, &c_b14, &b_ref(k, k), ldb, &a_ref(k + kb,
240  k), lda);
241  i__3 = *n - k - kb + 1;
242  template_blas_symm("Right", uplo, &i__3, &kb, &c_b16, &a_ref(k, k)
243  , lda, &b_ref(k + kb, k), ldb, &c_b14, &a_ref(
244  k + kb, k), lda);
245  i__3 = *n - k - kb + 1;
246  template_blas_syr2k(uplo, "No transpose", &i__3, &kb, &c_b19, &
247  a_ref(k + kb, k), lda, &b_ref(k + kb, k), ldb,
248  &c_b14, &a_ref(k + kb, k + kb), lda);
249  i__3 = *n - k - kb + 1;
250  template_blas_symm("Right", uplo, &i__3, &kb, &c_b16, &a_ref(k, k)
251  , lda, &b_ref(k + kb, k), ldb, &c_b14, &a_ref(
252  k + kb, k), lda);
253  i__3 = *n - k - kb + 1;
254  template_blas_trsm("Left", uplo, "No transpose", "Non-unit", &
255  i__3, &kb, &c_b14, &b_ref(k + kb, k + kb),
256  ldb, &a_ref(k + kb, k), lda);
257  }
258 /* L20: */
259  }
260  }
261  } else {
262  if (upper) {
263 
264 /* Compute U*A*U' */
265 
266  i__1 = *n;
267  i__2 = nb;
268  for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
269 /* Computing MIN */
270  i__3 = *n - k + 1;
271  kb = minMACRO(i__3,nb);
272 
273 /* Update the upper triangle of A(1:k+kb-1,1:k+kb-1) */
274 
275  i__3 = k - 1;
276  template_blas_trmm("Left", uplo, "No transpose", "Non-unit", &i__3, &
277  kb, &c_b14, &b[b_offset], ldb, &a_ref(1, k), lda);
278  i__3 = k - 1;
279  template_blas_symm("Right", uplo, &i__3, &kb, &c_b52, &a_ref(k, k),
280  lda, &b_ref(1, k), ldb, &c_b14, &a_ref(1, k), lda);
281  i__3 = k - 1;
282  template_blas_syr2k(uplo, "No transpose", &i__3, &kb, &c_b14, &a_ref(
283  1, k), lda, &b_ref(1, k), ldb, &c_b14, &a[
284  a_offset], lda);
285  i__3 = k - 1;
286  template_blas_symm("Right", uplo, &i__3, &kb, &c_b52, &a_ref(k, k),
287  lda, &b_ref(1, k), ldb, &c_b14, &a_ref(1, k), lda);
288  i__3 = k - 1;
289  template_blas_trmm("Right", uplo, "Transpose", "Non-unit", &i__3, &kb,
290  &c_b14, &b_ref(k, k), ldb, &a_ref(1, k), lda);
291  template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
292  ldb, info);
293 /* L30: */
294  }
295  } else {
296 
297 /* Compute L'*A*L */
298 
299  i__2 = *n;
300  i__1 = nb;
301  for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
302 /* Computing MIN */
303  i__3 = *n - k + 1;
304  kb = minMACRO(i__3,nb);
305 
306 /* Update the lower triangle of A(1:k+kb-1,1:k+kb-1) */
307 
308  i__3 = k - 1;
309  template_blas_trmm("Right", uplo, "No transpose", "Non-unit", &kb, &
310  i__3, &c_b14, &b[b_offset], ldb, &a_ref(k, 1),
311  lda);
312  i__3 = k - 1;
313  template_blas_symm("Left", uplo, &kb, &i__3, &c_b52, &a_ref(k, k),
314  lda, &b_ref(k, 1), ldb, &c_b14, &a_ref(k, 1), lda);
315  i__3 = k - 1;
316  template_blas_syr2k(uplo, "Transpose", &i__3, &kb, &c_b14, &a_ref(k,
317  1), lda, &b_ref(k, 1), ldb, &c_b14, &a[a_offset],
318  lda);
319  i__3 = k - 1;
320  template_blas_symm("Left", uplo, &kb, &i__3, &c_b52, &a_ref(k, k),
321  lda, &b_ref(k, 1), ldb, &c_b14, &a_ref(k, 1), lda);
322  i__3 = k - 1;
323  template_blas_trmm("Left", uplo, "Transpose", "Non-unit", &kb, &i__3,
324  &c_b14, &b_ref(k, k), ldb, &a_ref(k, 1), lda);
325  template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
326  ldb, info);
327 /* L40: */
328  }
329  }
330  }
331  }
332  return 0;
333 
334 /* End of DSYGST */
335 
336 } /* dsygst_ */
337 
338 #undef b_ref
339 #undef a_ref
340 
341 
342 #endif
ftnlen
int ftnlen
Definition: template_blas_common.h:42
template_lapack_sygst
int template_lapack_sygst(const integer *itype, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *info)
Definition: template_lapack_sygst.h:43
template_lapack_ilaenv
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition: template_lapack_common.cc:281
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
template_blas_trsm
int template_blas_trsm(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_trsm.h:43
template_lapack_sygs2.h
template_lapack_sygs2
int template_lapack_sygs2(const integer *itype, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *info)
Definition: template_lapack_sygs2.h:43
logical
bool logical
Definition: template_blas_common.h:41
template_blas_syr2k
int template_blas_syr2k(const char *uplo, const char *trans, 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_syr2k.h:42
template_blas_trmm
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
a_ref
#define a_ref(a_1, a_2)
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_symm
int template_blas_symm(const char *side, const char *uplo, const integer *m, const integer *n, 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_symm.h:42
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
b_ref
#define b_ref(a_1, a_2)
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45