ergo
template_lapack_lasr.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_LASR_HEADER
38 #define TEMPLATE_LAPACK_LASR_HEADER
39 
40 
41 template<class Treal>
42 int 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
mat::side
side
Definition: Matrix.h:75
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
a_ref
#define a_ref(a_1, a_2)
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
template_lapack_lasr
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