ergo
template_lapack_spgst.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_SPGST_HEADER
38 #define TEMPLATE_LAPACK_SPGST_HEADER
39 
40 #include "template_lapack_common.h"
41 
42 template<class Treal>
43 int 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
template_blas_tpmv
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
template_blas_dot
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
template_blas_spmv
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
template_blas_scal
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
template_lapack_spgst
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
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_common.h
template_blas_spr2
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
template_blas_tpsv
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
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
template_blas_axpy
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
integer
int integer
Definition: template_blas_common.h:40