src/gromacs/gmxlib/gmx_sort.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2012, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source 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
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <stdlib.h>
  41 #include "gmx_sort.h"
  42
  43
  44 static void
  45 qsort_swapfunc(void *        a,
  46                void *        b,
  47                size_t        n,
  48                int           swaptype)
  49 {
  50     int  *  ia;
  51     int  *  ib;
  52     int     itmp;
  53
  54     char *  ca;
  55     char *  cb;
  56     char    ctmp;
  57
  58     if (swaptype <= 1)
  59     {
  60         ia = (int *)a;
  61         ib = (int *)b;
  62         for (; n > 0; ia += 1, ib += 1, n -= sizeof(int))
  63         {
  64             itmp      = *ia;
  65             *ia       = *ib;
  66             *ib       = itmp;
  67         }
  68     }
  69     else
  70     {
  71         ca = (char *)a;
  72         cb = (char *)b;
  73         for (; n > 0; ca += 1, cb += 1, n -= 1)
  74         {
  75             ctmp       = *ca;
  76             *ca        = *cb;
  77             *cb        = ctmp;
  78         }
  79     }
  80 }
  81
  82
  83 static void *
  84 qsort_med3(void *        a,
  85            void *        b,
  86            void *        c,
  87            int          (*compar) (const void *a, const void *b))
  88 {
  89     if (compar(a, b) < 0)
  90     {
  91         if (compar(b, c) < 0)
  92         {
  93             return b;
  94         }
  95         else if (compar(a, c) < 0)
  96         {
  97             return c;
  98         }
  99         else
 100         {
 101             return a;
 102         }
 103     }
 104     else
 105     {
 106         if (compar(b, c) > 0)
 107         {
 108             return b;
 109         }
 110         else if (compar(a, c) > 0)
 111         {
 112             return c;
 113         }
 114         else
 115         {
 116             return a;
 117         }
 118     }
 119 }
 120
 121
 122 void
 123 gmx_qsort(void *           base,
 124           size_t           nmemb,
 125           size_t           size,
 126           int            (*compar)(const void *, const void *))
 127 {
 128 #define QSORT_EXCH(a, b, t) (t = a, a = b, b = t)
 129 #define QSORT_SWAP(a, b) swaptype != 0 ? qsort_swapfunc(a, b, size, swaptype) : \
 130     (void)QSORT_EXCH(*(int *)(a), *(int *)(b), t)
 131
 132     char  *pa, *pb, *pc, *pd, *pl, *pm, *pn, *pv, *cbase;
 133     int    r, swaptype;
 134     int    t, v;
 135     size_t s, st;
 136
 137     cbase = (char *)base;
 138
 139     swaptype = (size_t)(cbase - (char *)0) % sizeof(int) || size % sizeof(int) ? 2 : size == sizeof(int) ? 0 : 1;
 140
 141     if (nmemb < 7)
 142     {
 143         /* Insertion sort on smallest arrays */
 144         for (pm = cbase + size; pm < cbase + nmemb*size; pm += size)
 145         {
 146             for (pl = pm; (pl > cbase) && compar((void *)(pl-size), (void *) pl) > 0; pl -= size)
 147             {
 148                 QSORT_SWAP(pl, pl-size);
 149             }
 150         }
 151         return;
 152     }
 153
 154     /* Small arrays, middle element */
 155     pm = cbase + (nmemb/2)*size;
 156
 157     if (nmemb > 7)
 158     {
 159         pl = cbase;
 160         pn = cbase + (nmemb-1)*size;
 161         if (nmemb > 40)
 162         {
 163             /* Big arrays, pseudomedian of 9 */
 164             s  = (nmemb/8)*size;
 165             pl = (char *)qsort_med3((void *)pl, (void *)((size_t)pl+s), (void *)((size_t)pl+2*s), compar);
 166             pm = (char *)qsort_med3((void *)((size_t)pm-s), (void *)pm, (void *)((size_t)pm+s), compar);
 167             pn = (char *)qsort_med3((void *)((size_t)pn-2*s), (void *)((size_t)pn-s), (void *)pn, compar);
 168         }
 169         /* Mid-size, med of 3 */
 170         pm = (char *)qsort_med3((void *)pl, (void *)pm, (void *)pn, compar);
 171     }
 172
 173     /* pv points to partition value */
 174     if (swaptype != 0)
 175     {
 176         pv = cbase;
 177         QSORT_SWAP(pv, pm);
 178     }
 179     else
 180     {
 181         pv = (char*)(void*)&v;
 182         v  = *(int *)pm;
 183     }
 184
 185     pa = pb = cbase;
 186     pc = pd = cbase + (nmemb-1)*size;
 187
 188     for (;; )
 189     {
 190         while (pb <= pc && (r = compar((void *)pb, (void *) pv)) <= 0)
 191         {
 192             if (r == 0)
 193             {
 194                 QSORT_SWAP(pa, pb);
 195                 pa += size;
 196             }
 197             pb += size;
 198         }
 199         while (pc >= pb && (r = compar((void *)pc, (void *) pv)) >= 0)
 200         {
 201             if (r == 0)
 202             {
 203                 QSORT_SWAP(pc, pd);
 204                 pd -= size;
 205             }
 206             pc -= size;
 207         }
 208         if (pb > pc)
 209         {
 210             break;
 211         }
 212         QSORT_SWAP(pb, pc);
 213         pb += size;
 214         pc -= size;
 215     }
 216     pn = cbase + nmemb*size;
 217
 218     s  = pa-cbase;
 219     st = pb-pa;
 220     if (st < s)
 221     {
 222         s = st;
 223     }
 224
 225     if (s > 0)
 226     {
 227         qsort_swapfunc(cbase, pb-s, s, swaptype);
 228     }
 229
 230     s  = pd-pc;
 231     st = pn-pd-size;
 232     if (st < s)
 233     {
 234         s = st;
 235     }
 236
 237     if (s > 0)
 238     {
 239         qsort_swapfunc(pb, pn-s, s, swaptype);
 240     }
 241
 242     if ((s = pb-pa) > size)
 243     {
 244         gmx_qsort(cbase, s/size, size, compar);
 245     }
 246
 247     if ((s = pd-pc) > size)
 248     {
 249         gmx_qsort(pn-s, s/size, size, compar);
 250     }
 251
 252 #undef QSORT_EXCH
 253 #undef QSORT_SWAP
 254
 255     return;
 256 }