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