Merge changes If0c2f9d9,Ic2a8ad82 into release-4-5-patches

* changes:
  Fixed g_sham using more than three dimensions
  Integer multiplication that checks for overflow

diff --git a/include/maths.h b/include/maths.h
index ed7503205c0f96b1a07cb5fd53e43744d417dafe..c4eda03add0b05b48a7fb9774509f2d0e4b1aa71 100644 (file)
--- a/include/maths.h
+++ b/include/maths.h
@@ -165,6 +165,14 @@ gmx_log2(real x)
     return log( x ) * iclog2;
 }
 
+/*! /brief Multiply two large ints
+ *
+ *  Returns true when overflow did not occur.
+ */
+gmx_bool
+check_int_multiply_for_overflow(gmx_large_int_t a,
+                                gmx_large_int_t b,
+                                gmx_large_int_t *result);
 
 #ifdef __cplusplus
 }

diff --git a/include/types/simple.h b/include/types/simple.h
index 3146f9b73125b5d88e8d48f7849abc1da4e7c761..8279e135de40bb9c91bd7bef93d9318e63c67a55 100644 (file)
--- a/include/types/simple.h
+++ b/include/types/simple.h
@@ -174,6 +174,7 @@ typedef long long int gmx_large_int_t;
 #define gmx_large_int_pfmt "%lld"
 #define SIZEOF_GMX_LARGE_INT  8
 #define GMX_LARGE_INT_MAX     9223372036854775807LL
+#define GMX_LARGE_INT_MIN     (-GMX_LARGE_INT_MAX - 1LL)
 #define GMX_MPI_LARGE_INT MPI_LONG_LONG_INT
 
 #elif ( (defined LONG_MAX && LONG_MAX==9223372036854775807L) || (defined SIZEOF_LONG_INT && SIZEOF_LONG_INT==8) )
@@ -184,6 +185,7 @@ typedef long int gmx_large_int_t;
 #define gmx_large_int_pfmt "%ld"
 #define SIZEOF_GMX_LARGE_INT  8
 #define GMX_LARGE_INT_MAX     9223372036854775807LL
+#define GMX_LARGE_INT_MIN     (-GMX_LARGE_INT_MAX - 1LL)
 #define GMX_MPI_LARGE_INT MPI_LONG_INT
 
 #elif ( (defined INT_MAX && INT_MAX==9223372036854775807L) || (defined SIZEOF_INT && SIZEOF_INT==8) )
@@ -194,6 +196,7 @@ typedef int gmx_large_int_t;
 #define gmx_large_int_pfmt  "%d"
 #define SIZEOF_GMX_LARGE_INT  8
 #define GMX_LARGE_INT_MAX     9223372036854775807LL
+#define GMX_LARGE_INT_MIN     (-GMX_LARGE_INT_MAX - 1LL)
 #define GMX_MPI_LARGE_INT MPI_INT
 
 #elif ( (defined INT_MAX && INT_MAX==2147483647) || (defined SIZEOF_INT && SIZEOF_INT==4) )
@@ -204,6 +207,7 @@ typedef int gmx_large_int_t;
 #define gmx_large_int_pfmt "%d"
 #define SIZEOF_GMX_LARGE_INT  4
 #define GMX_LARGE_INT_MAX     2147483647
+#define GMX_LARGE_INT_MIN     (-GMX_LARGE_INT_MAX - 1)
 #define GMX_MPI_LARGE_INT MPI_INT
 
 #else

diff --git a/src/gmxlib/maths.c b/src/gmxlib/maths.c
index cfa55f6953fdb281e809482ec348d5d207c64025..81e49aae4a86685c5854b2fd51eadb69f551c32b 100644 (file)
--- a/src/gmxlib/maths.c
+++ b/src/gmxlib/maths.c
@@ -678,3 +678,33 @@ gmx_bool gmx_isfinite(real x)
 #endif
     return returnval;
 }
