ergo
template_lapack_org2l.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_ORG2L_HEADER
38#define TEMPLATE_LAPACK_ORG2L_HEADER
39
40
41template<class Treal>
42int template_lapack_org2l(const integer *m, const integer *n, const integer *k, Treal *
43 a, const integer *lda, const Treal *tau, Treal *work, 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 February 29, 1992
49
50
51 Purpose
52 =======
53
54 DORG2L generates an m by n real matrix Q with orthonormal columns,
55 which is defined as the last n columns of a product of k elementary
56 reflectors of order m
57
58 Q = H(k) . . . H(2) H(1)
59
60 as returned by DGEQLF.
61
62 Arguments
63 =========
64
65 M (input) INTEGER
66 The number of rows of the matrix Q. M >= 0.
67
68 N (input) INTEGER
69 The number of columns of the matrix Q. M >= N >= 0.
70
71 K (input) INTEGER
72 The number of elementary reflectors whose product defines the
73 matrix Q. N >= K >= 0.
74
75 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
76 On entry, the (n-k+i)-th column must contain the vector which
77 defines the elementary reflector H(i), for i = 1,2,...,k, as
78 returned by DGEQLF in the last k columns of its array
79 argument A.
80 On exit, the m by n matrix Q.
81
82 LDA (input) INTEGER
83 The first dimension of the array A. LDA >= max(1,M).
84
85 TAU (input) DOUBLE PRECISION array, dimension (K)
86 TAU(i) must contain the scalar factor of the elementary
87 reflector H(i), as returned by DGEQLF.
88
89 WORK (workspace) DOUBLE PRECISION array, dimension (N)
90
91 INFO (output) INTEGER
92 = 0: successful exit
93 < 0: if INFO = -i, the i-th argument has an illegal value
94
95 =====================================================================
96
97
98 Test the input arguments
99
100 Parameter adjustments */
101 /* Table of constant values */
102 integer c__1 = 1;
103
104 /* System generated locals */
105 integer a_dim1, a_offset, i__1, i__2, i__3;
106 Treal d__1;
107 /* Local variables */
108 integer i__, j, l;
109 integer ii;
110#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
111
112
113 a_dim1 = *lda;
114 a_offset = 1 + a_dim1 * 1;
115 a -= a_offset;
116 --tau;
117 --work;
118
119 /* Function Body */
120 *info = 0;
121 if (*m < 0) {
122 *info = -1;
123 } else if (*n < 0 || *n > *m) {
124 *info = -2;
125 } else if (*k < 0 || *k > *n) {
126 *info = -3;
127 } else if (*lda < maxMACRO(1,*m)) {
128 *info = -5;
129 }
130 if (*info != 0) {
131 i__1 = -(*info);
132 template_blas_erbla("ORG2L ", &i__1);
133 return 0;
134 }
135
136/* Quick return if possible */
137
138 if (*n <= 0) {
139 return 0;
140 }
141
142/* Initialise columns 1:n-k to columns of the unit matrix */
143
144 i__1 = *n - *k;
145 for (j = 1; j <= i__1; ++j) {
146 i__2 = *m;
147 for (l = 1; l <= i__2; ++l) {
148 a_ref(l, j) = 0.;
149/* L10: */
150 }
151 a_ref(*m - *n + j, j) = 1.;
152/* L20: */
153 }
154
155 i__1 = *k;
156 for (i__ = 1; i__ <= i__1; ++i__) {
157 ii = *n - *k + i__;
158
159/* Apply H(i) to A(1:m-k+i,1:n-k+i) from the left */
160
161 a_ref(*m - *n + ii, ii) = 1.;
162 i__2 = *m - *n + ii;
163 i__3 = ii - 1;
164 template_lapack_larf("Left", &i__2, &i__3, &a_ref(1, ii), &c__1, &tau[i__], &a[
165 a_offset], lda, &work[1]);
166 i__2 = *m - *n + ii - 1;
167 d__1 = -tau[i__];
168 template_blas_scal(&i__2, &d__1, &a_ref(1, ii), &c__1);
169 a_ref(*m - *n + ii, ii) = 1. - tau[i__];
170
171/* Set A(m-k+i+1:m,n-k+i) to zero */
172
173 i__2 = *m;
174 for (l = *m - *n + ii + 1; l <= i__2; ++l) {
175 a_ref(l, ii) = 0.;
176/* L30: */
177 }
178/* L40: */
179 }
180 return 0;
181
182/* End of DORG2L */
183
184} /* dorg2l_ */
185
186#undef a_ref
187
188
189#endif
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
int integer
Definition template_blas_common.h:40
#define maxMACRO(a, b)
Definition template_blas_common.h:45
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
int template_lapack_larf(const char *side, const integer *m, const integer *n, const Treal *v, const integer *incv, const Treal *tau, Treal *c__, const integer *ldc, Treal *work)
Definition template_lapack_larf.h:42
#define a_ref(a_1, a_2)
int template_lapack_org2l(const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *work, integer *info)
Definition template_lapack_org2l.h:42