ergo
template_lapack_gghrd.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_GGHRD_HEADER
38 #define TEMPLATE_LAPACK_GGHRD_HEADER
39 
40 
41 template<class Treal>
42 int 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
z___ref
#define z___ref(a_1, a_2)
template_lapack_laset
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
template_blas_rot
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
logical
bool logical
Definition: template_blas_common.h:41
q_ref
#define q_ref(a_1, a_2)
template_lapack_lartg
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:42
template_lapack_gghrd
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
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
TRUE_
#define TRUE_
Definition: template_lapack_common.h:42
integer
int integer
Definition: template_blas_common.h:40
b_ref
#define b_ref(a_1, a_2)
FALSE_
#define FALSE_
Definition: template_lapack_common.h:43
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
a_ref
#define a_ref(a_1, a_2)