+
+gmx_bool
+check_int_multiply_for_overflow(gmx_large_int_t a,
+                                gmx_large_int_t b,
+                                gmx_large_int_t *result)
+{
+    gmx_large_int_t sign = 1;
+    if((0 == a) || (0 == b))
+    {
+        *result = 0;
+        return TRUE;
+    }
+    if(a < 0)
+    {
+        a = -a;
+        sign = -sign;
+    }
+    if(b < 0)
+    {
+        b = -b;
+        sign = -sign;
+    }
+    if(GMX_LARGE_INT_MAX / b < a)
+    {
+        *result = (sign > 0) ? GMX_LARGE_INT_MAX : GMX_LARGE_INT_MIN;
+        return FALSE;
+    }
+    *result = sign * a * b;
+    return TRUE;
+}

diff --git a/src/tools/gmx_sham.c b/src/tools/gmx_sham.c
index 668cd47ef8440db80563fb04ea4e7f919b091692..6078f60e65d03dd910903bfa125cf28ecc27614c 100644 (file)
--- a/src/tools/gmx_sham.c
+++ b/src/tools/gmx_sham.c
@@ -67,9 +67,10 @@ static int index3(int *ibox,int x,int y,int z)
   return (ibox[2]*(ibox[1]*x+y)+z);
 }
 
-static int indexn(int ndim,int *ibox,int *nxyz) 
+static gmx_large_int_t indexn(int ndim, const int *ibox, const int *nxyz)
 {
-  int d,dd,k,kk;
+    gmx_large_int_t d, dd;
+    int k, kk;
   
   /* Compute index in 1-D array */
   d = 0;
@@ -180,7 +181,7 @@ static void normalize_p_e(int len,double *P,int *nbin,real *E,real pmin)
 }
 
 typedef struct {
-  int index;
+  gmx_large_int_t index;
   real ener;
 } t_minimum;
 
@@ -197,70 +198,183 @@ static int comp_minima(const void *a,const void *b)
     return 0;
 }
 
