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