ergo
template_lapack_orm2r.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 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_ORM2R_HEADER
38#define TEMPLATE_LAPACK_ORM2R_HEADER
39
40
41template<class Treal>
42int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n,
43 const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *
44 c__, const integer *ldc, Treal *work, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 February 29, 1992
50
51
52 Purpose
53 =======
54
55 DORM2R overwrites the general real m by n matrix C with
56
57 Q * C if SIDE = 'L' and TRANS = 'N', or
58
59 Q'* C if SIDE = 'L' and TRANS = 'T', or
60
61 C * Q if SIDE = 'R' and TRANS = 'N', or
62
63 C * Q' if SIDE = 'R' and TRANS = 'T',
64
65 where Q is a real orthogonal matrix defined as the product of k
66 elementary reflectors
67
68 Q = H(1) H(2) . . . H(k)
69
70 as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n
71 if SIDE = 'R'.
72
73 Arguments
74 =========
75
76 SIDE (input) CHARACTER*1
77 = 'L': apply Q or Q' from the Left
78 = 'R': apply Q or Q' from the Right
79
80 TRANS (input) CHARACTER*1
81 = 'N': apply Q (No transpose)
82 = 'T': apply Q' (Transpose)
83
84 M (input) INTEGER
85 The number of rows of the matrix C. M >= 0.
86
87 N (input) INTEGER
88 The number of columns of the matrix C. N >= 0.
89
90 K (input) INTEGER
91 The number of elementary reflectors whose product defines
92 the matrix Q.
93 If SIDE = 'L', M >= K >= 0;
94 if SIDE = 'R', N >= K >= 0.
95
96 A (input) DOUBLE PRECISION array, dimension (LDA,K)
97 The i-th column must contain the vector which defines the
98 elementary reflector H(i), for i = 1,2,...,k, as returned by
99 DGEQRF in the first k columns of its array argument A.
100 A is modified by the routine but restored on exit.
101
102 LDA (input) INTEGER
103 The leading dimension of the array A.
104 If SIDE = 'L', LDA >= max(1,M);
105 if SIDE = 'R', LDA >= max(1,N).
106
107 TAU (input) DOUBLE PRECISION array, dimension (K)
108 TAU(i) must contain the scalar factor of the elementary
109 reflector H(i), as returned by DGEQRF.
110
111 C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
112 On entry, the m by n matrix C.
113 On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
114
115 LDC (input) INTEGER
116 The leading dimension of the array C. LDC >= max(1,M).
117
118 WORK (workspace) DOUBLE PRECISION array, dimension
119 (N) if SIDE = 'L',
120 (M) if SIDE = 'R'
121
122 INFO (output) INTEGER
123 = 0: successful exit
124 < 0: if INFO = -i, the i-th argument had an illegal value
125
126 =====================================================================
127
128
129 Test the input arguments
130
131 Parameter adjustments */
132 /* Table of constant values */
133 integer c__1 = 1;
134
135 /* System generated locals */
136 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
137 /* Local variables */
138 logical left;
139 integer i__;
140 integer i1, i2, i3, ic, jc, mi, ni, nq;
141 logical notran;
142 Treal aii;
143#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
144#define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
145
146
147 a_dim1 = *lda;
148 a_offset = 1 + a_dim1 * 1;
149 a -= a_offset;
150 --tau;
151 c_dim1 = *ldc;
152 c_offset = 1 + c_dim1 * 1;
153 c__ -= c_offset;
154 --work;
155
156 /* Function Body */
157 *info = 0;
158 left = template_blas_lsame(side, "L");
159 notran = template_blas_lsame(trans, "N");
160
161/* NQ is the order of Q */
162
163 if (left) {
164 nq = *m;
165 } else {
166 nq = *n;
167 }
168 if (! left && ! template_blas_lsame(side, "R")) {
169 *info = -1;
170 } else if (! notran && ! template_blas_lsame(trans, "T")) {
171 *info = -2;
172 } else if (*m < 0) {
173 *info = -3;
174 } else if (*n < 0) {
175 *info = -4;
176 } else if (*k < 0 || *k > nq) {
177 *info = -5;
178 } else if (*lda < maxMACRO(1,nq)) {
179 *info = -7;
180 } else if (*ldc < maxMACRO(1,*m)) {
181 *info = -10;
182 }
183 if (*info != 0) {
184 i__1 = -(*info);
185 template_blas_erbla("ORM2R ", &i__1);
186 return 0;
187 }
188
189/* Quick return if possible */
190
191 if (*m == 0 || *n == 0 || *k == 0) {
192 return 0;
193 }
194
195 if ( ( left && ! notran ) || ( ! left && notran ) ) {
196 i1 = 1;
197 i2 = *k;
198 i3 = 1;
199 } else {
200 i1 = *k;
201 i2 = 1;
202 i3 = -1;
203 }
204
205 if (left) {
206 ni = *n;
207 jc = 1;
208 } else {
209 mi = *m;
210 ic = 1;
211 }
212
213 i__1 = i2;
214 i__2 = i3;
215 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
216 if (left) {
217
218/* H(i) is applied to C(i:m,1:n) */
219
220 mi = *m - i__ + 1;
221 ic = i__;
222 } else {
223
224/* H(i) is applied to C(1:m,i:n) */
225
226 ni = *n - i__ + 1;
227 jc = i__;
228 }
229
230/* Apply H(i) */
231
232 aii = a_ref(i__, i__);
233 a_ref(i__, i__) = 1.;
234 template_lapack_larf(side, &mi, &ni, &a_ref(i__, i__), &c__1, &tau[i__], &c___ref(
235 ic, jc), ldc, &work[1]);
236 a_ref(i__, i__) = aii;
237/* L10: */
238 }
239 return 0;
240
241/* End of DORM2R */
242
243} /* dorm2r_ */
244
245#undef c___ref
246#undef a_ref
247
248
249#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
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
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 c___ref(a_1, a_2)
#define a_ref(a_1, a_2)
int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *c__, const integer *ldc, Treal *work, integer *info)
Definition: template_lapack_orm2r.h:42