mt_gomory.c

Go to the documentation of this file.

/* MTgomory "multi tableau gomory cut" provides an implementation for gomory
 * cuts derived from multiple tableau rows in the spirit of the work of
 * Andersen et al (IPCO 2007), Cornuejols (es presented in George Nemhauser
 * Birthday Conference in Atlanta 2007) and Gomory (presented in the same
 * conference).
 *
 * Copyright (C) 2007 Daniel Espinoza.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but 
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA 
 * */
/* ========================================================================= */
#include <math.h>
#include "EGlib.h"
#include "mt_gomory.h"
#include "mt_tableau.h"
/** @file 
 * @ingroup MTgomory */
/** @addtogroup MTgomory */
/** @{ */
/* ========================================================================= */
int MTgomoryCut(const int nrows, 
                const int ncols,
                const double*const rowval,
                const int*const rowind,
                const int*const rowbeg,
                const double*const f,
                const double*const base,
                double*const cutval,
                double*const work)
{
  int rval = 0,redo,pos,sense;
  register int i,j;
  EGgcIt_t gc;
  double ratio = 0;
  #if MT_DEBUG_LVL +5 <= DEBUG
  FILE* fout = EGsfopen("cur_tableau.txt","w+");
  rval = MTwriteRTableau(fout, nrows, ncols, rowval, rowind, rowbeg, f);
  CHECKRVALG(rval,CLEANUP);
  fclose(fout);
  #endif
  memset(&gc,0,sizeof(EGgcIt_t));
  /* check for basic assumptions in the function */
  if(MT_DEBUG_LVL <= DEBUG)
  {
    for(i = nrows ; i-- ; )
    {
      TESTG((rval=(fabs(base[i])<MT_ZERO_TOL)), CLEANUP, 
            "coefficient %d in base is %.12lf<%.12lf, and is considered"
            " as zero!", i, base[i], MT_ZERO_TOL);
      TESTG((rval=(f[i]<=(-0x1p-1+MT_ZERO_TOL))), CLEANUP, 
            "%d-th fractional value is outside bounds %lf <= %lf", 
            i, f[i], (-0x1p-1+MT_ZERO_TOL));
      TESTG((rval=(f[i]>=(0x1p-1-MT_ZERO_TOL))), CLEANUP, 
            "%d-th fractional value is outside bounds %lf >= %lf", 
            i, f[i], (0x1p-1-MT_ZERO_TOL));
    }
  }
  /* allocate memory */
  EGgcItInit(&gc,nrows);
  /* compute work[j] = col[j]*c and initialize cut coefficient */
  for( j = ncols ; j-- ; )
  {
    cutval[j]=0.0;
    work[j]=0.0;
  }
  work[ncols] = 0;
  for( j = nrows ; j-- ; )
  {
    /* compute in work[ncols] = 1/2 *base -f* \sigma * base */
    MESSAGE(MT_VERB_LVL+5,"j=%d, work[ncols]=%.6lf", j, work[ncols]);
    sense = ((EGgcItGetTuple(&gc))[j]);
    if(sense) work[ncols] += (0.5-f[j])*base[j];
    else work[ncols] += (0.5+f[j])*base[j];
    /* compute in work [i] col_i * \sigma * base */
    for( i = rowbeg[j] ; i < rowbeg[j+1] ;i++)
    {
      if(MT_VERB_LVL) ratio = work[rowind[i]];
      if(sense) work[rowind[i]] += base[j]*rowval[i];
      else work[rowind[i]] -= base[j]*rowval[i];
      MESSAGE(MT_VERB_LVL+5, "updating for row %d work[%d] from %lf to "
              "%lf", j, rowind[i], ratio, work[rowind[i]]);
    }
  }
  /* now we start computing the coefficients */
  do{
    if(MT_VERB_LVL+7 <= DEBUG)
    {
      const int*const tuple = EGgcItGetTuple(&gc);
      const int sz = EGgcItGetSize(&gc);
      fprintf(stderr,"Sign Patter: ");
      for( i = 0 ; i < sz ; i++)
      {
        fprintf(stderr,"%s",tuple[i]?"+":"-");
      }
      fprintf(stderr," work[ncols]=%.6lf\n",work[ncols]);
    }
    TESTG((rval=(work[ncols]<0)), CLEANUP, "Imposible!, n/2-base*sigma*f"
          " (%.6lf) < 0",work[ncols]);
    for( i = ncols ; i-- ; )
    {
      ratio = work[i]/work[ncols];
      if((work[i]>MT_ZERO_TOL) && (cutval[i]<ratio))
      {
        cutval[i] = ratio;
        MESSAGE(MT_VERB_LVL+7,"cutval[%d]=%lf, %lf/%lf=%lf", i, cutval[i], 
                work[i], work[ncols], ratio);
      }
    }
    redo = EGgcItNext(&gc);
    if(!redo) break;
    /* now we should re-compute work[*] */
    pos = EGgcItGetChange(&gc);
    sense = ((EGgcItGetTuple(&gc))[pos]);
    if(sense)
    {
      work[ncols] -= 2*base[pos]*f[pos];
      for( i = rowbeg[pos] ; i < rowbeg[pos+1] ;i++)
        work[rowind[i]] += 2*base[pos]*rowval[i];
    }
    else
    {
      work[ncols] += 2*base[pos]*f[pos];
      for( i = rowbeg[pos] ; i < rowbeg[pos+1] ;i++)
        work[rowind[i]] -= 2*base[pos]*rowval[i];
    }
  }while(1);
  CLEANUP:
  EGgcItClear(&gc);
  return rval;
}

/* ========================================================================= */
/** @} */
/* end of mt_gomory.c */

---