src/gmxlib/gmx_lapack/slasrt2.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmx_lapack.h"
  36
  37 void F77_FUNC(slasrt2,SLASRT2)(const char *id,
  38               int *n,
  39               float *d__,
  40               int * key,
  41               int *info)
  42 {
  43     int i__1, i__2;
  44
  45     int i__, j;
  46     float d1, d2, d3;
  47     int dir;
  48     float tmp;
  49     int endd;
  50     int stack[64];
  51     float dmnmx;
  52     int start;
  53     int tmpkey, stkpnt;
  54
  55     --key;
  56     --d__;
  57
  58     *info = 0;
  59     dir = -1;
  60     if (*id=='D' || *id=='d')
  61         dir = 0;
  62     else if (*id=='I' || *id=='i')
  63         dir = 1;
  64
  65     if (dir == -1) {
  66         *info = -1;
  67     } else if (*n < 0) {
  68         *info = -2;
  69     }
  70     if (*info != 0) {
  71         i__1 = -(*info);
  72         return;
  73     }
  74
  75     if (*n <= 1) {
  76         return;
  77     }
  78
  79     stkpnt = 1;
  80     stack[0] = 1;
  81     stack[1] = *n;
  82 L10:
  83     start = stack[(stkpnt << 1) - 2];
  84     endd = stack[(stkpnt << 1) - 1];
  85     --stkpnt;
  86     if (endd - start > 0) {
  87
  88         if (dir == 0) {
  89
  90             i__1 = endd;
  91             for (i__ = start + 1; i__ <= i__1; ++i__) {
  92                 i__2 = start + 1;
  93                 for (j = i__; j >= i__2; --j) {
  94                     if (d__[j] > d__[j - 1]) {
  95                         dmnmx = d__[j];
  96                         d__[j] = d__[j - 1];
  97                         d__[j - 1] = dmnmx;
  98                         tmpkey = key[j];
  99                         key[j] = key[j - 1];
 100                         key[j - 1] = tmpkey;
 101                     } else {
 102                         break;
 103                     }
 104                 }
 105             }
 106
 107         } else {
 108
 109             i__1 = endd;
 110             for (i__ = start + 1; i__ <= i__1; ++i__) {
 111                 i__2 = start + 1;
 112                 for (j = i__; j >= i__2; --j) {
 113                     if (d__[j] < d__[j - 1]) {
 114                         dmnmx = d__[j];
 115                         d__[j] = d__[j - 1];
 116                         d__[j - 1] = dmnmx;
 117                         tmpkey = key[j];
 118                         key[j] = key[j - 1];
 119                         key[j - 1] = tmpkey;
 120                     } else {
 121                         break;
 122                     }
 123                 }
 124             }
 125
 126         }
 127
 128     } else if (endd - start > 20) {
 129
 130         d1 = d__[start];
 131         d2 = d__[endd];
 132         i__ = (start + endd) / 2;
 133         d3 = d__[i__];
 134         if (d1 < d2) {
 135             if (d3 < d1) {
 136                 dmnmx = d1;
 137             } else if (d3 < d2) {
 138                 dmnmx = d3;
 139             } else {
 140                 dmnmx = d2;
 141             }
 142         } else {
 143             if (d3 < d2) {
 144                 dmnmx = d2;
 145             } else if (d3 < d1) {
 146                 dmnmx = d3;
 147             } else {
 148                 dmnmx = d1;
 149             }
 150         }
 151
 152         if (dir == 0) {
 153
 154             i__ = start - 1;
 155             j = endd + 1;
 156 L60:
 157 L70:
 158             --j;
 159             if (d__[j] < dmnmx) {
 160                 goto L70;
 161             }
 162 L80:
 163             ++i__;
 164             if (d__[i__] > dmnmx) {
 165                 goto L80;
 166             }
 167             if (i__ < j) {
 168                 tmp = d__[i__];
 169                 d__[i__] = d__[j];
 170                 d__[j] = tmp;
 171                 tmpkey = key[j];
 172                 key[j] = key[i__];
 173                 key[i__] = tmpkey;
 174                 goto L60;
 175             }
 176             if (j - start > endd - j - 1) {
 177                 ++stkpnt;
 178                 stack[(stkpnt << 1) - 2] = start;
 179                 stack[(stkpnt << 1) - 1] = j;
 180                 ++stkpnt;
 181                 stack[(stkpnt << 1) - 2] = j + 1;
 182                 stack[(stkpnt << 1) - 1] = endd;
 183             } else {
 184                 ++stkpnt;
 185                 stack[(stkpnt << 1) - 2] = j + 1;
 186                 stack[(stkpnt << 1) - 1] = endd;
 187                 ++stkpnt;
 188                 stack[(stkpnt << 1) - 2] = start;
 189                 stack[(stkpnt << 1) - 1] = j;
 190             }
 191         } else {
 192
 193             i__ = start - 1;
 194             j = endd + 1;
 195 L90:
 196 L100:
 197             --j;
 198             if (d__[j] > dmnmx) {
 199                 goto L100;
 200             }
 201 L110:
 202             ++i__;
 203             if (d__[i__] < dmnmx) {
 204                 goto L110;
 205             }
 206             if (i__ < j) {
 207                 tmp = d__[i__];
 208                 d__[i__] = d__[j];
 209                 d__[j] = tmp;
 210                 tmpkey = key[j];
 211                 key[j] = key[i__];
 212                 key[i__] = tmpkey;
 213                 goto L90;
 214             }
 215             if (j - start > endd - j - 1) {
 216                 ++stkpnt;
 217                 stack[(stkpnt << 1) - 2] = start;
 218                 stack[(stkpnt << 1) - 1] = j;
 219                 ++stkpnt;
 220                 stack[(stkpnt << 1) - 2] = j + 1;
 221                 stack[(stkpnt << 1) - 1] = endd;
 222             } else {
 223                 ++stkpnt;
 224                 stack[(stkpnt << 1) - 2] = j + 1;
 225                 stack[(stkpnt << 1) - 1] = endd;
 226                 ++stkpnt;
 227                 stack[(stkpnt << 1) - 2] = start;
 228                 stack[(stkpnt << 1) - 1] = j;
 229             }
 230         }
 231     }
 232     if (stkpnt > 0) {
 233         goto L10;
 234     }
 235
 236     return;
 237 }