+static inline
+void print_minimum(FILE *fp, int num, const t_minimum *min)
+{
+    fprintf(fp,
+            "Minimum %d at index " gmx_large_int_pfmt " energy %10.3f\n",
+            num, min->index, min->ener);
+}
+
+static inline
+void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm)
+{
+    print_minimum(fp, num, min);
+    mm[num].index = min->index;
+    mm[num].ener  = min->ener;
+}
+
+static inline
+gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
+                                     int dimension_index,
+                                     int dimension_min,
+                                     int neighbour_index,
+                                     real *W)
+{
+    return ((dimension_index == dimension_min) ||
+            ((dimension_index > dimension_min) &&
+             (this_min->ener < W[neighbour_index])));
+    /* Note over/underflow within W cannot occur. */
+}
+
+static inline
+gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
+                                     int dimension_index,
+                                     int dimension_max,
+                                     int neighbour_index,
+                                     real *W)
+{
+    return ((dimension_index == dimension_max) ||
+            ((dimension_index < dimension_max) &&
+             (this_min->ener < W[neighbour_index])));
+    /* Note over/underflow within W cannot occur. */
+}
+
 static void pick_minima(const char *logfile,int *ibox,int ndim,int len,real W[])
 {
-  FILE *fp;
-  int  i,j,k,ijk,nmin;
-  gmx_bool bMin;
-  t_minimum *mm;
+    FILE *fp;
+    int  i,j,k,nmin;
+    t_minimum *mm, this_min;
+    int *this_point;
+    int loopmax, loopcounter;
   
-  snew(mm,len);
-  nmin = 0;
-  fp = ffopen(logfile,"w");
-  for(i=0; (i<ibox[0]); i++) {
-    for(j=0; (j<ibox[1]); j++) {
-      if (ndim == 3) {
-       for(k=0; (k<ibox[2]); k++) {
-         ijk    = index3(ibox,i,j,k);
-         bMin   = (((i == 0)       || ((i > 0)       && 
-                                       (W[ijk] < W[index3(ibox,i-1,j,k)]))) &&
-                   ((i == ibox[0]-1) || ((i < ibox[0]-1) && 
-                                       (W[ijk] < W[index3(ibox,i+1,j,k)]))) &&
-                   ((j == 0)       || ((j > 0)       && 
-                                       (W[ijk] < W[index3(ibox,i,j-1,k)]))) &&
-                   ((j == ibox[1]-1) || ((j < ibox[1]-1) && 
-                                       (W[ijk] < W[index3(ibox,i,j+1,k)]))) &&
-                   ((k == 0)       || ((k > 0)       && 
-                                       (W[ijk] < W[index3(ibox,i,j,k-1)]))) &&
-                   ((k == ibox[2]-1) || ((k < ibox[2]-1) && 
-                                       (W[ijk] < W[index3(ibox,i,j,k+1)]))));
-         if (bMin) {
-           fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
-                   nmin,ijk,W[ijk]);
-           mm[nmin].index = ijk;
-           mm[nmin].ener  = W[ijk];
-           nmin++;
-         }
-       }
-      }
-      else {
-       ijk    = index2(ibox,i,j);
-       bMin   = (((i == 0)       || ((i > 0)       && 
-                                     (W[ijk] < W[index2(ibox,i-1,j)]))) &&
-                 ((i == ibox[0]-1) || ((i < ibox[0]-1) && 
-                                     (W[ijk] < W[index2(ibox,i+1,j)]))) &&
-                 ((j == 0)       || ((j > 0)       && 
-                                     (W[ijk] < W[index2(ibox,i,j-1)]))) &&
-                 ((j == ibox[1]-1) || ((j < ibox[1]-1) && 
-                                     (W[ijk] < W[index2(ibox,i,j+1)]))));
-       if (bMin) {
-         fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
-                 nmin,ijk,W[ijk]);
-         mm[nmin].index = ijk;
-         mm[nmin].ener  = W[ijk];
-         nmin++;
-       }
-      }
+    snew(mm,len);
+    nmin = 0;
+    fp = ffopen(logfile,"w");
+    /* Loop over each element in the array of dimenion ndim seeking
+     * minima with respect to every dimension. Specialized loops for
+     * speed with ndim == 2 and ndim == 3. */
+    switch(ndim)
+    {
+    case 0:
+        /* This is probably impossible to reach anyway. */
+        break;
+    case 2:
+        for(i=0; (i<ibox[0]); i++) {
+            for(j=0; (j<ibox[1]); j++) {
+                /* Get the index of this point in the flat array */
+                this_min.index = index2(ibox,i,j);
+                this_min.ener = W[this_min.index];
+                if (is_local_minimum_from_below(&this_min, i, 0,         index2(ibox,i-1,j  ), W) &&
+                    is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox,i+1,j  ), W) &&
+                    is_local_minimum_from_below(&this_min, j, 0,         index2(ibox,i  ,j-1), W) &&
+                    is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox,i  ,j+1), W))
+                {
+                    add_minimum(fp, nmin, &this_min, mm);
+                    nmin++;
+                }
+            }
+        }
+        break;
+    case 3:
+        for(i=0; (i<ibox[0]); i++) {
+            for(j=0; (j<ibox[1]); j++) {
+                for(k=0; (k<ibox[2]); k++) {
+                    /* Get the index of this point in the flat array */
+                    this_min.index = index3(ibox,i,j,k);
+                    this_min.ener = W[this_min.index];
+                    if (is_local_minimum_from_below(&this_min, i, 0,         index3(ibox,i-1,j  ,k  ), W) &&
+                        is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox,i+1,j  ,k  ), W) &&
+                        is_local_minimum_from_below(&this_min, j, 0,         index3(ibox,i  ,j-1,k  ), W) &&
+                        is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox,i  ,j+1,k  ), W) &&
+                        is_local_minimum_from_below(&this_min, k, 0,         index3(ibox,i  ,j  ,k-1), W) &&
+                        is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox,i  ,j  ,k+1), W))
+                    {
+                        add_minimum(fp, nmin, &this_min, mm);
+                        nmin++;
+                    }
+                }
+            }
+        }
+        break;
+    default:
+        /* Note this treats ndim == 1 and ndim > 3 */
+
+        /* Set up an ndim-dimensional vector to loop over the points
+         * on the grid. (0,0,0, ... 0) is an acceptable place to
+         * start. */
+        snew(this_point, ndim);
+
+        /* Determine the number of points of the ndim-dimensional
+         * grid. */
+        loopmax = ibox[0];
+        for (i = 1; i < ndim; i++)
+        {
+            loopmax *= ibox[i];
+        }
+
+        loopcounter = 0;
+        while (loopmax > loopcounter)
+        {
+            gmx_bool bMin = TRUE;
+
+            /* Get the index of this_point in the flat array */
+            this_min.index = indexn(ndim, ibox, this_point);
+            this_min.ener = W[this_min.index];
+
+            /* Is this_point a minimum from above and below in each
+             * dimension? */
+            for (i = 0; bMin && (i < ndim); i++)
+            {
+                /* Save the index of this_point within the curent
+                 * dimension so we can change that index in the
+                 * this_point array for use with indexn(). */
+                int index = this_point[i];
+                this_point[i]--;
+                bMin = bMin &&
+                    is_local_minimum_from_below(&this_min, index, 0,         indexn(ndim, ibox, this_point), W);
+                this_point[i] += 2;
+                bMin = bMin &&
+                    is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
+                this_point[i]--;
+            }
+            if (bMin)
+            {
+                add_minimum(fp, nmin, &this_min, mm);
+                nmin++;
+            }
+
+            /* update global loop counter */
+            loopcounter++;
+
+            /* Avoid underflow of this_point[i] */
+            if (loopmax > loopcounter)
+            {
+                /* update this_point non-recursively */
+                i = ndim-1;
+                this_point[i]++;
+                while (ibox[i] == this_point[i])
+                {
+                    this_point[i] = 0;
+                    i--;
+                    /* this_point[i] cannot underflow because
+                     * loopmax > loopcounter. */
+                    this_point[i]++;
+                }
+            }
+        }
+
+        sfree(this_point);
+        break;
     }
