src/kernel/sorting.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_sorting_c = "$Id$";
  30
  31 #include <limits.h>
  32 #include "sysstuff.h"
  33 #include "smalloc.h"
  34 #include "sorting.h"
  35
  36 /*****************************************************************************
  37  *                                                                           *
  38  *                   Block sorting on coordinates                            *
  39  *                                                                           *
  40  *****************************************************************************/
  41
  42 static rvec *make_xblock(t_block *block,rvec x[])
  43 {
  44   int i,j,k,nr,n;
  45   rvec *xblock;
  46
  47   nr=block->nr;
  48   snew(xblock,nr);
  49   for (i=0; i<nr; i++)
  50     {
  51       for (j=0; j<DIM; j++) xblock[i][j]=0.0;
  52       for (j=block->index[i]; j<(int)(block->index[i+1]); j++)
  53         for (k=0; k<DIM; k++) xblock[i][k]+=x[j][k];
  54       n=block->index[i+1]-block->index[i];
  55       for (k=0; k<DIM; k++) xblock[i][k]/=n;
  56     }
  57   return xblock;
  58 }
  59
  60 static rvec *xblock; /* just global to bcomp1, used in qsort */
  61
  62 static int bomp1(const void *p1,const void *p2)
  63 {
  64   int i,i1,i2;
  65
  66   i1=*(int *)p1;
  67   i2=*(int *)p2;
  68   for (i=0; i<DIM; i++)
  69     if (xblock[i1][i]<xblock[i2][i]) return -1;
  70     else if (xblock[i1][i]>xblock[i2][i]) return 1;
  71   return 0;
  72 }
  73
  74 void sort_xblock(t_block *block,rvec x[],int renum[])
  75 {
  76   int i,nr,*invnum;
  77
  78   nr=block->nr;
  79   snew(invnum,nr);
  80   xblock=make_xblock(block,x);
  81   for (i=0; i<nr; i++) invnum[i]=i;
  82   qsort((void *)invnum,nr,(size_t)sizeof(invnum[0]),bomp1);
  83   for (i=0; i<nr; i++) renum[invnum[i]]=i;
  84   sfree(xblock);
  85   sfree(invnum);
  86 }
  87
  88 /*****************************************************************************
  89  *                                                                           *
  90  *                        Bond list sorting                                  *
  91  *                                                                           *
  92  *****************************************************************************/
  93
  94 static int bcomp2(const void *p1,const void *p2)
  95 {
  96   int done;
  97
  98   if ((((atom_id *)p1)[0])!=(((atom_id *)p2)[0]))
  99     done=((((atom_id *)p1)[0])-(((atom_id *)p2)[0]));
 100   else
 101     done=((((atom_id *)p1)[1])-(((atom_id *)p2)[1]));
 102 #ifdef DEBUG
 103   printf("bcomp2: [%d,%d] with [%d,%d] result %d\n",
 104           ((atom_id *)p1)[0],((atom_id *)p1)[1],
 105           ((atom_id *)p2)[0],((atom_id *)p2)[1],done);
 106 #endif
 107   return done;
 108 }
 109
 110 void sort_bond_list(t_bond bonds[],int nr)
 111 {
 112   qsort((void *)bonds,nr,(size_t)sizeof(bonds[0]),bcomp2);
 113 }