src/gmxlib/gmx_lapack/dlasrt.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
  38 F77_FUNC(dlasrt,DLASRT)(const char *id,
  39         int *n,
  40         double *d__,
  41         int *info)
  42 {
  43     int i__1, i__2;
  44
  45     int i__, j;
  46     double d1, d2, d3;
  47     int dir;
  48     double tmp;
  49     int endd;
  50     int stack[64];
  51     double dmnmx;
  52     int start;
  53     int stkpnt;
  54
  55     --d__;
  56
  57     *info = 0;
  58     dir = -1;
  59     if (*id=='D' || *id=='d')
  60         dir = 0;
  61     else if (*id=='I' || *id=='i')
  62         dir = 1;
  63
  64     if (dir == -1) {
  65         *info = -1;
  66     } else if (*n < 0) {
  67         *info = -2;
  68     }
  69     if (*info != 0) {
  70         i__1 = -(*info);
  71         return;
  72     }
  73     if (*n <= 1) {
  74         return;
  75     }
  76
  77     stkpnt = 1;
  78     stack[0] = 1;
  79     stack[1] = *n;
  80 L10:
  81     start = stack[(stkpnt << 1) - 2];
  82     endd = stack[(stkpnt << 1) - 1];
  83     --stkpnt;
  84     if (endd - start <= 20 && endd - start > 0) {
  85
  86
  87         if (dir == 0) {
  88
  89             i__1 = endd;
  90             for (i__ = start + 1; i__ <= i__1; ++i__) {
  91                 i__2 = start + 1;
  92                 for (j = i__; j >= i__2; --j) {
  93                     if (d__[j] > d__[j - 1]) {
  94                         dmnmx = d__[j];
  95                         d__[j] = d__[j - 1];
  96                         d__[j - 1] = dmnmx;
  97                     } else {
  98                         goto L30;
  99                     }
 100                 }
 101 L30:
 102                 ;
 103             }
 104
 105         } else {
 106
 107             i__1 = endd;
 108             for (i__ = start + 1; i__ <= i__1; ++i__) {
 109                 i__2 = start + 1;
 110                 for (j = i__; j >= i__2; --j) {
 111                     if (d__[j] < d__[j - 1]) {
 112                         dmnmx = d__[j];
 113                         d__[j] = d__[j - 1];
 114                         d__[j - 1] = dmnmx;
 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                 goto L60;
 169             }
 170             if (j - start > endd - j - 1) {
 171                 ++stkpnt;
 172                 stack[(stkpnt << 1) - 2] = start;
 173                 stack[(stkpnt << 1) - 1] = j;
 174                 ++stkpnt;
 175                 stack[(stkpnt << 1) - 2] = j + 1;
 176                 stack[(stkpnt << 1) - 1] = endd;
 177             } else {
 178                 ++stkpnt;
 179                 stack[(stkpnt << 1) - 2] = j + 1;
 180                 stack[(stkpnt << 1) - 1] = endd;
 181                 ++stkpnt;
 182                 stack[(stkpnt << 1) - 2] = start;
 183                 stack[(stkpnt << 1) - 1] = j;
 184             }
 185         } else {
 186
 187             i__ = start - 1;
 188             j = endd + 1;
 189 L90:
 190 L100:
 191             --j;
 192             if (d__[j] > dmnmx) {
 193                 goto L100;
 194             }
 195 L110:
 196             ++i__;
 197             if (d__[i__] < dmnmx) {
 198                 goto L110;
 199             }
 200             if (i__ < j) {
 201                 tmp = d__[i__];
 202                 d__[i__] = d__[j];
 203                 d__[j] = tmp;
 204                 goto L90;
 205             }
 206             if (j - start > endd - j - 1) {
 207                 ++stkpnt;
 208                 stack[(stkpnt << 1) - 2] = start;
 209                 stack[(stkpnt << 1) - 1] = j;
 210                 ++stkpnt;
 211                 stack[(stkpnt << 1) - 2] = j + 1;
 212                 stack[(stkpnt << 1) - 1] = endd;
 213             } else {
 214                 ++stkpnt;
 215                 stack[(stkpnt << 1) - 2] = j + 1;
 216                 stack[(stkpnt << 1) - 1] = endd;
 217                 ++stkpnt;
 218                 stack[(stkpnt << 1) - 2] = start;
 219                 stack[(stkpnt << 1) - 1] = j;
 220             }
 221         }
 222     }
 223     if (stkpnt > 0) {
 224         goto L10;
 225     }
 226     return;
 227
 228 }