src/gromacs/utility/qsort_threadsafe.c

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