src/gromacs/utility/qsort_threadsafe.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2012,2014,2015, 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     if (size == 0)
 135     {
 136         return;
 137     }
 138
 139     cbase = (char *)base;
 140
 141     swaptype = (size_t)(cbase - (char *)0) % sizeof(int) || size % sizeof(int) ? 2 : size == sizeof(int) ? 0 : 1;
 142
 143     if (nmemb < 7)
 144     {
 145         /* Insertion sort on smallest arrays */
 146         for (pm = cbase + size; pm < cbase + nmemb*size; pm += size)
 147         {
 148             for (pl = pm; (pl > cbase) && compar((void *)(pl-size), (void *) pl) > 0; pl -= size)
 149             {
 150                 QSORT_SWAP(pl, pl-size);
 151             }
 152         }
 153         return;
 154     }
 155
 156     /* Small arrays, middle element */
 157     pm = cbase + (nmemb/2)*size;
 158
 159     if (nmemb > 7)
 160     {
 161         pl = cbase;
 162         pn = cbase + (nmemb-1)*size;
 163         if (nmemb > 40)
 164         {
 165             /* Big arrays, pseudomedian of 9 */
 166             s  = (nmemb/8)*size;
 167             pl = (char *)qsort_med3((void *)pl, (void *)((size_t)pl+s), (void *)((size_t)pl+2*s), compar);
 168             pm = (char *)qsort_med3((void *)((size_t)pm-s), (void *)pm, (void *)((size_t)pm+s), compar);
 169             pn = (char *)qsort_med3((void *)((size_t)pn-2*s), (void *)((size_t)pn-s), (void *)pn, compar);
 170         }
 171         /* Mid-size, med of 3 */
 172         pm = (char *)qsort_med3((void *)pl, (void *)pm, (void *)pn, compar);
 173     }
 174
 175     /* pv points to partition value */
 176     if (swaptype != 0)
 177     {
 178         pv = cbase;
 179         QSORT_SWAP(pv, pm);
 180     }
 181     else
 182     {
 183         v  = *(int *)pm;
 184         pv = (char*)(void*)&v;
 185     }
 186
 187     pa = pb = cbase;
 188     pc = pd = cbase + (nmemb-1)*size;
 189
 190     for (;; )
 191     {
 192         while (pb <= pc && (r = compar((void *)pb, (void *) pv)) <= 0)
 193         {
 194             if (r == 0)
 195             {
 196                 QSORT_SWAP(pa, pb);
 197                 pa += size;
 198             }
 199             pb += size;
 200         }
 201         while (pc >= pb && (r = compar((void *)pc, (void *) pv)) >= 0)
 202         {
 203             if (r == 0)
 204             {
 205                 QSORT_SWAP(pc, pd);
 206                 pd -= size;
 207             }
 208             pc -= size;
 209         }
 210         if (pb > pc)
 211         {
 212             break;
 213         }
 214         QSORT_SWAP(pb, pc);
 215         pb += size;
 216         pc -= size;
 217     }
 218     pn = cbase + nmemb*size;
 219
 220     s  = pa-cbase;
 221     st = pb-pa;
 222     if (st < s)
 223     {
 224         s = st;
 225     }
 226
 227     if (s > 0)
 228     {
 229         qsort_swapfunc(cbase, pb-s, s, swaptype);
 230     }
 231
 232     s  = pd-pc;
 233     st = pn-pd-size;
 234     if (st < s)
 235     {
 236         s = st;
 237     }
 238
 239     if (s > 0)
 240     {
 241         qsort_swapfunc(pb, pn-s, s, swaptype);
 242     }
 243
 244     if ((s = pb-pa) > size)
 245     {
 246         gmx_qsort(cbase, s/size, size, compar);
 247     }
 248
 249     if ((s = pd-pc) > size)
 250     {
 251         gmx_qsort(pn-s, s/size, size, compar);
 252     }
 253
 254 #undef QSORT_EXCH
 255 #undef QSORT_SWAP
 256
 257     return;
 258 }