src/gmxlib/gmx_lapack/ilasrt2.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, 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
  38 F77_FUNC(ilasrt2,ILASRT2)(const char *id,
  39          int *n,
  40          int *d__,
  41          int *key,
  42          int *info)
  43 {
  44     int i__1, i__2;
  45     int i__, j, d1, d2, d3, dir, tmp, endd;
  46     int stack[64], dmnmx, start;
  47     int tmpkey, stkpnt;
  48
  49     --key;
  50     --d__;
  51
  52     *info = 0;
  53     dir = -1;
  54     if (*id=='D' || *id=='d')
  55         dir = 0;
  56     else if (*id=='I' || *id=='i')
  57         dir = 1;
  58
  59     if (dir == -1) {
  60         *info = -1;
  61     } else if (*n < 0) {
  62         *info = -2;
  63     }
  64     if (*info != 0) {
  65         return;
  66     }
  67
  68     if (*n <= 1) {
  69         return;
  70     }
  71
  72     stkpnt = 1;
  73     stack[0] = 1;
  74     stack[1] = *n;
  75 L10:
  76     start = stack[(stkpnt << 1) - 2];
  77     endd = stack[(stkpnt << 1) - 1];
  78     --stkpnt;
  79     if (endd - start > 0) {
  80
  81         if (dir == 0) {
  82
  83             i__1 = endd;
  84             for (i__ = start + 1; i__ <= i__1; ++i__) {
  85                 i__2 = start + 1;
  86                 for (j = i__; j >= i__2; --j) {
  87                     if (d__[j] > d__[j - 1]) {
  88                         dmnmx = d__[j];
  89                         d__[j] = d__[j - 1];
  90                         d__[j - 1] = dmnmx;
  91                         tmpkey = key[j];
  92                         key[j] = key[j - 1];
  93                         key[j - 1] = tmpkey;
  94                     } else {
  95                         goto L30;
  96                     }
  97                 }
  98 L30:
  99                 ;
 100             }
 101
 102         } else {
 103
 104             i__1 = endd;
 105             for (i__ = start + 1; i__ <= i__1; ++i__) {
 106                 i__2 = start + 1;
 107                 for (j = i__; j >= i__2; --j) {
 108                     if (d__[j] < d__[j - 1]) {
 109                         dmnmx = d__[j];
 110                         d__[j] = d__[j - 1];
 111                         d__[j - 1] = dmnmx;
 112                         tmpkey = key[j];
 113                         key[j] = key[j - 1];
 114                         key[j - 1] = tmpkey;
 115                     } else {
 116                         goto L50;
 117                     }
 118                 }
 119 L50:
 120                 ;
 121             }
 122
 123         }
 124
 125     } else if (endd - start > 20) {
 126
 127         d1 = d__[start];
 128         d2 = d__[endd];
 129         i__ = (start + endd) / 2;
 130         d3 = d__[i__];
 131         if (d1 < d2) {
 132             if (d3 < d1) {
 133                 dmnmx = d1;
 134             } else if (d3 < d2) {
 135                 dmnmx = d3;
 136             } else {
 137                 dmnmx = d2;
 138             }
 139         } else {
 140             if (d3 < d2) {
 141                 dmnmx = d2;
 142             } else if (d3 < d1) {
 143                 dmnmx = d3;
 144             } else {
 145                 dmnmx = d1;
 146             }
 147         }
 148
 149         if (dir == 0) {
 150
 151             i__ = start - 1;
 152             j = endd + 1;
 153 L60:
 154 L70:
 155             --j;
 156             if (d__[j] < dmnmx) {
 157                 goto L70;
 158             }
 159 L80:
 160             ++i__;
 161             if (d__[i__] > dmnmx) {
 162                 goto L80;
 163             }
 164             if (i__ < j) {
 165                 tmp = d__[i__];
 166                 d__[i__] = d__[j];
 167                 d__[j] = tmp;
 168                 tmpkey = key[j];
 169                 key[j] = key[i__];
 170                 key[i__] = tmpkey;
 171                 goto L60;
 172             }
 173             if (j - start > endd - j - 1) {
 174                 ++stkpnt;
 175                 stack[(stkpnt << 1) - 2] = start;
 176                 stack[(stkpnt << 1) - 1] = j;
 177                 ++stkpnt;
 178                 stack[(stkpnt << 1) - 2] = j + 1;
 179                 stack[(stkpnt << 1) - 1] = endd;
 180             } else {
 181                 ++stkpnt;
 182                 stack[(stkpnt << 1) - 2] = j + 1;
 183                 stack[(stkpnt << 1) - 1] = endd;
 184                 ++stkpnt;
 185                 stack[(stkpnt << 1) - 2] = start;
 186                 stack[(stkpnt << 1) - 1] = j;
 187             }
 188         } else {
 189
 190             i__ = start - 1;
 191             j = endd + 1;
 192 L90:
 193 L100:
 194             --j;
 195             if (d__[j] > dmnmx) {
 196                 goto L100;
 197             }
 198 L110:
 199             ++i__;
 200             if (d__[i__] < dmnmx) {
 201                 goto L110;
 202             }
 203             if (i__ < j) {
 204                 tmp = d__[i__];
 205                 d__[i__] = d__[j];
 206                 d__[j] = tmp;
 207                 tmpkey = key[j];
 208                 key[j] = key[i__];
 209                 key[i__] = tmpkey;
 210                 goto L90;
 211             }
 212             if (j - start > endd - j - 1) {
 213                 ++stkpnt;
 214                 stack[(stkpnt << 1) - 2] = start;
 215                 stack[(stkpnt << 1) - 1] = j;
 216                 ++stkpnt;
 217                 stack[(stkpnt << 1) - 2] = j + 1;
 218                 stack[(stkpnt << 1) - 1] = endd;
 219             } else {
 220                 ++stkpnt;
 221                 stack[(stkpnt << 1) - 2] = j + 1;
 222                 stack[(stkpnt << 1) - 1] = endd;
 223                 ++stkpnt;
 224                 stack[(stkpnt << 1) - 2] = start;
 225                 stack[(stkpnt << 1) - 1] = j;
 226             }
 227         }
 228     }
 229     if (stkpnt > 0) {
 230         goto L10;
 231     }
 232
 233
 234     return;
 235 }