ergo
template_lapack_lasr.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_LASR_HEADER
38#define TEMPLATE_LAPACK_LASR_HEADER
39
40
41template<class Treal>
42int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m,
43 const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer *
44 lda)
45{
46/* -- LAPACK auxiliary routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 October 31, 1992
50
51
52 Purpose
53 =======
54
55 DLASR performs the transformation
56
57 A := P*A, when SIDE = 'L' or 'l' ( Left-hand side )
58
59 A := A*P', when SIDE = 'R' or 'r' ( Right-hand side )
60
61 where A is an m by n real matrix and P is an orthogonal matrix,
62 consisting of a sequence of plane rotations determined by the
63 parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'
64 and z = n when SIDE = 'R' or 'r' ):
65
66 When DIRECT = 'F' or 'f' ( Forward sequence ) then
67
68 P = P( z - 1 )*...*P( 2 )*P( 1 ),
69
70 and when DIRECT = 'B' or 'b' ( Backward sequence ) then
71
72 P = P( 1 )*P( 2 )*...*P( z - 1 ),
73
74 where P( k ) is a plane rotation matrix for the following planes:
75
76 when PIVOT = 'V' or 'v' ( Variable pivot ),
77 the plane ( k, k + 1 )
78
79 when PIVOT = 'T' or 't' ( Top pivot ),
80 the plane ( 1, k + 1 )
81
82 when PIVOT = 'B' or 'b' ( Bottom pivot ),
83 the plane ( k, z )
84
85 c( k ) and s( k ) must contain the cosine and sine that define the
86 matrix P( k ). The two by two plane rotation part of the matrix
87 P( k ), R( k ), is assumed to be of the form
88
89 R( k ) = ( c( k ) s( k ) ).
90 ( -s( k ) c( k ) )
91
92 This version vectorises across rows of the array A when SIDE = 'L'.
93
94 Arguments
95 =========
96
97 SIDE (input) CHARACTER*1
98 Specifies whether the plane rotation matrix P is applied to
99 A on the left or the right.
100 = 'L': Left, compute A := P*A
101 = 'R': Right, compute A:= A*P'
102
103 DIRECT (input) CHARACTER*1
104 Specifies whether P is a forward or backward sequence of
105 plane rotations.
106 = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )
107 = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )
108
109 PIVOT (input) CHARACTER*1
110 Specifies the plane for which P(k) is a plane rotation
111 matrix.
112 = 'V': Variable pivot, the plane (k,k+1)
113 = 'T': Top pivot, the plane (1,k+1)
114 = 'B': Bottom pivot, the plane (k,z)
115
116 M (input) INTEGER
117 The number of rows of the matrix A. If m <= 1, an immediate
118 return is effected.
119
120 N (input) INTEGER
121 The number of columns of the matrix A. If n <= 1, an
122 immediate return is effected.
123
124 C, S (input) DOUBLE PRECISION arrays, dimension
125 (M-1) if SIDE = 'L'
126 (N-1) if SIDE = 'R'
127 c(k) and s(k) contain the cosine and sine that define the
128 matrix P(k). The two by two plane rotation part of the
129 matrix P(k), R(k), is assumed to be of the form
130 R( k ) = ( c( k ) s( k ) ).
131 ( -s( k ) c( k ) )
132
133 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
134 The m by n matrix A. On exit, A is overwritten by P*A if
135 SIDE = 'R' or by A*P' if SIDE = 'L'.
136
137 LDA (input) INTEGER
138 The leading dimension of the array A. LDA >= max(1,M).
139
140 =====================================================================
141
142
143 Test the input parameters
144
145 Parameter adjustments */
146 /* System generated locals */
147 integer a_dim1, a_offset, i__1, i__2;
148 /* Local variables */
149 integer info;
150 Treal temp;
151 integer i__, j;
152 Treal ctemp, stemp;
153#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
154
155 --c__;
156 --s;
157 a_dim1 = *lda;
158 a_offset = 1 + a_dim1 * 1;
159 a -= a_offset;
160
161 /* Function Body */
162 info = 0;
163 if (! (template_blas_lsame(side, "L") || template_blas_lsame(side, "R"))) {
164 info = 1;
165 } else if (! (template_blas_lsame(pivot, "V") || template_blas_lsame(pivot,
166 "T") || template_blas_lsame(pivot, "B"))) {
167 info = 2;
168 } else if (! (template_blas_lsame(direct, "F") || template_blas_lsame(direct,
169 "B"))) {
170 info = 3;
171 } else if (*m < 0) {
172 info = 4;
173 } else if (*n < 0) {
174 info = 5;
175 } else if (*lda < maxMACRO(1,*m)) {
176 info = 9;
177 }
178 if (info != 0) {
179 template_blas_erbla("LASR ", &info);
180 return 0;
181 }
182
183/* Quick return if possible */
184
185 if (*m == 0 || *n == 0) {
186 return 0;
187 }
188 if (template_blas_lsame(side, "L")) {
189
190/* Form P * A */
191
192 if (template_blas_lsame(pivot, "V")) {
193 if (template_blas_lsame(direct, "F")) {
194 i__1 = *m - 1;
195 for (j = 1; j <= i__1; ++j) {
196 ctemp = c__[j];
197 stemp = s[j];
198 if (ctemp != 1. || stemp != 0.) {
199 i__2 = *n;
200 for (i__ = 1; i__ <= i__2; ++i__) {
201 temp = a_ref(j + 1, i__);
202 a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
203 j, i__);
204 a_ref(j, i__) = stemp * temp + ctemp * a_ref(j,
205 i__);
206/* L10: */
207 }
208 }
209/* L20: */
210 }
211 } else if (template_blas_lsame(direct, "B")) {
212 for (j = *m - 1; j >= 1; --j) {
213 ctemp = c__[j];
214 stemp = s[j];
215 if (ctemp != 1. || stemp != 0.) {
216 i__1 = *n;
217 for (i__ = 1; i__ <= i__1; ++i__) {
218 temp = a_ref(j + 1, i__);
219 a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
220 j, i__);
221 a_ref(j, i__) = stemp * temp + ctemp * a_ref(j,
222 i__);
223/* L30: */
224 }
225 }
226/* L40: */
227 }
228 }
229 } else if (template_blas_lsame(pivot, "T")) {
230 if (template_blas_lsame(direct, "F")) {
231 i__1 = *m;
232 for (j = 2; j <= i__1; ++j) {
233 ctemp = c__[j - 1];
234 stemp = s[j - 1];
235 if (ctemp != 1. || stemp != 0.) {
236 i__2 = *n;
237 for (i__ = 1; i__ <= i__2; ++i__) {
238 temp = a_ref(j, i__);
239 a_ref(j, i__) = ctemp * temp - stemp * a_ref(1,
240 i__);
241 a_ref(1, i__) = stemp * temp + ctemp * a_ref(1,
242 i__);
243/* L50: */
244 }
245 }
246/* L60: */
247 }
248 } else if (template_blas_lsame(direct, "B")) {
249 for (j = *m; j >= 2; --j) {
250 ctemp = c__[j - 1];
251 stemp = s[j - 1];
252 if (ctemp != 1. || stemp != 0.) {
253 i__1 = *n;
254 for (i__ = 1; i__ <= i__1; ++i__) {
255 temp = a_ref(j, i__);
256 a_ref(j, i__) = ctemp * temp - stemp * a_ref(1,
257 i__);
258 a_ref(1, i__) = stemp * temp + ctemp * a_ref(1,
259 i__);
260/* L70: */
261 }
262 }
263/* L80: */
264 }
265 }
266 } else if (template_blas_lsame(pivot, "B")) {
267 if (template_blas_lsame(direct, "F")) {
268 i__1 = *m - 1;
269 for (j = 1; j <= i__1; ++j) {
270 ctemp = c__[j];
271 stemp = s[j];
272 if (ctemp != 1. || stemp != 0.) {
273 i__2 = *n;
274 for (i__ = 1; i__ <= i__2; ++i__) {
275 temp = a_ref(j, i__);
276 a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp *
277 temp;
278 a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp *
279 temp;
280/* L90: */
281 }
282 }
283/* L100: */
284 }
285 } else if (template_blas_lsame(direct, "B")) {
286 for (j = *m - 1; j >= 1; --j) {
287 ctemp = c__[j];
288 stemp = s[j];
289 if (ctemp != 1. || stemp != 0.) {
290 i__1 = *n;
291 for (i__ = 1; i__ <= i__1; ++i__) {
292 temp = a_ref(j, i__);
293 a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp *
294 temp;
295 a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp *
296 temp;
297/* L110: */
298 }
299 }
300/* L120: */
301 }
302 }
303 }
304 } else if (template_blas_lsame(side, "R")) {
305
306/* Form A * P' */
307
308 if (template_blas_lsame(pivot, "V")) {
309 if (template_blas_lsame(direct, "F")) {
310 i__1 = *n - 1;
311 for (j = 1; j <= i__1; ++j) {
312 ctemp = c__[j];
313 stemp = s[j];
314 if (ctemp != 1. || stemp != 0.) {
315 i__2 = *m;
316 for (i__ = 1; i__ <= i__2; ++i__) {
317 temp = a_ref(i__, j + 1);
318 a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
319 i__, j);
320 a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__,
321 j);
322/* L130: */
323 }
324 }
325/* L140: */
326 }
327 } else if (template_blas_lsame(direct, "B")) {
328 for (j = *n - 1; j >= 1; --j) {
329 ctemp = c__[j];
330 stemp = s[j];
331 if (ctemp != 1. || stemp != 0.) {
332 i__1 = *m;
333 for (i__ = 1; i__ <= i__1; ++i__) {
334 temp = a_ref(i__, j + 1);
335 a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
336 i__, j);
337 a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__,
338 j);
339/* L150: */
340 }
341 }
342/* L160: */
343 }
344 }
345 } else if (template_blas_lsame(pivot, "T")) {
346 if (template_blas_lsame(direct, "F")) {
347 i__1 = *n;
348 for (j = 2; j <= i__1; ++j) {
349 ctemp = c__[j - 1];
350 stemp = s[j - 1];
351 if (ctemp != 1. || stemp != 0.) {
352 i__2 = *m;
353 for (i__ = 1; i__ <= i__2; ++i__) {
354 temp = a_ref(i__, j);
355 a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__,
356 1);
357 a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__,
358 1);
359/* L170: */
360 }
361 }
362/* L180: */
363 }
364 } else if (template_blas_lsame(direct, "B")) {
365 for (j = *n; j >= 2; --j) {
366 ctemp = c__[j - 1];
367 stemp = s[j - 1];
368 if (ctemp != 1. || stemp != 0.) {
369 i__1 = *m;
370 for (i__ = 1; i__ <= i__1; ++i__) {
371 temp = a_ref(i__, j);
372 a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__,
373 1);
374 a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__,
375 1);
376/* L190: */
377 }
378 }
379/* L200: */
380 }
381 }
382 } else if (template_blas_lsame(pivot, "B")) {
383 if (template_blas_lsame(direct, "F")) {
384 i__1 = *n - 1;
385 for (j = 1; j <= i__1; ++j) {
386 ctemp = c__[j];
387 stemp = s[j];
388 if (ctemp != 1. || stemp != 0.) {
389 i__2 = *m;
390 for (i__ = 1; i__ <= i__2; ++i__) {
391 temp = a_ref(i__, j);
392 a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp *
393 temp;
394 a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp *
395 temp;
396/* L210: */
397 }
398 }
399/* L220: */
400 }
401 } else if (template_blas_lsame(direct, "B")) {
402 for (j = *n - 1; j >= 1; --j) {
403 ctemp = c__[j];
404 stemp = s[j];
405 if (ctemp != 1. || stemp != 0.) {
406 i__1 = *m;
407 for (i__ = 1; i__ <= i__1; ++i__) {
408 temp = a_ref(i__, j);
409 a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp *
410 temp;
411 a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp *
412 temp;
413/* L230: */
414 }
415 }
416/* L240: */
417 }
418 }
419 }
420 }
421
422 return 0;
423
424/* End of DLASR */
425
426} /* dlasr_ */
427
428#undef a_ref
429
430
431#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
#define a_ref(a_1, a_2)
int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m, const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer *lda)
Definition: template_lapack_lasr.h:42