ergo
template_lapack_gghrd.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_GGHRD_HEADER
38#define TEMPLATE_LAPACK_GGHRD_HEADER
39
40
41template<class Treal>
42int template_lapack_gghrd(const char *compq, const char *compz, const integer *n, const integer *
43 ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b,
44 const integer *ldb, Treal *q, const integer *ldq, Treal *z__, const integer *
45 ldz, integer *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 DGGHRD reduces a pair of real matrices (A,B) to generalized upper
57 Hessenberg form using orthogonal transformations, where A is a
58 general matrix and B is upper triangular: Q' * A * Z = H and
59 Q' * B * Z = T, where H is upper Hessenberg, T is upper triangular,
60 and Q and Z are orthogonal, and ' means transpose.
61
62 The orthogonal matrices Q and Z are determined as products of Givens
63 rotations. They may either be formed explicitly, or they may be
64 postmultiplied into input matrices Q1 and Z1, so that
65
66 Q1 * A * Z1' = (Q1*Q) * H * (Z1*Z)'
67 Q1 * B * Z1' = (Q1*Q) * T * (Z1*Z)'
68
69 Arguments
70 =========
71
72 COMPQ (input) CHARACTER*1
73 = 'N': do not compute Q;
74 = 'I': Q is initialized to the unit matrix, and the
75 orthogonal matrix Q is returned;
76 = 'V': Q must contain an orthogonal matrix Q1 on entry,
77 and the product Q1*Q is returned.
78
79 COMPZ (input) CHARACTER*1
80 = 'N': do not compute Z;
81 = 'I': Z is initialized to the unit matrix, and the
82 orthogonal matrix Z is returned;
83 = 'V': Z must contain an orthogonal matrix Z1 on entry,
84 and the product Z1*Z is returned.
85
86 N (input) INTEGER
87 The order of the matrices A and B. N >= 0.
88
89 ILO (input) INTEGER
90 IHI (input) INTEGER
91 It is assumed that A is already upper triangular in rows and
92 columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set
93 by a previous call to DGGBAL; otherwise they should be set
94 to 1 and N respectively.
95 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
96
97 A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
98 On entry, the N-by-N general matrix to be reduced.
99 On exit, the upper triangle and the first subdiagonal of A
100 are overwritten with the upper Hessenberg matrix H, and the
101 rest is set to zero.
102
103 LDA (input) INTEGER
104 The leading dimension of the array A. LDA >= max(1,N).
105
106 B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
107 On entry, the N-by-N upper triangular matrix B.
108 On exit, the upper triangular matrix T = Q' B Z. The
109 elements below the diagonal are set to zero.
110
111 LDB (input) INTEGER
112 The leading dimension of the array B. LDB >= max(1,N).
113
114 Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
115 If COMPQ='N': Q is not referenced.
116 If COMPQ='I': on entry, Q need not be set, and on exit it
117 contains the orthogonal matrix Q, where Q'
118 is the product of the Givens transformations
119 which are applied to A and B on the left.
120 If COMPQ='V': on entry, Q must contain an orthogonal matrix
121 Q1, and on exit this is overwritten by Q1*Q.
122
123 LDQ (input) INTEGER
124 The leading dimension of the array Q.
125 LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
126
127 Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
128 If COMPZ='N': Z is not referenced.
129 If COMPZ='I': on entry, Z need not be set, and on exit it
130 contains the orthogonal matrix Z, which is
131 the product of the Givens transformations
132 which are applied to A and B on the right.
133 If COMPZ='V': on entry, Z must contain an orthogonal matrix
134 Z1, and on exit this is overwritten by Z1*Z.
135
136 LDZ (input) INTEGER
137 The leading dimension of the array Z.
138 LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
139
140 INFO (output) INTEGER
141 = 0: successful exit.
142 < 0: if INFO = -i, the i-th argument had an illegal value.
143
144 Further Details
145 ===============
146
147 This routine reduces A to Hessenberg and B to triangular form by
148 an unblocked reduction, as described in _Matrix_Computations_,
149 by Golub and Van Loan (Johns Hopkins Press.)
150
151 =====================================================================
152
153
154 Decode COMPQ
155
156 Parameter adjustments */
157 /* Table of constant values */
158 Treal c_b10 = 0.;
159 Treal c_b11 = 1.;
160 integer c__1 = 1;
161
162 /* System generated locals */
163 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
164 z_offset, i__1, i__2, i__3;
165 /* Local variables */
166 integer jcol;
167 Treal temp;
168 integer jrow;
169 Treal c__, s;
170 integer icompq, icompz;
171 logical ilq, ilz;
172#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
173#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
174#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
175#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
176
177
178 a_dim1 = *lda;
179 a_offset = 1 + a_dim1 * 1;
180 a -= a_offset;
181 b_dim1 = *ldb;
182 b_offset = 1 + b_dim1 * 1;
183 b -= b_offset;
184 q_dim1 = *ldq;
185 q_offset = 1 + q_dim1 * 1;
186 q -= q_offset;
187 z_dim1 = *ldz;
188 z_offset = 1 + z_dim1 * 1;
189 z__ -= z_offset;
190
191 /* Initialization added by Elias to get rid of compiler warnings. */
192 ilq = ilz = 0;
193 /* Function Body */
194 if (template_blas_lsame(compq, "N")) {
195 ilq = FALSE_;
196 icompq = 1;
197 } else if (template_blas_lsame(compq, "V")) {
198 ilq = TRUE_;
199 icompq = 2;
200 } else if (template_blas_lsame(compq, "I")) {
201 ilq = TRUE_;
202 icompq = 3;
203 } else {
204 icompq = 0;
205 }
206
207/* Decode COMPZ */
208
209 if (template_blas_lsame(compz, "N")) {
210 ilz = FALSE_;
211 icompz = 1;
212 } else if (template_blas_lsame(compz, "V")) {
213 ilz = TRUE_;
214 icompz = 2;
215 } else if (template_blas_lsame(compz, "I")) {
216 ilz = TRUE_;
217 icompz = 3;
218 } else {
219 icompz = 0;
220 }
221
222/* Test the input parameters. */
223
224 *info = 0;
225 if (icompq <= 0) {
226 *info = -1;
227 } else if (icompz <= 0) {
228 *info = -2;
229 } else if (*n < 0) {
230 *info = -3;
231 } else if (*ilo < 1) {
232 *info = -4;
233 } else if (*ihi > *n || *ihi < *ilo - 1) {
234 *info = -5;
235 } else if (*lda < maxMACRO(1,*n)) {
236 *info = -7;
237 } else if (*ldb < maxMACRO(1,*n)) {
238 *info = -9;
239 } else if ( ( ilq && *ldq < *n ) || *ldq < 1) {
240 *info = -11;
241 } else if ( ( ilz && *ldz < *n ) || *ldz < 1) {
242 *info = -13;
243 }
244 if (*info != 0) {
245 i__1 = -(*info);
246 template_blas_erbla("GGHRD ", &i__1);
247 return 0;
248 }
249
250/* Initialize Q and Z if desired. */
251
252 if (icompq == 3) {
253 template_lapack_laset("Full", n, n, &c_b10, &c_b11, &q[q_offset], ldq);
254 }
255 if (icompz == 3) {
256 template_lapack_laset("Full", n, n, &c_b10, &c_b11, &z__[z_offset], ldz);
257 }
258
259/* Quick return if possible */
260
261 if (*n <= 1) {
262 return 0;
263 }
264
265/* Zero out lower triangle of B */
266
267 i__1 = *n - 1;
268 for (jcol = 1; jcol <= i__1; ++jcol) {
269 i__2 = *n;
270 for (jrow = jcol + 1; jrow <= i__2; ++jrow) {
271 b_ref(jrow, jcol) = 0.;
272/* L10: */
273 }
274/* L20: */
275 }
276
277/* Reduce A and B */
278
279 i__1 = *ihi - 2;
280 for (jcol = *ilo; jcol <= i__1; ++jcol) {
281
282 i__2 = jcol + 2;
283 for (jrow = *ihi; jrow >= i__2; --jrow) {
284
285/* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */
286
287 temp = a_ref(jrow - 1, jcol);
288 template_lapack_lartg(&temp, &a_ref(jrow, jcol), &c__, &s, &a_ref(jrow - 1,
289 jcol));
290 a_ref(jrow, jcol) = 0.;
291 i__3 = *n - jcol;
292 template_blas_rot(&i__3, &a_ref(jrow - 1, jcol + 1), lda, &a_ref(jrow, jcol +
293 1), lda, &c__, &s);
294 i__3 = *n + 2 - jrow;
295 template_blas_rot(&i__3, &b_ref(jrow - 1, jrow - 1), ldb, &b_ref(jrow, jrow -
296 1), ldb, &c__, &s);
297 if (ilq) {
298 template_blas_rot(n, &q_ref(1, jrow - 1), &c__1, &q_ref(1, jrow), &c__1, &
299 c__, &s);
300 }
301
302/* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */
303
304 temp = b_ref(jrow, jrow);
305 template_lapack_lartg(&temp, &b_ref(jrow, jrow - 1), &c__, &s, &b_ref(jrow,
306 jrow));
307 b_ref(jrow, jrow - 1) = 0.;
308 template_blas_rot(ihi, &a_ref(1, jrow), &c__1, &a_ref(1, jrow - 1), &c__1, &
309 c__, &s);
310 i__3 = jrow - 1;
311 template_blas_rot(&i__3, &b_ref(1, jrow), &c__1, &b_ref(1, jrow - 1), &c__1, &
312 c__, &s);
313 if (ilz) {
314 template_blas_rot(n, &z___ref(1, jrow), &c__1, &z___ref(1, jrow - 1), &
315 c__1, &c__, &s);
316 }
317/* L30: */
318 }
319/* L40: */
320 }
321
322 return 0;
323
324/* End of DGGHRD */
325
326} /* dgghrd_ */
327
328#undef z___ref
329#undef q_ref
330#undef b_ref
331#undef a_ref
332
333
334#endif
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
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
int template_blas_rot(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy, const Treal *c__, const Treal *s)
Definition: template_blas_rot.h:42
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)
#define z___ref(a_1, a_2)
int template_lapack_gghrd(const char *compq, const char *compz, const integer *n, const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b, const integer *ldb, Treal *q, const integer *ldq, Treal *z__, const integer *ldz, integer *info)
Definition: template_lapack_gghrd.h:42
#define q_ref(a_1, a_2)
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:42
int template_lapack_laset(const char *uplo, const integer *m, const integer *n, const Treal *alpha, const Treal *beta, Treal *a, const integer *lda)
Definition: template_lapack_laset.h:42