template_lapack_lasrt.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7