ergo
template_lapack_spgst.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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_SPGST_HEADER
38#define TEMPLATE_LAPACK_SPGST_HEADER
39
41
42template<class Treal>
43int template_lapack_spgst(const integer *itype, const char *uplo, const integer *n,
44 Treal *ap, const Treal *bp, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 March 31, 1993
50
51
52 Purpose
53 =======
54
55 DSPGST reduces a real symmetric-definite generalized eigenproblem
56 to standard form, using packed storage.
57
58 If ITYPE = 1, the problem is A*x = lambda*B*x,
59 and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
60
61 If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
62 B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
63
64 B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
65
66 Arguments
67 =========
68
69 ITYPE (input) INTEGER
70 = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
71 = 2 or 3: compute U*A*U**T or L**T*A*L.
72
73 UPLO (input) CHARACTER
74 = 'U': Upper triangle of A is stored and B is factored as
75 U**T*U;
76 = 'L': Lower triangle of A is stored and B is factored as
77 L*L**T.
78
79 N (input) INTEGER
80 The order of the matrices A and B. N >= 0.
81
82 AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
83 On entry, the upper or lower triangle of the symmetric matrix
84 A, packed columnwise in a linear array. The j-th column of A
85 is stored in the array AP as follows:
86 if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
87 if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
88
89 On exit, if INFO = 0, the transformed matrix, stored in the
90 same format as A.
91
92 BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
93 The triangular factor from the Cholesky factorization of B,
94 stored in the same format as A, as returned by DPPTRF.
95
96 INFO (output) INTEGER
97 = 0: successful exit
98 < 0: if INFO = -i, the i-th argument had an illegal value
99
100 =====================================================================
101
102
103 Test the input parameters.
104
105 Parameter adjustments */
106 /* Table of constant values */
107 integer c__1 = 1;
108 Treal c_b9 = -1.;
109 Treal c_b11 = 1.;
110
111 /* System generated locals */
112 integer i__1, i__2;
113 Treal d__1;
114 /* Local variables */
115 integer j, k;
116 logical upper;
117 integer j1, k1;
118 integer jj, kk;
119 Treal ct;
120 Treal ajj;
121 integer j1j1;
122 Treal akk;
123 integer k1k1;
124 Treal bjj, bkk;
125
126
127 --bp;
128 --ap;
129
130 /* Function Body */
131 *info = 0;
132 upper = template_blas_lsame(uplo, "U");
133 if (*itype < 1 || *itype > 3) {
134 *info = -1;
135 } else if (! upper && ! template_blas_lsame(uplo, "L")) {
136 *info = -2;
137 } else if (*n < 0) {
138 *info = -3;
139 }
140 if (*info != 0) {
141 i__1 = -(*info);
142 template_blas_erbla("SPGST ", &i__1);
143 return 0;
144 }
145
146 if (*itype == 1) {
147 if (upper) {
148
149/* Compute inv(U')*A*inv(U)
150
151 J1 and JJ are the indices of A(1,j) and A(j,j) */
152
153 jj = 0;
154 i__1 = *n;
155 for (j = 1; j <= i__1; ++j) {
156 j1 = jj + 1;
157 jj += j;
158
159/* Compute the j-th column of the upper triangle of A */
160
161 bjj = bp[jj];
162 template_blas_tpsv(uplo, "Transpose", "Nonunit", &j, &bp[1], &ap[j1], &
163 c__1);
164 i__2 = j - 1;
165 template_blas_spmv(uplo, &i__2, &c_b9, &ap[1], &bp[j1], &c__1, &c_b11, &
166 ap[j1], &c__1);
167 i__2 = j - 1;
168 d__1 = 1. / bjj;
169 template_blas_scal(&i__2, &d__1, &ap[j1], &c__1);
170 i__2 = j - 1;
171 ap[jj] = (ap[jj] - template_blas_dot(&i__2, &ap[j1], &c__1, &bp[j1], &
172 c__1)) / bjj;
173/* L10: */
174 }
175 } else {
176
177/* Compute inv(L)*A*inv(L')
178
179 KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) */
180
181 kk = 1;
182 i__1 = *n;
183 for (k = 1; k <= i__1; ++k) {
184 k1k1 = kk + *n - k + 1;
185
186/* Update the lower triangle of A(k:n,k:n) */
187
188 akk = ap[kk];
189 bkk = bp[kk];
190/* Computing 2nd power */
191 d__1 = bkk;
192 akk /= d__1 * d__1;
193 ap[kk] = akk;
194 if (k < *n) {
195 i__2 = *n - k;
196 d__1 = 1. / bkk;
197 template_blas_scal(&i__2, &d__1, &ap[kk + 1], &c__1);
198 ct = akk * -.5;
199 i__2 = *n - k;
200 template_blas_axpy(&i__2, &ct, &bp[kk + 1], &c__1, &ap[kk + 1], &c__1)
201 ;
202 i__2 = *n - k;
203 template_blas_spr2(uplo, &i__2, &c_b9, &ap[kk + 1], &c__1, &bp[kk + 1]
204 , &c__1, &ap[k1k1]);
205 i__2 = *n - k;
206 template_blas_axpy(&i__2, &ct, &bp[kk + 1], &c__1, &ap[kk + 1], &c__1)
207 ;
208 i__2 = *n - k;
209 template_blas_tpsv(uplo, "No transpose", "Non-unit", &i__2, &bp[k1k1],
210 &ap[kk + 1], &c__1);
211 }
212 kk = k1k1;
213/* L20: */
214 }
215 }
216 } else {
217 if (upper) {
218
219/* Compute U*A*U'
220
221 K1 and KK are the indices of A(1,k) and A(k,k) */
222
223 kk = 0;
224 i__1 = *n;
225 for (k = 1; k <= i__1; ++k) {
226 k1 = kk + 1;
227 kk += k;
228
229/* Update the upper triangle of A(1:k,1:k) */
230
231 akk = ap[kk];
232 bkk = bp[kk];
233 i__2 = k - 1;
234 template_blas_tpmv(uplo, "No transpose", "Non-unit", &i__2, &bp[1], &ap[
235 k1], &c__1);
236 ct = akk * .5;
237 i__2 = k - 1;
238 template_blas_axpy(&i__2, &ct, &bp[k1], &c__1, &ap[k1], &c__1);
239 i__2 = k - 1;
240 template_blas_spr2(uplo, &i__2, &c_b11, &ap[k1], &c__1, &bp[k1], &c__1, &
241 ap[1]);
242 i__2 = k - 1;
243 template_blas_axpy(&i__2, &ct, &bp[k1], &c__1, &ap[k1], &c__1);
244 i__2 = k - 1;
245 template_blas_scal(&i__2, &bkk, &ap[k1], &c__1);
246/* Computing 2nd power */
247 d__1 = bkk;
248 ap[kk] = akk * (d__1 * d__1);
249/* L30: */
250 }
251 } else {
252
253/* Compute L'*A*L
254
255 JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) */
256
257 jj = 1;
258 i__1 = *n;
259 for (j = 1; j <= i__1; ++j) {
260 j1j1 = jj + *n - j + 1;
261
262/* Compute the j-th column of the lower triangle of A */
263
264 ajj = ap[jj];
265 bjj = bp[jj];
266 i__2 = *n - j;
267 ap[jj] = ajj * bjj + template_blas_dot(&i__2, &ap[jj + 1], &c__1, &bp[jj
268 + 1], &c__1);
269 i__2 = *n - j;
270 template_blas_scal(&i__2, &bjj, &ap[jj + 1], &c__1);
271 i__2 = *n - j;
272 template_blas_spmv(uplo, &i__2, &c_b11, &ap[j1j1], &bp[jj + 1], &c__1, &
273 c_b11, &ap[jj + 1], &c__1);
274 i__2 = *n - j + 1;
275 template_blas_tpmv(uplo, "Transpose", "Non-unit", &i__2, &bp[jj], &ap[jj],
276 &c__1);
277 jj = j1j1;
278/* L40: */
279 }
280 }
281 }
282 return 0;
283
284/* End of DSPGST */
285
286} /* dspgst_ */
287
288#endif
int template_blas_axpy(const integer *n, const Treal *da, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition template_blas_axpy.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
bool logical
Definition template_blas_common.h:41
Treal template_blas_dot(const integer *n, const Treal *dx, const integer *incx, const Treal *dy, const integer *incy)
Definition template_blas_dot.h:43
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
int template_blas_spmv(const char *uplo, const integer *n, const Treal *alpha, Treal *ap, const Treal *x, const integer *incx, const Treal *beta, Treal *y, const integer *incy)
Definition template_blas_spmv.h:42
int template_blas_spr2(const char *uplo, const integer *n, const Treal *alpha, const Treal *x, const integer *incx, const Treal *y, const integer *incy, Treal *ap)
Definition template_blas_spr2.h:43
int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition template_blas_tpmv.h:42
int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition template_blas_tpsv.h:42
int template_lapack_spgst(const integer *itype, const char *uplo, const integer *n, Treal *ap, const Treal *bp, integer *info)
Definition template_lapack_spgst.h:43