ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_LASRT_HEADER 00036 #define TEMPLATE_LAPACK_LASRT_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer * 00041 info) 00042 { 00043 /* -- LAPACK routine (version 3.0) -- 00044 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00045 Courant Institute, Argonne National Lab, and Rice University 00046 September 30, 1994 00047 00048 00049 Purpose 00050 ======= 00051 00052 Sort the numbers in D in increasing order (if ID = 'I') or 00053 in decreasing order (if ID = 'D' ). 00054 00055 Use Quick Sort, reverting to Insertion sort on arrays of 00056 size <= 20. Dimension of STACK limits N to about 2**32. 00057 00058 Arguments 00059 ========= 00060 00061 ID (input) CHARACTER*1 00062 = 'I': sort D in increasing order; 00063 = 'D': sort D in decreasing order. 00064 00065 N (input) INTEGER 00066 The length of the array D. 00067 00068 D (input/output) DOUBLE PRECISION array, dimension (N) 00069 On entry, the array to be sorted. 00070 On exit, D has been sorted into increasing order 00071 (D(1) <= ... <= D(N) ) or into decreasing order 00072 (D(1) >= ... >= D(N) ), depending on ID. 00073 00074 INFO (output) INTEGER 00075 = 0: successful exit 00076 < 0: if INFO = -i, the i-th argument had an illegal value 00077 00078 ===================================================================== 00079 00080 00081 Test the input paramters. 00082 00083 Parameter adjustments */ 00084 /* System generated locals */ 00085 integer i__1, i__2; 00086 /* Local variables */ 00087 integer endd, i__, j; 00088 integer stack[64] /* was [2][32] */; 00089 Treal dmnmx, d1, d2, d3; 00090 integer start; 00091 integer stkpnt, dir; 00092 Treal tmp; 00093 #define stack_ref(a_1,a_2) stack[(a_2)*2 + a_1 - 3] 00094 00095 --d__; 00096 00097 /* Function Body */ 00098 *info = 0; 00099 dir = -1; 00100 if (template_blas_lsame(id, "D")) { 00101 dir = 0; 00102 } else if (template_blas_lsame(id, "I")) { 00103 dir = 1; 00104 } 00105 if (dir == -1) { 00106 *info = -1; 00107 } else if (*n < 0) { 00108 *info = -2; 00109 } 00110 if (*info != 0) { 00111 i__1 = -(*info); 00112 template_blas_erbla("LASRT ", &i__1); 00113 return 0; 00114 } 00115 00116 /* Quick return if possible */ 00117 00118 if (*n <= 1) { 00119 return 0; 00120 } 00121 00122 stkpnt = 1; 00123 stack_ref(1, 1) = 1; 00124 stack_ref(2, 1) = *n; 00125 L10: 00126 start = stack_ref(1, stkpnt); 00127 endd = stack_ref(2, stkpnt); 00128 --stkpnt; 00129 if (endd - start <= 20 && endd - start > 0) { 00130 00131 /* Do Insertion sort on D( START:ENDD ) */ 00132 00133 if (dir == 0) { 00134 00135 /* Sort into decreasing order */ 00136 00137 i__1 = endd; 00138 for (i__ = start + 1; i__ <= i__1; ++i__) { 00139 i__2 = start + 1; 00140 for (j = i__; j >= i__2; --j) { 00141 if (d__[j] > d__[j - 1]) { 00142 dmnmx = d__[j]; 00143 d__[j] = d__[j - 1]; 00144 d__[j - 1] = dmnmx; 00145 } else { 00146 goto L30; 00147 } 00148 /* L20: */ 00149 } 00150 L30: 00151 ; 00152 } 00153 00154 } else { 00155 00156 /* Sort into increasing order */ 00157 00158 i__1 = endd; 00159 for (i__ = start + 1; i__ <= i__1; ++i__) { 00160 i__2 = start + 1; 00161 for (j = i__; j >= i__2; --j) { 00162 if (d__[j] < d__[j - 1]) { 00163 dmnmx = d__[j]; 00164 d__[j] = d__[j - 1]; 00165 d__[j - 1] = dmnmx; 00166 } else { 00167 goto L50; 00168 } 00169 /* L40: */ 00170 } 00171 L50: 00172 ; 00173 } 00174 00175 } 00176 00177 } else if (endd - start > 20) { 00178 00179 /* Partition D( START:ENDD ) and stack parts, largest one first 00180 00181 Choose partition entry as median of 3 */ 00182 00183 d1 = d__[start]; 00184 d2 = d__[endd]; 00185 i__ = (start + endd) / 2; 00186 d3 = d__[i__]; 00187 if (d1 < d2) { 00188 if (d3 < d1) { 00189 dmnmx = d1; 00190 } else if (d3 < d2) { 00191 dmnmx = d3; 00192 } else { 00193 dmnmx = d2; 00194 } 00195 } else { 00196 if (d3 < d2) { 00197 dmnmx = d2; 00198 } else if (d3 < d1) { 00199 dmnmx = d3; 00200 } else { 00201 dmnmx = d1; 00202 } 00203 } 00204 00205 if (dir == 0) { 00206 00207 /* Sort into decreasing order */ 00208 00209 i__ = start - 1; 00210 j = endd + 1; 00211 L60: 00212 L70: 00213 --j; 00214 if (d__[j] < dmnmx) { 00215 goto L70; 00216 } 00217 L80: 00218 ++i__; 00219 if (d__[i__] > dmnmx) { 00220 goto L80; 00221 } 00222 if (i__ < j) { 00223 tmp = d__[i__]; 00224 d__[i__] = d__[j]; 00225 d__[j] = tmp; 00226 goto L60; 00227 } 00228 if (j - start > endd - j - 1) { 00229 ++stkpnt; 00230 stack_ref(1, stkpnt) = start; 00231 stack_ref(2, stkpnt) = j; 00232 ++stkpnt; 00233 stack_ref(1, stkpnt) = j + 1; 00234 stack_ref(2, stkpnt) = endd; 00235 } else { 00236 ++stkpnt; 00237 stack_ref(1, stkpnt) = j + 1; 00238 stack_ref(2, stkpnt) = endd; 00239 ++stkpnt; 00240 stack_ref(1, stkpnt) = start; 00241 stack_ref(2, stkpnt) = j; 00242 } 00243 } else { 00244 00245 /* Sort into increasing order */ 00246 00247 i__ = start - 1; 00248 j = endd + 1; 00249 L90: 00250 L100: 00251 --j; 00252 if (d__[j] > dmnmx) { 00253 goto L100; 00254 } 00255 L110: 00256 ++i__; 00257 if (d__[i__] < dmnmx) { 00258 goto L110; 00259 } 00260 if (i__ < j) { 00261 tmp = d__[i__]; 00262 d__[i__] = d__[j]; 00263 d__[j] = tmp; 00264 goto L90; 00265 } 00266 if (j - start > endd - j - 1) { 00267 ++stkpnt; 00268 stack_ref(1, stkpnt) = start; 00269 stack_ref(2, stkpnt) = j; 00270 ++stkpnt; 00271 stack_ref(1, stkpnt) = j + 1; 00272 stack_ref(2, stkpnt) = endd; 00273 } else { 00274 ++stkpnt; 00275 stack_ref(1, stkpnt) = j + 1; 00276 stack_ref(2, stkpnt) = endd; 00277 ++stkpnt; 00278 stack_ref(1, stkpnt) = start; 00279 stack_ref(2, stkpnt) = j; 00280 } 00281 } 00282 } 00283 if (stkpnt > 0) { 00284 goto L10; 00285 } 00286 return 0; 00287 00288 /* End of DLASRT */ 00289 00290 } /* dlasrt_ */ 00291 00292 #undef stack_ref 00293 00294 00295 #endif