ergo
template_lapack_gghrd.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_GGHRD_HEADER
36 #define TEMPLATE_LAPACK_GGHRD_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_gghrd(const char *compq, const char *compz, const integer *n, const integer *
41  ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b,
42  const integer *ldb, Treal *q, const integer *ldq, Treal *z__, const integer *
43  ldz, integer *info)
44 {
45 /* -- LAPACK routine (version 3.0) --
46  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47  Courant Institute, Argonne National Lab, and Rice University
48  September 30, 1994
49 
50 
51  Purpose
52  =======
53 
54  DGGHRD reduces a pair of real matrices (A,B) to generalized upper
55  Hessenberg form using orthogonal transformations, where A is a
56  general matrix and B is upper triangular: Q' * A * Z = H and
57  Q' * B * Z = T, where H is upper Hessenberg, T is upper triangular,
58  and Q and Z are orthogonal, and ' means transpose.
59 
60  The orthogonal matrices Q and Z are determined as products of Givens
61  rotations. They may either be formed explicitly, or they may be
62  postmultiplied into input matrices Q1 and Z1, so that
63 
64  Q1 * A * Z1' = (Q1*Q) * H * (Z1*Z)'
65  Q1 * B * Z1' = (Q1*Q) * T * (Z1*Z)'
66 
67  Arguments
68  =========
69 
70  COMPQ (input) CHARACTER*1
71  = 'N': do not compute Q;
72  = 'I': Q is initialized to the unit matrix, and the
73  orthogonal matrix Q is returned;
74  = 'V': Q must contain an orthogonal matrix Q1 on entry,
75  and the product Q1*Q is returned.
76 
77  COMPZ (input) CHARACTER*1
78  = 'N': do not compute Z;
79  = 'I': Z is initialized to the unit matrix, and the
80  orthogonal matrix Z is returned;
81  = 'V': Z must contain an orthogonal matrix Z1 on entry,
82  and the product Z1*Z is returned.
83 
84  N (input) INTEGER
85  The order of the matrices A and B. N >= 0.
86 
87  ILO (input) INTEGER
88  IHI (input) INTEGER
89  It is assumed that A is already upper triangular in rows and
90  columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set
91  by a previous call to DGGBAL; otherwise they should be set
92  to 1 and N respectively.
93  1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
94 
95  A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
96  On entry, the N-by-N general matrix to be reduced.
97  On exit, the upper triangle and the first subdiagonal of A
98  are overwritten with the upper Hessenberg matrix H, and the
99  rest is set to zero.
100 
101  LDA (input) INTEGER
102  The leading dimension of the array A. LDA >= max(1,N).
103 
104  B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
105  On entry, the N-by-N upper triangular matrix B.
106  On exit, the upper triangular matrix T = Q' B Z. The
107  elements below the diagonal are set to zero.
108 
109  LDB (input) INTEGER
110  The leading dimension of the array B. LDB >= max(1,N).
111 
112  Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
113  If COMPQ='N': Q is not referenced.
114  If COMPQ='I': on entry, Q need not be set, and on exit it
115  contains the orthogonal matrix Q, where Q'
116  is the product of the Givens transformations
117  which are applied to A and B on the left.
118  If COMPQ='V': on entry, Q must contain an orthogonal matrix
119  Q1, and on exit this is overwritten by Q1*Q.
120 
121  LDQ (input) INTEGER
122  The leading dimension of the array Q.
123  LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
124 
125  Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
126  If COMPZ='N': Z is not referenced.
127  If COMPZ='I': on entry, Z need not be set, and on exit it
128  contains the orthogonal matrix Z, which is
129  the product of the Givens transformations
130  which are applied to A and B on the right.
131  If COMPZ='V': on entry, Z must contain an orthogonal matrix
132  Z1, and on exit this is overwritten by Z1*Z.
133 
134  LDZ (input) INTEGER
135  The leading dimension of the array Z.
136  LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
137 
138  INFO (output) INTEGER
139  = 0: successful exit.
140  < 0: if INFO = -i, the i-th argument had an illegal value.
141 
142  Further Details
143  ===============
144 
145  This routine reduces A to Hessenberg and B to triangular form by
146  an unblocked reduction, as described in _Matrix_Computations_,
147  by Golub and Van Loan (Johns Hopkins Press.)
148 
149  =====================================================================
150 
151 
152  Decode COMPQ
153 
154  Parameter adjustments */
155  /* Table of constant values */
156  Treal c_b10 = 0.;
157  Treal c_b11 = 1.;
158  integer c__1 = 1;
159 
160  /* System generated locals */
161  integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
162  z_offset, i__1, i__2, i__3;
163  /* Local variables */
164  integer jcol;
165  Treal temp;
166  integer jrow;
167  Treal c__, s;
168  integer icompq, icompz;
169  logical ilq, ilz;
170 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
171 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
172 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
173 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
174 
175 
176  a_dim1 = *lda;
177  a_offset = 1 + a_dim1 * 1;
178  a -= a_offset;
179  b_dim1 = *ldb;
180  b_offset = 1 + b_dim1 * 1;
181  b -= b_offset;
182  q_dim1 = *ldq;
183  q_offset = 1 + q_dim1 * 1;
184  q -= q_offset;
185  z_dim1 = *ldz;
186  z_offset = 1 + z_dim1 * 1;
187  z__ -= z_offset;
188 
189  /* Initialization added by Elias to get rid of compiler warnings. */
190  ilq = ilz = 0;
191  /* Function Body */
192  if (template_blas_lsame(compq, "N")) {
193  ilq = FALSE_;
194  icompq = 1;
195  } else if (template_blas_lsame(compq, "V")) {
196  ilq = TRUE_;
197  icompq = 2;
198  } else if (template_blas_lsame(compq, "I")) {
199  ilq = TRUE_;
200  icompq = 3;
201  } else {
202  icompq = 0;
203  }
204 
205 /* Decode COMPZ */
206 
207  if (template_blas_lsame(compz, "N")) {
208  ilz = FALSE_;
209  icompz = 1;
210  } else if (template_blas_lsame(compz, "V")) {
211  ilz = TRUE_;
212  icompz = 2;
213  } else if (template_blas_lsame(compz, "I")) {
214  ilz = TRUE_;
215  icompz = 3;
216  } else {
217  icompz = 0;
218  }
219 
220 /* Test the input parameters. */
221 
222  *info = 0;
223  if (icompq <= 0) {
224  *info = -1;
225  } else if (icompz <= 0) {
226  *info = -2;
227  } else if (*n < 0) {
228  *info = -3;
229  } else if (*ilo < 1) {
230  *info = -4;
231  } else if (*ihi > *n || *ihi < *ilo - 1) {
232  *info = -5;
233  } else if (*lda < maxMACRO(1,*n)) {
234  *info = -7;
235  } else if (*ldb < maxMACRO(1,*n)) {
236  *info = -9;
237  } else if ( ( ilq && *ldq < *n ) || *ldq < 1) {
238  *info = -11;
239  } else if ( ( ilz && *ldz < *n ) || *ldz < 1) {
240  *info = -13;
241  }
242  if (*info != 0) {
243  i__1 = -(*info);
244  template_blas_erbla("GGHRD ", &i__1);
245  return 0;
246  }
247 
248 /* Initialize Q and Z if desired. */
249 
250  if (icompq == 3) {
251  template_lapack_laset("Full", n, n, &c_b10, &c_b11, &q[q_offset], ldq);
252  }
253  if (icompz == 3) {
254  template_lapack_laset("Full", n, n, &c_b10, &c_b11, &z__[z_offset], ldz);
255  }
256 
257 /* Quick return if possible */
258 
259  if (*n <= 1) {
260  return 0;
261  }
262 
263 /* Zero out lower triangle of B */
264 
265  i__1 = *n - 1;
266  for (jcol = 1; jcol <= i__1; ++jcol) {
267  i__2 = *n;
268  for (jrow = jcol + 1; jrow <= i__2; ++jrow) {
269  b_ref(jrow, jcol) = 0.;
270 /* L10: */
271  }
272 /* L20: */
273  }
274 
275 /* Reduce A and B */
276 
277  i__1 = *ihi - 2;
278  for (jcol = *ilo; jcol <= i__1; ++jcol) {
279 
280  i__2 = jcol + 2;
281  for (jrow = *ihi; jrow >= i__2; --jrow) {
282 
283 /* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */
284 
285  temp = a_ref(jrow - 1, jcol);
286  template_lapack_lartg(&temp, &a_ref(jrow, jcol), &c__, &s, &a_ref(jrow - 1,
287  jcol));
288  a_ref(jrow, jcol) = 0.;
289  i__3 = *n - jcol;
290  template_blas_rot(&i__3, &a_ref(jrow - 1, jcol + 1), lda, &a_ref(jrow, jcol +
291  1), lda, &c__, &s);
292  i__3 = *n + 2 - jrow;
293  template_blas_rot(&i__3, &b_ref(jrow - 1, jrow - 1), ldb, &b_ref(jrow, jrow -
294  1), ldb, &c__, &s);
295  if (ilq) {
296  template_blas_rot(n, &q_ref(1, jrow - 1), &c__1, &q_ref(1, jrow), &c__1, &
297  c__, &s);
298  }
299 
300 /* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */
301 
302  temp = b_ref(jrow, jrow);
303  template_lapack_lartg(&temp, &b_ref(jrow, jrow - 1), &c__, &s, &b_ref(jrow,
304  jrow));
305  b_ref(jrow, jrow - 1) = 0.;
306  template_blas_rot(ihi, &a_ref(1, jrow), &c__1, &a_ref(1, jrow - 1), &c__1, &
307  c__, &s);
308  i__3 = jrow - 1;
309  template_blas_rot(&i__3, &b_ref(1, jrow), &c__1, &b_ref(1, jrow - 1), &c__1, &
310  c__, &s);
311  if (ilz) {
312  template_blas_rot(n, &z___ref(1, jrow), &c__1, &z___ref(1, jrow - 1), &
313  c__1, &c__, &s);
314  }
315 /* L30: */
316  }
317 /* L40: */
318  }
319 
320  return 0;
321 
322 /* End of DGGHRD */
323 
324 } /* dgghrd_ */
325 
326 #undef z___ref
327 #undef q_ref
328 #undef b_ref
329 #undef a_ref
330 
331 
332 #endif
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:40
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
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:40
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:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.h:40
#define z___ref(a_1, a_2)
#define b_ref(a_1, a_2)
#define a_ref(a_1, a_2)
bool logical
Definition: template_blas_common.h:39
#define TRUE_
Definition: template_lapack_common.h:40
#define FALSE_
Definition: template_lapack_common.h:41
#define q_ref(a_1, a_2)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44