ergo
template_lapack_lasr.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LASR_HEADER
36 #define TEMPLATE_LAPACK_LASR_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m,
41  const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer *
42  lda)
43 {
44 /* -- LAPACK auxiliary routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  October 31, 1992
48 
49 
50  Purpose
51  =======
52 
53  DLASR performs the transformation
54 
55  A := P*A, when SIDE = 'L' or 'l' ( Left-hand side )
56 
57  A := A*P', when SIDE = 'R' or 'r' ( Right-hand side )
58 
59  where A is an m by n real matrix and P is an orthogonal matrix,
60  consisting of a sequence of plane rotations determined by the
61  parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'
62  and z = n when SIDE = 'R' or 'r' ):
63 
64  When DIRECT = 'F' or 'f' ( Forward sequence ) then
65 
66  P = P( z - 1 )*...*P( 2 )*P( 1 ),
67 
68  and when DIRECT = 'B' or 'b' ( Backward sequence ) then
69 
70  P = P( 1 )*P( 2 )*...*P( z - 1 ),
71 
72  where P( k ) is a plane rotation matrix for the following planes:
73 
74  when PIVOT = 'V' or 'v' ( Variable pivot ),
75  the plane ( k, k + 1 )
76 
77  when PIVOT = 'T' or 't' ( Top pivot ),
78  the plane ( 1, k + 1 )
79 
80  when PIVOT = 'B' or 'b' ( Bottom pivot ),
81  the plane ( k, z )
82 
83  c( k ) and s( k ) must contain the cosine and sine that define the
84  matrix P( k ). The two by two plane rotation part of the matrix
85  P( k ), R( k ), is assumed to be of the form
86 
87  R( k ) = ( c( k ) s( k ) ).
88  ( -s( k ) c( k ) )
89 
90  This version vectorises across rows of the array A when SIDE = 'L'.
91 
92  Arguments
93  =========
94 
95  SIDE (input) CHARACTER*1
96  Specifies whether the plane rotation matrix P is applied to
97  A on the left or the right.
98  = 'L': Left, compute A := P*A
99  = 'R': Right, compute A:= A*P'
100 
101  DIRECT (input) CHARACTER*1
102  Specifies whether P is a forward or backward sequence of
103  plane rotations.
104  = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )
105  = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )
106 
107  PIVOT (input) CHARACTER*1
108  Specifies the plane for which P(k) is a plane rotation
109  matrix.
110  = 'V': Variable pivot, the plane (k,k+1)
111  = 'T': Top pivot, the plane (1,k+1)
112  = 'B': Bottom pivot, the plane (k,z)
113 
114  M (input) INTEGER
115  The number of rows of the matrix A. If m <= 1, an immediate
116  return is effected.
117 
118  N (input) INTEGER
119  The number of columns of the matrix A. If n <= 1, an
120  immediate return is effected.
121 
122  C, S (input) DOUBLE PRECISION arrays, dimension
123  (M-1) if SIDE = 'L'
124  (N-1) if SIDE = 'R'
125  c(k) and s(k) contain the cosine and sine that define the
126  matrix P(k). The two by two plane rotation part of the
127  matrix P(k), R(k), is assumed to be of the form
128  R( k ) = ( c( k ) s( k ) ).
129  ( -s( k ) c( k ) )
130 
131  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
132  The m by n matrix A. On exit, A is overwritten by P*A if
133  SIDE = 'R' or by A*P' if SIDE = 'L'.
134 
135  LDA (input) INTEGER
136  The leading dimension of the array A. LDA >= max(1,M).
137 
138  =====================================================================
139 
140 
141  Test the input parameters
142 
143  Parameter adjustments */
144  /* System generated locals */
145  integer a_dim1, a_offset, i__1, i__2;
146  /* Local variables */
147  integer info;
148  Treal temp;
149  integer i__, j;
150  Treal ctemp, stemp;
151 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
152 
153  --c__;
154  --s;
155  a_dim1 = *lda;
156  a_offset = 1 + a_dim1 * 1;
157  a -= a_offset;
158 
159  /* Function Body */
160  info = 0;
161  if (! (template_blas_lsame(side, "L") || template_blas_lsame(side, "R"))) {
162  info = 1;
163  } else if (! (template_blas_lsame(pivot, "V") || template_blas_lsame(pivot,
164  "T") || template_blas_lsame(pivot, "B"))) {
165  info = 2;
166  } else if (! (template_blas_lsame(direct, "F") || template_blas_lsame(direct,
167  "B"))) {
168  info = 3;
169  } else if (*m < 0) {
170  info = 4;
171  } else if (*n < 0) {
172  info = 5;
173  } else if (*lda < maxMACRO(1,*m)) {
174  info = 9;
175  }
176  if (info != 0) {
177  template_blas_erbla("LASR ", &info);
178  return 0;
179  }
180 
181 /* Quick return if possible */
182 
183  if (*m == 0 || *n == 0) {
184  return 0;
185  }
186  if (template_blas_lsame(side, "L")) {
187 
188 /* Form P * A */
189 
190  if (template_blas_lsame(pivot, "V")) {
191  if (template_blas_lsame(direct, "F")) {
192  i__1 = *m - 1;
193  for (j = 1; j <= i__1; ++j) {
194  ctemp = c__[j];
195  stemp = s[j];
196  if (ctemp != 1. || stemp != 0.) {
197  i__2 = *n;
198  for (i__ = 1; i__ <= i__2; ++i__) {
199  temp = a_ref(j + 1, i__);
200  a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
201  j, i__);
202  a_ref(j, i__) = stemp * temp + ctemp * a_ref(j,
203  i__);
204 /* L10: */
205  }
206  }
207 /* L20: */
208  }
209  } else if (template_blas_lsame(direct, "B")) {
210  for (j = *m - 1; j >= 1; --j) {
211  ctemp = c__[j];
212  stemp = s[j];
213  if (ctemp != 1. || stemp != 0.) {
214  i__1 = *n;
215  for (i__ = 1; i__ <= i__1; ++i__) {
216  temp = a_ref(j + 1, i__);
217  a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
218  j, i__);
219  a_ref(j, i__) = stemp * temp + ctemp * a_ref(j,
220  i__);
221 /* L30: */
222  }
223  }
224 /* L40: */
225  }
226  }
227  } else if (template_blas_lsame(pivot, "T")) {
228  if (template_blas_lsame(direct, "F")) {
229  i__1 = *m;
230  for (j = 2; j <= i__1; ++j) {
231  ctemp = c__[j - 1];
232  stemp = s[j - 1];
233  if (ctemp != 1. || stemp != 0.) {
234  i__2 = *n;
235  for (i__ = 1; i__ <= i__2; ++i__) {
236  temp = a_ref(j, i__);
237  a_ref(j, i__) = ctemp * temp - stemp * a_ref(1,
238  i__);
239  a_ref(1, i__) = stemp * temp + ctemp * a_ref(1,
240  i__);
241 /* L50: */
242  }
243  }
244 /* L60: */
245  }
246  } else if (template_blas_lsame(direct, "B")) {
247  for (j = *m; j >= 2; --j) {
248  ctemp = c__[j - 1];
249  stemp = s[j - 1];
250  if (ctemp != 1. || stemp != 0.) {
251  i__1 = *n;
252  for (i__ = 1; i__ <= i__1; ++i__) {
253  temp = a_ref(j, i__);
254  a_ref(j, i__) = ctemp * temp - stemp * a_ref(1,
255  i__);
256  a_ref(1, i__) = stemp * temp + ctemp * a_ref(1,
257  i__);
258 /* L70: */
259  }
260  }
261 /* L80: */
262  }
263  }
264  } else if (template_blas_lsame(pivot, "B")) {
265  if (template_blas_lsame(direct, "F")) {
266  i__1 = *m - 1;
267  for (j = 1; j <= i__1; ++j) {
268  ctemp = c__[j];
269  stemp = s[j];
270  if (ctemp != 1. || stemp != 0.) {
271  i__2 = *n;
272  for (i__ = 1; i__ <= i__2; ++i__) {
273  temp = a_ref(j, i__);
274  a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp *
275  temp;
276  a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp *
277  temp;
278 /* L90: */
279  }
280  }
281 /* L100: */
282  }
283  } else if (template_blas_lsame(direct, "B")) {
284  for (j = *m - 1; j >= 1; --j) {
285  ctemp = c__[j];
286  stemp = s[j];
287  if (ctemp != 1. || stemp != 0.) {
288  i__1 = *n;
289  for (i__ = 1; i__ <= i__1; ++i__) {
290  temp = a_ref(j, i__);
291  a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp *
292  temp;
293  a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp *
294  temp;
295 /* L110: */
296  }
297  }
298 /* L120: */
299  }
300  }
301  }
302  } else if (template_blas_lsame(side, "R")) {
303 
304 /* Form A * P' */
305 
306  if (template_blas_lsame(pivot, "V")) {
307  if (template_blas_lsame(direct, "F")) {
308  i__1 = *n - 1;
309  for (j = 1; j <= i__1; ++j) {
310  ctemp = c__[j];
311  stemp = s[j];
312  if (ctemp != 1. || stemp != 0.) {
313  i__2 = *m;
314  for (i__ = 1; i__ <= i__2; ++i__) {
315  temp = a_ref(i__, j + 1);
316  a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
317  i__, j);
318  a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__,
319  j);
320 /* L130: */
321  }
322  }
323 /* L140: */
324  }
325  } else if (template_blas_lsame(direct, "B")) {
326  for (j = *n - 1; j >= 1; --j) {
327  ctemp = c__[j];
328  stemp = s[j];
329  if (ctemp != 1. || stemp != 0.) {
330  i__1 = *m;
331  for (i__ = 1; i__ <= i__1; ++i__) {
332  temp = a_ref(i__, j + 1);
333  a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
334  i__, j);
335  a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__,
336  j);
337 /* L150: */
338  }
339  }
340 /* L160: */
341  }
342  }
343  } else if (template_blas_lsame(pivot, "T")) {
344  if (template_blas_lsame(direct, "F")) {
345  i__1 = *n;
346  for (j = 2; j <= i__1; ++j) {
347  ctemp = c__[j - 1];
348  stemp = s[j - 1];
349  if (ctemp != 1. || stemp != 0.) {
350  i__2 = *m;
351  for (i__ = 1; i__ <= i__2; ++i__) {
352  temp = a_ref(i__, j);
353  a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__,
354  1);
355  a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__,
356  1);
357 /* L170: */
358  }
359  }
360 /* L180: */
361  }
362  } else if (template_blas_lsame(direct, "B")) {
363  for (j = *n; j >= 2; --j) {
364  ctemp = c__[j - 1];
365  stemp = s[j - 1];
366  if (ctemp != 1. || stemp != 0.) {
367  i__1 = *m;
368  for (i__ = 1; i__ <= i__1; ++i__) {
369  temp = a_ref(i__, j);
370  a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__,
371  1);
372  a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__,
373  1);
374 /* L190: */
375  }
376  }
377 /* L200: */
378  }
379  }
380  } else if (template_blas_lsame(pivot, "B")) {
381  if (template_blas_lsame(direct, "F")) {
382  i__1 = *n - 1;
383  for (j = 1; j <= i__1; ++j) {
384  ctemp = c__[j];
385  stemp = s[j];
386  if (ctemp != 1. || stemp != 0.) {
387  i__2 = *m;
388  for (i__ = 1; i__ <= i__2; ++i__) {
389  temp = a_ref(i__, j);
390  a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp *
391  temp;
392  a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp *
393  temp;
394 /* L210: */
395  }
396  }
397 /* L220: */
398  }
399  } else if (template_blas_lsame(direct, "B")) {
400  for (j = *n - 1; j >= 1; --j) {
401  ctemp = c__[j];
402  stemp = s[j];
403  if (ctemp != 1. || stemp != 0.) {
404  i__1 = *m;
405  for (i__ = 1; i__ <= i__1; ++i__) {
406  temp = a_ref(i__, j);
407  a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp *
408  temp;
409  a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp *
410  temp;
411 /* L230: */
412  }
413  }
414 /* L240: */
415  }
416  }
417  }
418  }
419 
420  return 0;
421 
422 /* End of DLASR */
423 
424 } /* dlasr_ */
425 
426 #undef a_ref
427 
428 
429 #endif