ergo
template_lapack_lasrt.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_LASRT_HEADER
38#define TEMPLATE_LAPACK_LASRT_HEADER
39
40
41template<class Treal>
42int 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;
127L10:
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 }
152L30:
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 }
173L50:
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;
213L60:
214L70:
215 --j;
216 if (d__[j] < dmnmx) {
217 goto L70;
218 }
219L80:
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;
251L90:
252L100:
253 --j;
254 if (d__[j] > dmnmx) {
255 goto L100;
256 }
257L110:
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
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 stack_ref(a_1, a_2)
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:42