ergo
template_lapack_orgtr.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_ORGTR_HEADER
38#define TEMPLATE_LAPACK_ORGTR_HEADER
39
40
41template<class Treal>
42int 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_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
int ftnlen
Definition template_blas_common.h:42
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
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
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
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 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
#define a_ref(a_1, a_2)