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