ergo
template_lapack_orgtr.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_ORGTR_HEADER
36 #define TEMPLATE_LAPACK_ORGTR_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_orgtr(const char *uplo, const integer *n, Treal *a, const integer *
41  lda, const Treal *tau, Treal *work, const integer *lwork, integer *info)
42 {
43 /* -- LAPACK routine (version 3.0) --
44  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
45  Courant Institute, Argonne National Lab, and Rice University
46  June 30, 1999
47 
48 
49  Purpose
50  =======
51 
52  DORGTR generates a real orthogonal matrix Q which is defined as the
53  product of n-1 elementary reflectors of order N, as returned by
54  DSYTRD:
55 
56  if UPLO = 'U', Q = H(n-1) . . . H(2) H(1),
57 
58  if UPLO = 'L', Q = H(1) H(2) . . . H(n-1).
59 
60  Arguments
61  =========
62 
63  UPLO (input) CHARACTER*1
64  = 'U': Upper triangle of A contains elementary reflectors
65  from DSYTRD;
66  = 'L': Lower triangle of A contains elementary reflectors
67  from DSYTRD.
68 
69  N (input) INTEGER
70  The order of the matrix Q. N >= 0.
71 
72  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
73  On entry, the vectors which define the elementary reflectors,
74  as returned by DSYTRD.
75  On exit, the N-by-N orthogonal matrix Q.
76 
77  LDA (input) INTEGER
78  The leading dimension of the array A. LDA >= max(1,N).
79 
80  TAU (input) DOUBLE PRECISION array, dimension (N-1)
81  TAU(i) must contain the scalar factor of the elementary
82  reflector H(i), as returned by DSYTRD.
83 
84  WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
85  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
86 
87  LWORK (input) INTEGER
88  The dimension of the array WORK. LWORK >= max(1,N-1).
89  For optimum performance LWORK >= (N-1)*NB, where NB is
90  the optimal blocksize.
91 
92  If LWORK = -1, then a workspace query is assumed; the routine
93  only calculates the optimal size of the WORK array, returns
94  this value as the first entry of the WORK array, and no error
95  message related to LWORK is issued by XERBLA.
96 
97  INFO (output) INTEGER
98  = 0: successful exit
99  < 0: if INFO = -i, the i-th argument had an illegal value
100 
101  =====================================================================
102 
103 
104  Test the input arguments
105 
106  Parameter adjustments */
107  /* Table of constant values */
108  integer c__1 = 1;
109  integer c_n1 = -1;
110 
111  /* System generated locals */
112  integer a_dim1, a_offset, i__1, i__2, i__3;
113  /* Local variables */
114  integer i__, j;
115  integer iinfo;
116  logical upper;
117  integer nb;
118  integer lwkopt;
119  logical lquery;
120 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
121 
122 
123  a_dim1 = *lda;
124  a_offset = 1 + a_dim1 * 1;
125  a -= a_offset;
126  --tau;
127  --work;
128 
129  /* Initialization added by Elias to get rid of compiler warnings. */
130  lwkopt = 0;
131  /* Function Body */
132  *info = 0;
133  lquery = *lwork == -1;
134  upper = template_blas_lsame(uplo, "U");
135  if (! upper && ! template_blas_lsame(uplo, "L")) {
136  *info = -1;
137  } else if (*n < 0) {
138  *info = -2;
139  } else if (*lda < maxMACRO(1,*n)) {
140  *info = -4;
141  } else /* if(complicated condition) */ {
142 /* Computing MAX */
143  i__1 = 1, i__2 = *n - 1;
144  if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
145  *info = -7;
146  }
147  }
148 
149  if (*info == 0) {
150  if (upper) {
151  i__1 = *n - 1;
152  i__2 = *n - 1;
153  i__3 = *n - 1;
154  nb = template_lapack_ilaenv(&c__1, "DORGQL", " ", &i__1, &i__2, &i__3, &c_n1, (
155  ftnlen)6, (ftnlen)1);
156  } else {
157  i__1 = *n - 1;
158  i__2 = *n - 1;
159  i__3 = *n - 1;
160  nb = template_lapack_ilaenv(&c__1, "DORGQR", " ", &i__1, &i__2, &i__3, &c_n1, (
161  ftnlen)6, (ftnlen)1);
162  }
163 /* Computing MAX */
164  i__1 = 1, i__2 = *n - 1;
165  lwkopt = maxMACRO(i__1,i__2) * nb;
166  work[1] = (Treal) lwkopt;
167  }
168 
169  if (*info != 0) {
170  i__1 = -(*info);
171  template_blas_erbla("ORGTR ", &i__1);
172  return 0;
173  } else if (lquery) {
174  return 0;
175  }
176 
177 /* Quick return if possible */
178 
179  if (*n == 0) {
180  work[1] = 1.;
181  return 0;
182  }
183 
184  if (upper) {
185 
186 /* Q was determined by a call to DSYTRD with UPLO = 'U'
187 
188  Shift the vectors which define the elementary reflectors one
189  column to the left, and set the last row and column of Q to
190  those of the unit matrix */
191 
192  i__1 = *n - 1;
193  for (j = 1; j <= i__1; ++j) {
194  i__2 = j - 1;
195  for (i__ = 1; i__ <= i__2; ++i__) {
196  a_ref(i__, j) = a_ref(i__, j + 1);
197 /* L10: */
198  }
199  a_ref(*n, j) = 0.;
200 /* L20: */
201  }
202  i__1 = *n - 1;
203  for (i__ = 1; i__ <= i__1; ++i__) {
204  a_ref(i__, *n) = 0.;
205 /* L30: */
206  }
207  a_ref(*n, *n) = 1.;
208 
209 /* Generate Q(1:n-1,1:n-1) */
210 
211  i__1 = *n - 1;
212  i__2 = *n - 1;
213  i__3 = *n - 1;
214  template_lapack_orgql(&i__1, &i__2, &i__3, &a[a_offset], lda, &tau[1], &work[1],
215  lwork, &iinfo);
216 
217  } else {
218 
219 /* Q was determined by a call to DSYTRD with UPLO = 'L'.
220 
221  Shift the vectors which define the elementary reflectors one
222  column to the right, and set the first row and column of Q to
223  those of the unit matrix */
224 
225  for (j = *n; j >= 2; --j) {
226  a_ref(1, j) = 0.;
227  i__1 = *n;
228  for (i__ = j + 1; i__ <= i__1; ++i__) {
229  a_ref(i__, j) = a_ref(i__, j - 1);
230 /* L40: */
231  }
232 /* L50: */
233  }
234  a_ref(1, 1) = 1.;
235  i__1 = *n;
236  for (i__ = 2; i__ <= i__1; ++i__) {
237  a_ref(i__, 1) = 0.;
238 /* L60: */
239  }
240  if (*n > 1) {
241 
242 /* Generate Q(2:n,2:n) */
243 
244  i__1 = *n - 1;
245  i__2 = *n - 1;
246  i__3 = *n - 1;
248  &i__1,
249  &i__2,
250  &i__3,
251  &a_ref(2, 2),
252  lda,
253  &tau[1],
254  &work[1],
255  lwork,
256  &iinfo
257  );
258  }
259  }
260  work[1] = (Treal) lwkopt;
261  return 0;
262 
263 /* End of DORGTR */
264 
265 } /* dorgtr_ */
266 
267 #undef a_ref
268 
269 
270 #endif
int template_lapack_orgqr(const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_orgqr.h:41
int integer
Definition: template_blas_common.h:38
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition: template_lapack_common.cc:279
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_orgtr(const char *uplo, const integer *n, Treal *a, const integer *lda, const Treal *tau, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_orgtr.h:40
bool logical
Definition: template_blas_common.h:39
#define a_ref(a_1, a_2)
int ftnlen
Definition: template_blas_common.h:40
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44
int template_lapack_orgql(const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_orgql.h:40