ergo
template_lapack_lasrt.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_LASRT_HEADER
36 #define TEMPLATE_LAPACK_LASRT_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *
41  info)
42 {
43 /* -- LAPACK routine (version 3.0) --
44  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
45  Courant Institute, Argonne National Lab, and Rice University
46  September 30, 1994
47 
48 
49  Purpose
50  =======
51 
52  Sort the numbers in D in increasing order (if ID = 'I') or
53  in decreasing order (if ID = 'D' ).
54 
55  Use Quick Sort, reverting to Insertion sort on arrays of
56  size <= 20. Dimension of STACK limits N to about 2**32.
57 
58  Arguments
59  =========
60 
61  ID (input) CHARACTER*1
62  = 'I': sort D in increasing order;
63  = 'D': sort D in decreasing order.
64 
65  N (input) INTEGER
66  The length of the array D.
67 
68  D (input/output) DOUBLE PRECISION array, dimension (N)
69  On entry, the array to be sorted.
70  On exit, D has been sorted into increasing order
71  (D(1) <= ... <= D(N) ) or into decreasing order
72  (D(1) >= ... >= D(N) ), depending on ID.
73 
74  INFO (output) INTEGER
75  = 0: successful exit
76  < 0: if INFO = -i, the i-th argument had an illegal value
77 
78  =====================================================================
79 
80 
81  Test the input paramters.
82 
83  Parameter adjustments */
84  /* System generated locals */
85  integer i__1, i__2;
86  /* Local variables */
87  integer endd, i__, j;
88  integer stack[64] /* was [2][32] */;
89  Treal dmnmx, d1, d2, d3;
90  integer start;
91  integer stkpnt, dir;
92  Treal tmp;
93 #define stack_ref(a_1,a_2) stack[(a_2)*2 + a_1 - 3]
94 
95  --d__;
96 
97  /* Function Body */
98  *info = 0;
99  dir = -1;
100  if (template_blas_lsame(id, "D")) {
101  dir = 0;
102  } else if (template_blas_lsame(id, "I")) {
103  dir = 1;
104  }
105  if (dir == -1) {
106  *info = -1;
107  } else if (*n < 0) {
108  *info = -2;
109  }
110  if (*info != 0) {
111  i__1 = -(*info);
112  template_blas_erbla("LASRT ", &i__1);
113  return 0;
114  }
115 
116 /* Quick return if possible */
117 
118  if (*n <= 1) {
119  return 0;
120  }
121 
122  stkpnt = 1;
123  stack_ref(1, 1) = 1;
124  stack_ref(2, 1) = *n;
125 L10:
126  start = stack_ref(1, stkpnt);
127  endd = stack_ref(2, stkpnt);
128  --stkpnt;
129  if (endd - start <= 20 && endd - start > 0) {
130 
131 /* Do Insertion sort on D( START:ENDD ) */
132 
133  if (dir == 0) {
134 
135 /* Sort into decreasing order */
136 
137  i__1 = endd;
138  for (i__ = start + 1; i__ <= i__1; ++i__) {
139  i__2 = start + 1;
140  for (j = i__; j >= i__2; --j) {
141  if (d__[j] > d__[j - 1]) {
142  dmnmx = d__[j];
143  d__[j] = d__[j - 1];
144  d__[j - 1] = dmnmx;
145  } else {
146  goto L30;
147  }
148 /* L20: */
149  }
150 L30:
151  ;
152  }
153 
154  } else {
155 
156 /* Sort into increasing order */
157 
158  i__1 = endd;
159  for (i__ = start + 1; i__ <= i__1; ++i__) {
160  i__2 = start + 1;
161  for (j = i__; j >= i__2; --j) {
162  if (d__[j] < d__[j - 1]) {
163  dmnmx = d__[j];
164  d__[j] = d__[j - 1];
165  d__[j - 1] = dmnmx;
166  } else {
167  goto L50;
168  }
169 /* L40: */
170  }
171 L50:
172  ;
173  }
174 
175  }
176 
177  } else if (endd - start > 20) {
178 
179 /* Partition D( START:ENDD ) and stack parts, largest one first
180 
181  Choose partition entry as median of 3 */
182 
183  d1 = d__[start];
184  d2 = d__[endd];
185  i__ = (start + endd) / 2;
186  d3 = d__[i__];
187  if (d1 < d2) {
188  if (d3 < d1) {
189  dmnmx = d1;
190  } else if (d3 < d2) {
191  dmnmx = d3;
192  } else {
193  dmnmx = d2;
194  }
195  } else {
196  if (d3 < d2) {
197  dmnmx = d2;
198  } else if (d3 < d1) {
199  dmnmx = d3;
200  } else {
201  dmnmx = d1;
202  }
203  }
204 
205  if (dir == 0) {
206 
207 /* Sort into decreasing order */
208 
209  i__ = start - 1;
210  j = endd + 1;
211 L60:
212 L70:
213  --j;
214  if (d__[j] < dmnmx) {
215  goto L70;
216  }
217 L80:
218  ++i__;
219  if (d__[i__] > dmnmx) {
220  goto L80;
221  }
222  if (i__ < j) {
223  tmp = d__[i__];
224  d__[i__] = d__[j];
225  d__[j] = tmp;
226  goto L60;
227  }
228  if (j - start > endd - j - 1) {
229  ++stkpnt;
230  stack_ref(1, stkpnt) = start;
231  stack_ref(2, stkpnt) = j;
232  ++stkpnt;
233  stack_ref(1, stkpnt) = j + 1;
234  stack_ref(2, stkpnt) = endd;
235  } else {
236  ++stkpnt;
237  stack_ref(1, stkpnt) = j + 1;
238  stack_ref(2, stkpnt) = endd;
239  ++stkpnt;
240  stack_ref(1, stkpnt) = start;
241  stack_ref(2, stkpnt) = j;
242  }
243  } else {
244 
245 /* Sort into increasing order */
246 
247  i__ = start - 1;
248  j = endd + 1;
249 L90:
250 L100:
251  --j;
252  if (d__[j] > dmnmx) {
253  goto L100;
254  }
255 L110:
256  ++i__;
257  if (d__[i__] < dmnmx) {
258  goto L110;
259  }
260  if (i__ < j) {
261  tmp = d__[i__];
262  d__[i__] = d__[j];
263  d__[j] = tmp;
264  goto L90;
265  }
266  if (j - start > endd - j - 1) {
267  ++stkpnt;
268  stack_ref(1, stkpnt) = start;
269  stack_ref(2, stkpnt) = j;
270  ++stkpnt;
271  stack_ref(1, stkpnt) = j + 1;
272  stack_ref(2, stkpnt) = endd;
273  } else {
274  ++stkpnt;
275  stack_ref(1, stkpnt) = j + 1;
276  stack_ref(2, stkpnt) = endd;
277  ++stkpnt;
278  stack_ref(1, stkpnt) = start;
279  stack_ref(2, stkpnt) = j;
280  }
281  }
282  }
283  if (stkpnt > 0) {
284  goto L10;
285  }
286  return 0;
287 
288 /* End of DLASRT */
289 
290 } /* dlasrt_ */
291 
292 #undef stack_ref
293 
294 
295 #endif
int integer
Definition: template_blas_common.h:38
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44
#define stack_ref(a_1, a_2)