-  }
-  qsort(mm,nmin,sizeof(mm[0]),comp_minima);
-  fprintf(fp,"Minima sorted after energy\n");
-  for(i=0; (i<nmin); i++) {
-    fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
-           i,mm[i].index,mm[i].ener);
-  }
-  ffclose(fp);
-  sfree(mm);
+    qsort(mm,nmin,sizeof(mm[0]),comp_minima);
+    fprintf(fp,"Minima sorted after energy\n");
+    for(i=0; (i<nmin); i++)
+    {
+        print_minimum(fp, i, &mm[i]);
+    }
+    ffclose(fp);
+    sfree(mm);
 }
 
 static void do_sham(const char *fn,const char *ndx,
@@ -778,6 +892,7 @@ int gmx_sham(int argc,char *argv[])
   double   *av,*sig,cum1,cum2,cum3,cum4,db;
   const char     *fn_ge,*fn_ene;
   output_env_t oenv;
+  gmx_large_int_t num_grid_points;
     
   t_filenm fnm[] = { 
     { efXVG, "-f",    "graph",    ffREAD   },
@@ -851,10 +966,10 @@ int gmx_sham(int argc,char *argv[])
   if (fn_ene && et_val)
     ehisto(opt2fn("-histo",NFILE,fnm),e_n,et_val,oenv);
 
-  snew(idim,3);
-  snew(ibox,3);
-  snew(rmin,3);
-  snew(rmax,3);
+  snew(idim,max(3,nset));
+  snew(ibox,max(3,nset));
+  snew(rmin,max(3,nset));
+  snew(rmax,max(3,nset));
   for(i=0; (i<min(3,nset)); i++) {
     idim[i] = nrdim[i];
     ibox[i] = nrbox[i];
@@ -868,6 +983,22 @@ int gmx_sham(int argc,char *argv[])
     rmax[i] = xmax[2];
   }
 
+  /* Check that the grid size is manageable. */
+  num_grid_points = ibox[0];
+  for(i = 1; i < nset; i++)
+  {
+      gmx_large_int_t result;
+      if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
+      {
+          gmx_fatal(FARGS,
+                    "The number of dimensions and grid points is too large for this tool\n"
+                    "to handle with what it knows about the architecture upon which it\n"
+                    "is running. Use a different machine or consult the GROMACS mailing list.");
+      }
+      num_grid_points = result;
+  }
+  /* The number of grid points fits in a gmx_large_int_t. */
+
   do_sham(opt2fn("-dist",NFILE,fnm),opt2fn("-bin",NFILE,fnm),
          opt2fn("-lp",NFILE,fnm),
          opt2fn("-ls",NFILE,fnm),opt2fn("-lsh",NFILE,fnm),