COding » MTgomory » src

mt_cplex_cbk.c

Go to the documentation of this file.
00001 /* MTgomory "multi tableau gomory cut" provides an implementation for gomory
00002  * cuts derived from multiple tableau rows in the spirit of the work of
00003  * Andersen et al (IPCO 2007), Cornuejols (es presented in George Nemhauser
00004  * Birthday Conference in Atlanta 2007) and Gomory (presented in the same
00005  * conference).
00006  *
00007  * Copyright (C) 2007 Daniel Espinoza.
00008  * 
00009  * This library is free software; you can redistribute it and/or modify it
00010  * under the terms of the GNU Lesser General Public License as published by the
00011  * Free Software Foundation; either version 2.1 of the License, or (at your
00012  * option) any later version.
00013  *
00014  * This library is distributed in the hope that it will be useful, but 
00015  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
00016  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
00017  * License for more details.
00018  *
00019  * You should have received a copy of the GNU Lesser General Public License
00020  * along with this library; if not, write to the Free Software Foundation,
00021  * Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA 
00022  * */
00023 /* ========================================================================= */
00024 #include <math.h>
00025 #include "cplex.h"
00026 #include "mt_cplex_cbk.h"
00027 #include "mt_gomory.h"
00028 #include "EGlib.h"
00029 /** @file 
00030  * @ingroup MTgomory */
00031 /** @addtogroup MTgomory */
00032 /** @{ */
00033 /* ========================================================================= */
00034 int MTcheckTabrow(
00035     const MTlp_t*const lp,
00036     const double*tableau,
00037     int*const status,
00038     double*const ratio,
00039     int*const etnz)
00040 {
00041   const int nrows = lp->nrows;
00042   const int ncols = lp->ncols;
00043   const char*const ctype = lp->ctype;
00044   const char*const rtype = lp->rtype;
00045   const int*const cstat = lp->cstat;
00046   const int*const rstat = lp->rstat;
00047   register int j;
00048   int nbasic, fail_tableau, tnz=0, rval=0;
00049   double absval, absmin=DBL_MAX, absmax=0, dtmp;
00050   *status = 0;
00051   /* process structural variables */
00052   for( j = ncols, nbasic = 0, fail_tableau = 0, tnz = 0 ; j-- ; )
00053   {
00054     dtmp = tableau[j];
00055     absval = EGabs(dtmp);
00056     if(absval < MT_ZERO_TOL) continue;
00057     if(absval > absmax) absmax = absval;
00058     if(absval < absmin) absmin = absval;
00059     if(cstat[j] == CPX_BASIC)
00060     {
00061       rval = (ctype[j] != CPX_BINARY && ctype[j] != CPX_INTEGER);
00062       #if MT_CCBK_USE_NAMES
00063       WARNINGL(MT_VERB_LVL+5, rval, "Imposible!, variable %s is basic,"
00064           " but type %c", lp->colnames[j], ctype[j]);
00065       #else
00066       WARNINGL(MT_VERB_LVL+5, rval, "Imposible!, variable %d is basic,"
00067           " but type %c", j, ctype[j]);
00068       #endif
00069       if(rval)
00070       {
00071         MTccbk_fail_tableau_btype++;
00072         fail_tableau = 1;
00073         rval = 0;
00074         WARNINGL(MT_VERB_LVL,1, "Found faulty tableau, |basic_coeff| = %.9lg, zero tol %lg ", absval, MT_ZERO_TOL);
00075         if(MT_VERB_LVL < DEBUG)
00076         {
00077           MTcplex_display_tableau(stderr,lp,tableau);
00078           fprintf(stderr,"\n");
00079         }
00080         break;
00081       }
00082       rval = ((absval < 1 - MT_ZERO_TOL) || (absval > 1+MT_ZERO_TOL));
00083       #if MT_CCBK_USE_NAMES
00084       WARNINGL(MT_VERB_LVL+5, rval, "Imposible!, basic variable %s"
00085           " has tableau value %.3lg/%.3lg", lp->colnames[j], 
00086           dtmp, tableau[j]);
00087       #else
00088       WARNINGL(MT_VERB_LVL+5, rval, "Imposible!, basic variable %d"
00089           " has tableau value %.3lg/%.3lg", j, dtmp, tableau[j]);
00090       #endif
00091       if(rval)
00092       {
00093         MTccbk_fail_tableau_coeff++;
00094         fail_tableau = 1;
00095         rval = 0;
00096         WARNINGL(MT_VERB_LVL,1, "Found faulty tableau, |basic_coeff| = %.9lg, zero tol %lg ", absval, MT_ZERO_TOL);
00097         if(MT_VERB_LVL < DEBUG)
00098         {
00099           MTcplex_display_tableau(stderr,lp,tableau);
00100           fprintf(stderr,"\n");
00101         }
00102         break;
00103       }
00104       nbasic++;
00105     }
00106     if(cstat[j] == CPX_BASIC) continue;
00107     /* integer variables only have fractional value in the tableau */
00108     if(ctype[j] == CPX_BINARY || ctype[j] == CPX_INTEGER)
00109     {
00110       dtmp -= floor(dtmp);
00111       if(dtmp >= 1 - MT_ZERO_TOL) dtmp = 0;
00112     }
00113     /* count tableau non-zeros */
00114     if(EGabs(dtmp) < MT_ZERO_TOL) continue;
00115     MESSAGE(MT_VERB_LVL+9,"Var %d in tableau",j);
00116     tnz++;
00117   }
00118   if(fail_tableau) goto CLEANUP;;
00119   if( nbasic != 1)
00120   {
00121     MTccbk_fail_tableau_nbasic++;
00122     fail_tableau = 1;
00123     WARNINGL(MT_VERB_LVL,1, "Found faulty tableau, nbasic = %d, expected 1", nbasic);
00124     if(MT_VERB_LVL < DEBUG)
00125     {
00126       MTcplex_display_tableau(stderr,lp,tableau);
00127       fprintf(stderr,"\n");
00128     }
00129   }
00130   if(fail_tableau) goto CLEANUP;
00131   /* process logical variables */
00132   tableau += ncols;
00133   for( j = nrows ; j-- ; )
00134   {
00135     dtmp = tableau[j];
00136     absval = EGabs(dtmp);
00137     if(absval < MT_ZERO_TOL) continue;
00138     if(absval > absmax) absmax = absval;
00139     if(absval < absmin) absmin = absval;
00140     if(rstat[j] == CPX_BASIC)
00141     {
00142       MTccbk_fail_tableau_nbasic++;
00143       nbasic++;
00144       WARNINGL(MT_VERB_LVL,1, "Found faulty tableau, nbasic = %d, expected 1", nbasic);
00145       if(MT_VERB_LVL<DEBUG)
00146       {
00147         MTcplex_display_tableau(stderr,lp,tableau-ncols);
00148         fprintf(stderr,"\n");
00149       }
00150       fail_tableau = 1;
00151       break;
00152     }
00153     if(rtype[j] == CPX_INTEGER || rtype[j] == CPX_BINARY)
00154     {
00155       dtmp -= floor(dtmp);
00156       if(dtmp >= 1- MT_ZERO_TOL) dtmp = 0;
00157     }
00158     if(EGabs(dtmp) < MT_ZERO_TOL) continue;
00159     MESSAGE(MT_VERB_LVL+9,"Logical %d in tableau",j);
00160     tnz++;
00161   }
00162   tableau -= ncols;
00163   if(fail_tableau) goto CLEANUP;
00164   /* now we check that the tableau at hand is not too bad numerically */
00165   if(absmin/absmax < MT_CCBK_MIN_TABLEAU_RATIO)
00166   {
00167     MTccbk_fail_tableau_ratio++; 
00168     WARNINGL(MT_VERB_LVL, 1, "Tableau row with ratio %.9lf < %.9lf", 
00169         absmin/absmax, MT_CCBK_MIN_TABLEAU_RATIO);
00170     fail_tableau = 1;
00171     goto CLEANUP;
00172   }
00173   CLEANUP:
00174   if(fail_tableau) *status = 1;
00175   if(ratio) *ratio = absmin/absmax;
00176   if(etnz) *etnz = tnz;
00177   return rval;
00178 }
00179 /* ========================================================================= */
00180 /** @brief simple handler to update f according to the complementation of
00181  * variables */
00182 #define UPDATE_F_AND_VAL(__f,__cstat,__bound,__val) do{\
00183   if((__cstat) == CPX_AT_LOWER){\
00184     (__f) -= (__val) * (__bound);}\
00185   else if((__cstat) == CPX_AT_UPPER){\
00186     (__f) -= (__val) * (__bound);\
00187     (__val) = -1*(__val);}\
00188   else{TESTG((rval=1),CLEANUP,"Imposible!, cstat is %d", (__cstat));}}while(0)  
00189 /* ========================================================================= */
00190 int MTcompressTabRow(
00191     const MTlp_t*const lp,
00192     const MTsol_t*const sol,
00193     double*rowval,
00194     int*rowind,
00195     double*const rhs,
00196     int*const nz)
00197 {
00198   const double*const x = sol->x;
00199   const double*const slack = sol->slack;
00200   const int nrows = lp->nrows;
00201   const int ncols = lp->ncols;
00202   const char*const ctype = lp->ctype;
00203   const char*const rtype = lp->rtype;
00204   const int*const cstat = lp->cstat;
00205   const int*const rstat = lp->rstat;
00206   register int j;
00207   double dtmp;  
00208   int tnz=0, rval=0;
00209   for( j = tnz = 0; j < ncols ; j++)
00210   {
00211     dtmp = rowval[j];
00212     if(cstat[j] == CPX_BASIC) continue;
00213     if(ctype[j] == CPX_BINARY || ctype[j] == CPX_INTEGER)
00214     {
00215       dtmp -= floor(dtmp);
00216       if(dtmp >= 1 - MT_ZERO_TOL) dtmp = 0;
00217     }
00218     if(EGabs(dtmp) < MT_ZERO_TOL) continue; 
00219     if( (ctype[j] == CPX_BINARY || ctype[j] == CPX_INTEGER) && dtmp >0.5)
00220       dtmp -= 1.0;
00221     UPDATE_F_AND_VAL(*rhs, cstat[j], x[j], dtmp);
00222     MESSAGE(MT_VERB_LVL+9,"Var %d in tableau",j);
00223     rowind[tnz] = j;
00224     rowval[tnz++] = dtmp;
00225   }
00226   for( j = 0 ; j < nrows ; j++)
00227   {
00228     dtmp = rowval[j+ncols];
00229     if(rtype[j] == CPX_BINARY || rtype[j] == CPX_INTEGER)
00230     {
00231       dtmp -= floor(dtmp);
00232       if(dtmp >= 1 - MT_ZERO_TOL) dtmp = 0;
00233     }
00234     if(EGabs(dtmp) < MT_ZERO_TOL) continue; 
00235     if( (rtype[j] == CPX_BINARY || rtype[j] == CPX_INTEGER) && dtmp >0.5)
00236       dtmp -= 1.0;
00237     UPDATE_F_AND_VAL(*rhs, rstat[j], slack[j], dtmp);
00238     MESSAGE(MT_VERB_LVL+9,"Logical %d in tableau",j);
00239     rowind[tnz] = j+ncols;
00240     rowval[tnz++] = dtmp;
00241   }
00242   if(MT_VERB_LVL+8 <= DEBUG)
00243   {
00244     fprintf(stderr,"Compress (and complemented) tableau:");
00245     MTcplex_display_compress_tableau(stderr, lp, rowval, rowind, tnz);
00246     fprintf(stderr," = %.3lf\n", *rhs); 
00247   }
00248   CLEANUP:
00249   *nz = tnz;
00250   return rval;
00251 }
00252 /* ========================================================================= */
00253 /** @brief resize an lp structure 
00254  * @param MTlp pointer to an #MTlp_t structure.
00255  * @param MTrsz new row-space (it will only grow).
00256  * @param MTcsz new col space (it will only grow).
00257  * @param MTtsz new non-zero space (it will only grow).
00258  * */
00259 static inline void MTlp_rsz(MTlp_t*const MTlp,
00260     const int rsz,
00261     const int csz,
00262     const int tsz)
00263 {
00264   if(MTlp->rspace < rsz)
00265   {
00266     MTlp->rspace = rsz;
00267     EGfree(MTlp->sense);
00268     EGfree(MTlp->rhs);
00269     EGfree(MTlp->matbeg);
00270     EGfree(MTlp->rtype);
00271     EGfree(MTlp->rstat);
00272     #if MT_CCBK_USE_NAMES
00273     EGfree(MTlp->rownames);
00274     MTlp->rownames = EGsMalloc(char*,rsz);
00275     #endif
00276     MTlp->sense = EGsMalloc(char,rsz);
00277     MTlp->rhs = EGsMalloc(double,rsz);
00278     MTlp->matbeg = EGsMalloc(int,rsz+1);
00279     MTlp->rtype = EGsMalloc(char,rsz);
00280     MTlp->rstat = EGsMalloc(int,rsz);
00281   }
00282   if(MTlp->cspace < csz)
00283   {
00284     MTlp->cspace = csz;
00285     EGfree(MTlp->ctype);
00286     EGfree(MTlp->cstat);
00287     EGfree(MTlp->ub);
00288     EGfree(MTlp->lb);
00289     #if MT_CCBK_USE_NAMES
00290     EGfree(MTlp->colnames);
00291     MTlp->colnames = EGsMalloc(char*,csz);
00292     #endif
00293     MTlp->ctype = EGsMalloc(char,csz);
00294     MTlp->cstat = EGsMalloc(int,csz);
00295     MTlp->lb = EGsMalloc(double,csz);
00296     MTlp->ub = EGsMalloc(double,csz);
00297   }
00298   if(MTlp->nzspace < tsz)
00299   {
00300     MTlp->nzspace = tsz;
00301     EGfree(MTlp->matind);
00302     EGfree(MTlp->matval);
00303     MTlp->matind = EGsMalloc(int,tsz);
00304     MTlp->matval = EGsMalloc(double,tsz);
00305   }
00306   return;
00307 }
00308 /* ========================================================================= */
00309 int MTcompute_integer_vars(const MTlp_t*const lp,
00310     const double*const x,
00311     int*const intvars,
00312     int*const nint,
00313     const double intgap)
00314 {
00315   const char*const ctype = lp->ctype;
00316   register int i = lp->ncols;
00317   int n=0;
00318   double dtmp=0;
00319   for( ; i-- ; )
00320   {
00321     /* discard non-integer variables */
00322     if(ctype[i] != CPX_BINARY && ctype[i] != CPX_INTEGER) continue;
00323     dtmp = fabs(x[i] - round(x[i]));
00324     if(dtmp <= intgap)
00325     {
00326       intvars[n++] = i;
00327       #if MT_CCBK_USE_NAMES
00328       MESSAGE(MT_VERB_LVL+9,"%s=%.9lf is integer-valued", 
00329           lp->colnames[i],x[i]);
00330       #else
00331       MESSAGE(MT_VERB_LVL+9,"X[%d]=%.9lf is integer-valued",i,x[i]);
00332       #endif
00333     }
00334   }
00335   *nint = n;
00336   return 0;
00337 }
00338 /* ========================================================================= */
00339 int MTcompute_fractional_vars(const MTlp_t*const lp,
00340     const double*const x,
00341     int*const fracvars,
00342     int*const nfrac,
00343     const double min_gap)
00344 {
00345   const char*const ctype = lp->ctype;
00346   register int i = lp->ncols;
00347   int n=0;
00348   double dtmp=0;
00349   for( ; i-- ; )
00350   {
00351     /* discard non-integer variables */
00352     if(ctype[i] != CPX_BINARY && ctype[i] != CPX_INTEGER) continue;
00353     dtmp = fabs(x[i] - round(x[i]));
00354     if(dtmp > min_gap )
00355     {
00356       fracvars[n++] = i;
00357       #if MT_CCBK_USE_NAMES
00358       MESSAGE(MT_VERB_LVL+9,"%s=%.9lf is integer-fractional", 
00359           lp->colnames[i],x[i]);
00360       #else
00361       MESSAGE(MT_VERB_LVL+9,"X[%d]=%.9lf is integer-fractional",i,x[i]);
00362       #endif
00363     }
00364   }
00365   *nfrac = n;
00366   return 0;
00367 }
00368 /* ========================================================================= */
00369 void MTdisplay_lp(const MTlp_t*const lp,
00370     const MTsol_t*const sol,
00371     FILE*stream)
00372 {
00373   const double*const x = sol->x;
00374   const double*const slack = sol->slack;
00375   register int i = lp->ncols;
00376   for( ; i-- ; )
00377   {
00378     #if MT_CCBK_USE_NAMES
00379     fprintf(stream,"%s[%lf,%lf] = %.6lf %d %c\n", lp->colnames[i],
00380         lp->lb[i], lp->ub[i], x[i], lp->cstat[i],
00381         lp->ctype[i]);
00382     #else
00383     fprintf(stream,"x_%d[%lf,%lf] = %.6lf %d %c\n", i, lp->lb[i],
00384         lp->ub[i], x[i], lp->cstat[i], 
00385         lp->ctype[i]);
00386     #endif
00387   }
00388   for( i = lp->nrows ; i-- ; )
00389   {
00390     #if MT_CCBK_USE_NAMES
00391     fprintf(stderr,"%s = %.6lf %d %c\n", lp->rownames[i],
00392         slack[i], lp->rstat[i], lp->rtype[i]);
00393     #else
00394     fprintf(stderr,"s_%d = %.6lf %d %c\n", i, slack[i], 
00395         lp->rstat[i], lp->rtype[i]);
00396     #endif
00397   }
00398 }
00399 /* ========================================================================= */
00400 char __mt_cplex_errbuf[4096];
00401 /* ========================================================================= */
00402 void MTcompute_f(const MTlp_t*const lp, 
00403     const double*const tableau,
00404     const double*const x,
00405     const double*const slack,
00406     double*const f)
00407 {
00408   const double*const stab = tableau + lp->ncols;
00409   double rhs = 0;
00410   register int i;
00411   for( i = lp->ncols ; i-- ; )
00412   {
00413     if(EGabs(tableau[i]) < MT_ZERO_TOL) continue;
00414     rhs += x[i]*tableau[i];
00415   }
00416   for( i = lp->nrows ; i-- ; )
00417   {
00418     if(EGabs(stab[i]) < MT_ZERO_TOL) continue;
00419     rhs += slack[i]*stab[i];
00420   }
00421   *f = rhs;
00422 }
00423 /* ========================================================================= */
00424 int MTcplex_binv_to_tableau(const MTlp_t*const lp, 
00425     const double*const binvrow, 
00426     double*const tableau)
00427 {
00428   int rval=0;
00429   const int nrows = lp->nrows;
00430   const int ncols = lp->ncols;
00431   const int nvars = nrows+ncols;
00432   register int i,j;
00433   for( i = nvars ; i-- ; ) tableau[i]=0;
00434   for( i = nrows ; i-- ; )
00435   {
00436     if(EGabs(binvrow[i])>MT_ZERO_TOL)
00437     {
00438       /* set-up structural values */
00439       for( j = lp->matbeg[i]; j < lp->matbeg[i+1] ; j++)
00440       {
00441         tableau[lp->matind[j]] += lp->matval[j]*binvrow[i];
00442       }
00443       /* set-up logical values */
00444       switch(lp->sense[i])
00445       {
00446         case 'L':
00447         case 'l':
00448         case 'R':
00449         case 'r':
00450           tableau[ncols+i] += binvrow[i];
00451           break;
00452         case 'G':
00453         case 'g':
00454           tableau[ncols+i] -= binvrow[i];
00455           break;
00456         case 'E':
00457         case 'e':
00458           break;
00459         default:
00460           TESTG((rval=1),CLEANUP,"Imposible!");
00461           break;
00462       }
00463     }
00464   }
00465   CLEANUP:
00466   return rval;
00467 }
00468 /* ========================================================================= */
00469 int MTlp_load_problem(CPXCENVptr env, 
00470     MTlp_t*const lp, 
00471     CPXLPptr CPXlp, 
00472     void*cbdata, 
00473     int wherefrom)
00474 {
00475   int rval = 0,i1,i2,isint;
00476   double dtmp;
00477   lp->ncols = CPXgetnumcols(env,CPXlp);
00478   lp->nrows = CPXgetnumrows(env,CPXlp);
00479   MTlp_rsz(lp,lp->nrows,lp->ncols,lp->nz);
00480   CPXgetrows( env, CPXlp, &rval, lp->matbeg, lp->matind, lp->matval, 
00481       0, &(lp->nz), 0, lp->nrows-1);
00482   lp->nz = -lp->nz;
00483   MTlp_rsz(lp,lp->nrows,lp->ncols,lp->nz);
00484   lp->matbeg[lp->nrows]=lp->nz;
00485   MESSAGE(MT_VERB_LVL+5, "NZ: %d COLS: %d ROWS: %d matbeg %p, matind %p,"
00486       " matval %p", lp->nz, lp->ncols, lp->nrows, lp->matbeg, 
00487       lp->matind, lp->matval);
00488   /* get rows */
00489   rval = CPXgetrows(env, CPXlp, &i1, lp->matbeg, lp->matind, 
00490       lp->matval, lp->nz, &i2, 0, lp->nrows-1);
00491   MTcplexCHECKRVALG(env,rval,CLEANUP);
00492   /* get sense and rhs */
00493   rval = CPXgetsense(env, CPXlp, lp->sense, 0 , lp->nrows-1);
00494   MTcplexCHECKRVALG(env,rval,CLEANUP);
00495   for( i1 = lp->nrows ; i1-- ; )
00496   {
00497     TESTG((rval=(lp->sense[i1] == 'R')), CLEANUP, 
00498         "We don't support problems with range constraints"); 
00499   }
00500   rval = CPXgetrhs(env, CPXlp, lp->rhs, 0, lp->nrows-1);
00501   MTcplexCHECKRVALG(env,rval,CLEANUP);
00502   /* get bounds */
00503   rval = CPXgetlb(env, CPXlp, lp->lb, 0, lp->ncols-1);
00504   MTcplexCHECKRVALG(env,rval,CLEANUP);
00505   rval = CPXgetub(env, CPXlp, lp->ub, 0, lp->ncols-1);
00506   MTcplexCHECKRVALG(env,rval,CLEANUP);
00507   /* get ctype and rtype */
00508   rval = CPXgetcallbackctype(env,cbdata,wherefrom,lp->ctype,0,lp->ncols-1);
00509   MTcplexCHECKRVALG(env,rval,CLEANUP);
00510   /* now we loop in each constraint to see if the asociated slack is integer or
00511    * not */
00512   for( i1 = lp->nrows ; i1-- ; )
00513   {
00514     isint=1;
00515     for( i2 = lp->matbeg[i1] ; i2 < lp->matbeg[i1+1] ; i2++)
00516     {
00517       dtmp = lp->matval[i2];
00518       if(EGabs(dtmp)<MT_ZERO_TOL) continue;
00519       switch(lp->ctype[lp->matind[i2]])
00520       {
00521         case CPX_BINARY:
00522         case CPX_INTEGER:
00523           break;
00524         case CPX_CONTINUOUS:
00525           isint=0;
00526           goto BREAK;
00527           break;
00528         default:
00529           TESTG( (rval=1), CLEANUP, "Unexpected var type %c for var %d",
00530               lp->ctype[lp->matind[i2]], lp->matind[i2]);
00531           break;      
00532       }
00533       dtmp -= floor(dtmp);
00534       if(dtmp>MT_ZERO_TOL)
00535       {
00536         isint=0;
00537         break;
00538       }
00539     }
00540   BREAK:
00541     lp->rtype[i1] = ((char)(isint ? CPX_INTEGER: CPX_CONTINUOUS));
00542   }
00543   /* we obtain the basis for the current lp */
00544   rval = CPXgetbase(env,CPXlp,lp->cstat,lp->rstat);
00545   MTcplexCHECKRVALG(env,rval,CLEANUP);
00546   /* if we use column/row names, we load them now */
00547   #if MT_CCBK_USE_NAMES
00548   CPXgetcolname(env, CPXlp, lp->colnames, lp->ccspace, 0, &(lp->ccsz),
00549       0, lp->ncols-1);
00550   EGfree(lp->ccspace);
00551   lp->ccsz = -lp->ccsz;
00552   lp->ccspace = EGsMalloc(char,lp->ccsz);
00553   rval = CPXgetcolname( env, CPXlp, lp->colnames, lp->ccspace, lp->ccsz, 
00554       &i1, 0, lp->ncols-1);
00555   MTcplexCHECKRVALG(env,rval,CLEANUP);
00556   CPXgetrowname(env, CPXlp, lp->rownames, lp->rcspace, 0, &(lp->rcsz),
00557       0, lp->nrows-1);
00558   EGfree(lp->rcspace);
00559   lp->rcsz = -lp->rcsz;
00560   lp->rcspace = EGsMalloc(char,lp->rcsz);
00561   rval = CPXgetrowname( env, CPXlp, lp->rownames, lp->rcspace, lp->rcsz, 
00562       &i1, 0, lp->nrows-1);
00563   MTcplexCHECKRVALG(env,rval,CLEANUP);
00564   #endif
00565   CLEANUP:
00566   return rval;
00567 }
00568 /* ========================================================================= */
00569 int MTccbk_fail_tableau_ratio = 0;
00570 int MTccbk_fail_tableau_coeff = 0;
00571 int MTccbk_fail_tableau_frac = 0;
00572 int MTccbk_fail_tableau_btype = 0;
00573 int MTccbk_fail_tableau_nbasic = 0;
00574 int MTccbk_fail_violation = 0;
00575 int MTccbk_fail_ratio = 0;
00576 double MTccbk_discarded_violation = 0;
00577 double MTccbk_discarded_ratio = DBL_MAX;
00578 /* ========================================================================= */
00579 int MTgomory_ccbk(CPXCENVptr env,
00580     void*cbdata,
00581     int wherefrom,
00582     void*cbhandle,
00583     int*useraction_p)
00584 {
00585   MTgomory_ccbk_t*const data=(MTgomory_ccbk_t*)cbhandle;
00586   register int i;
00587   int rval = 0, action = 0;
00588   CPXLPptr lp = 0;
00589   int ncols,nrows; 
00590   double dtmp;
00591   double *base, *f;
00592   int cur_tabrow = 0;
00593   /* set default action */ 
00594   MESSAGE(MT_VERB_LVL+100,"Entering");
00595   *useraction_p = CPX_CALLBACK_DEFAULT; 
00596   /* process callback info this control the call back behavior. */
00597   rval = MTccbk_info_process(&(data->info), env, cbdata, wherefrom, &action);
00598   CHECKRVALG(rval,CLEANUP);
00599   if(action) goto CLEANUP;
00600   /* get the LP */
00601   rval = CPXgetcallbacknodelp (env, cbdata, wherefrom, &lp); 
00602   MTcplexCHECKRVALG(env,rval,CLEANUP);
00603   /* if we are verbose, write the problem to a file */
00604   if(MT_VERB_LVL +5 <= DEBUG)
00605   {
00606     rval = CPXwriteprob(env,lp,"node_lp.lp.gz","LP");
00607     MTcplexCHECKRVALG(env,rval,CLEANUP);
00608   }
00609   /* and get its description */
00610   rval = MTlp_load_problem(env, &(data->lp), lp, cbdata, wherefrom);
00611   CHECKRVALG(rval,CLEANUP);
00612   /* get numrows and numcols */
00613   ncols = CPXgetnumcols(env,lp);
00614   nrows = CPXgetnumrows(env,lp);
00615   MESSAGE(MT_VERB_LVL+3,"Problem has %d vars and %d constr",ncols,nrows);
00616   /* resize */
00617   rval = EGmax(nrows,ncols);
00618   MTgomory_ccbk_rsz(data,nrows,ncols+nrows,rval);
00619   f = data->f;
00620   base = data->base;
00621   /* now we get the current X and slack */
00622   rval = MTget_sol(&(data->lp), env, lp, cbdata, wherefrom, &(data->sol));
00623   CHECKRVALG(rval,CLEANUP);
00624   /* if we are debugging, display variable information */
00625   if(MT_VERB_LVL + 10 <= DEBUG)
00626     MTdisplay_lp(&(data->lp), &(data->sol), stderr);
00627   /* get the best MT_CCBK_MAX_VARS tableaus */
00628   rval = MTget_best_k_tableau_rows(&(data->lp), &(data->sol), data->info.max_rows,
00629       env, lp, &(data->tb), data->f);
00630   CHECKRVALG(rval,CLEANUP);
00631   if(!data->tb.nrows) goto CLEANUP;
00632   /* now we should set base appropiatelly */
00633   for( i = data->tb.nrows ; i-- ;  )
00634   {
00635     base[i] = 1.0;
00636     if(MT_VERB_LVL+5 <= DEBUG)
00637     {
00638       fprintf(stderr,"f[%d] = %7.4lf base[%d] = %7.4lf\n",i,f[i],i,base[i]);
00639     }
00640   }
00641   /* and now we should call the cut-generator routine */
00642   cur_tabrow = data->tb.nrows;
00643   REDO:
00644   rval = MTgomoryCut(cur_tabrow, data->tb.ncols, data->tb.rowval, 
00645       data->tb.rowind, data->tb.rowbeg, data->f, base, 
00646       data->extcut, data->work);
00647   CHECKRVALG(rval,CLEANUP);
00648   /* now we conver the extended raw cut to complemented compact form */
00649   rval = MTraw_to_complemented_cut(&(data->lp), &(data->sol), data->extcut, &(data->cut));
00650   CHECKRVALG(rval,CLEANUP);
00651   /* if we are debugging to the max, we solve the current IP
00652    * with objective value the found cut, and check that everything 
00653    * is OK with the value */
00654   if(MT_DEBUG_LVL +10 <= DEBUG)
00655   {
00656     rval = MTcplex_test_cut(&(data->cut),&(data->lp));
00657     CHECKRVALG(rval,CLEANUP);
00658   }
00659   /* if ratio is too small, don't add the cut */
00660   dtmp = data->cut.max_abs/data->cut.min_abs;
00661   if(dtmp > MT_CCBK_MAX_RATIO)
00662   {
00663     MESSAGE(MT_VERB_LVL+1, "Discarded by ratio (%lf<= %lf)", 
00664         dtmp, MT_CCBK_MAX_RATIO);
00665     MTccbk_fail_ratio++;
00666     if(MTccbk_discarded_ratio > dtmp) 
00667       MTccbk_discarded_ratio = dtmp;
00668     goto CLEANUP;
00669   }
00670   /* If violation is too small, don't add the cut */
00671   if(data->cut.violation < MT_CCBK_MIN_VIO)
00672   { 
00673     if(MTccbk_discarded_violation < data->cut.violation) 
00674       MTccbk_discarded_violation = data->cut.violation;
00675     MTccbk_fail_violation++;
00676     MESSAGE(MT_VERB_LVL+1, "Discarded by violation (%lf<= %lf)", 
00677         data->cut.violation, MT_CCBK_MIN_VIO);
00678     goto CLEANUP;
00679   }
00680   /* now we are in shape to add the cut to CPLEX */
00681   rval = CPXcutcallbackaddlocal( env, cbdata, wherefrom, data->cut.nz, 
00682       data->cut.rhs, data->cut.sense , 
00683       data->cut.cutind, data->cut.cutval);
00684   MTcplexCHECKRVALG(env,rval,CLEANUP);
00685   data->info.added_cuts++;
00686   data->info.node_cuts++;
00687   MESSAGE(MT_VERB_LVL, "Node %zd LP:(%d,%d,%d) Constr: %d nz, %.2lf rhs,"
00688       " ratio %.2lf, violation (>0) %lf", data->info.node_cnt, 
00689       data->lp.nrows, data->lp.ncols, data->lp.nz, data->cut.nz, data->cut.rhs, 
00690       data->cut.max_abs/data->cut.min_abs, data->cut.violation);
00691   if(data->info.cut_style && cur_tabrow>2) 
00692   {
00693     cur_tabrow--;
00694     goto REDO;
00695   }
00696   CLEANUP:
00697   MESSAGE(MT_VERB_LVL+100,"Done, status %d",rval);
00698   return rval;
00699 }
00700 /* ========================================================================= */
00701 int MTread_cplex_options(CPXENVptr env, FILE * inputfile)
00702 {
00703   const char delim[3] = " :";
00704   const char coment[3] = "%#";
00705   char str[2048],
00706        *argv[1024];
00707   int argc,
00708       rval=0,
00709       todo = 1,
00710       itmp;
00711   double tmpd;
00712 
00713   while(todo)
00714   {
00715     rval = EGioReadLine (str, (size_t) 2048, inputfile);
00716     CHECKRVALG(rval,CLEANUP);
00717     TESTG((rval=(feof(inputfile)||ferror(inputfile))), CLEANUP, 
00718         "file ended prematurelly");
00719     EGioNParse (str, 1024, delim, coment, &argc, argv);
00720     /* skip empty lines */
00721     if (!argc) continue;
00722     /* for the moment, we discard this information */
00723     if (strncmp (argv[0], "NAME", (size_t) 33) == 0) continue;
00724     /* for the moment, we discard this information */
00725     else if (strncmp (argv[0], "COMMENT", (size_t) 33) == 0) continue;
00726     /* for the moment, check type equal to GEOMIN */
00727     else if (strncmp (argv[0], "TYPE", (size_t) 33) == 0)
00728     {
00729       TESTG((rval=(argc!=2)), CLEANUP, "TYPE has not 2 tokens");
00730       TESTG((rval=(strncmp (argv[1], "OPTIONFILE", (size_t) 33))), CLEANUP, 
00731           "File TYPE is not OPTIONFILE, instead %s", argv[1]);
00732     }
00733     else if (strncmp (argv[0], "MIRCUTS", (size_t)33)==0)
00734     {
00735       TESTG((rval=(argc!=2)), CLEANUP, "MIRCUTS has not 2 tokens");
00736       rval = CPXsetintparam (env, CPX_PARAM_MIRCUTS, atoi(argv[1]));
00737       MTcplexCHECKRVALG(env,rval,CLEANUP);
00738       fprintf(stdout,"Setting: MIRCUTS %d\n",atoi(argv[1]));
00739     }
00740     else if (strncmp (argv[0], "GUBCOVERS", (size_t)33)==0)
00741     {
00742       TESTG((rval=(argc!=2)), CLEANUP, "GUBCOVERS has not 2 tokens");
00743       rval = CPXsetintparam (env, CPX_PARAM_GUBCOVERS, atoi(argv[1]));
00744       MTcplexCHECKRVALG(env,rval,CLEANUP);
00745       fprintf(stdout,"Setting: GUBCOVERS %d\n",atoi(argv[1]));
00746     }
00747     else if (strncmp (argv[0], "IMPLBD", (size_t)33)==0)
00748     {
00749       TESTG((rval=(argc!=2)), CLEANUP, "IMPLBD has not 2 tokens");
00750       rval = CPXsetintparam (env, CPX_PARAM_IMPLBD, atoi(argv[1]));
00751       MTcplexCHECKRVALG(env,rval,CLEANUP);
00752       fprintf(stdout,"Setting: IMPLBD %d\n",atoi(argv[1]));
00753     }
00754     else if (strncmp (argv[0], "COVERS", (size_t)33)==0)
00755     {
00756       TESTG((rval=(argc!=2)), CLEANUP, "COVERS has not 2 tokens");
00757       rval = CPXsetintparam (env, CPX_PARAM_COVERS, atoi(argv[1]));
00758       MTcplexCHECKRVALG(env,rval,CLEANUP);
00759       fprintf(stdout,"Setting: COVERS %d\n",atoi(argv[1]));
00760     }
00761     else if (strncmp (argv[0], "ZEROHALFCUTS", (size_t)33)==0)
00762     {
00763       TESTG((rval=(argc!=2)), CLEANUP, "ZEROHALFCUTS has not 2 tokens");
00764       if(CPX_VERSION < 1100)
00765       {
00766         fprintf(stdout,"WARNING NotSetting ZEROHALFCUTS CPLEX VERSION %.2lg < 11.00",((double)CPX_VERSION)/100.0);
00767       }
00768       else
00769       {
00770         #ifdef CPX_PARAM_ZEROHALFCUTS
00771         rval = CPXsetintparam (env, CPX_PARAM_ZEROHALFCUTS, atoi(argv[1]));
00772         MTcplexCHECKRVALG(env,rval,CLEANUP);
00773         fprintf(stdout,"Setting: ZEROHALFCUTS %d\n",atoi(argv[1]));
00774         #endif
00775       }
00776     }
00777     else if (strncmp (argv[0], "CLIQUES", (size_t)33)==0)
00778     {
00779       TESTG((rval=(argc!=2)), CLEANUP, "CLIQUES has not 2 tokens");
00780       rval = CPXsetintparam (env, CPX_PARAM_CLIQUES, atoi(argv[1]));
00781       MTcplexCHECKRVALG(env,rval,CLEANUP);
00782       fprintf(stdout,"Setting: CLIQUES %d\n",atoi(argv[1]));
00783     }
00784     else if (strncmp (argv[0], "DISJCUTS", (size_t)33)==0)
00785     {
00786       TESTG((rval=(argc!=2)), CLEANUP, "DISJCUTS has not 2 tokens");
00787       rval = CPXsetintparam (env, CPX_PARAM_DISJCUTS, atoi(argv[1]));
00788       MTcplexCHECKRVALG(env,rval,CLEANUP);
00789       fprintf(stdout,"Setting: DISJCUTS %d\n",atoi(argv[1]));
00790     }
00791     else if (strncmp (argv[0], "FLOWPATHS", (size_t)33)==0)
00792     {
00793       TESTG((rval=(argc!=2)), CLEANUP, "FLOWPATHS has not 2 tokens");
00794       rval = CPXsetintparam (env, CPX_PARAM_FLOWPATHS, atoi(argv[1]));
00795       MTcplexCHECKRVALG(env,rval,CLEANUP);
00796       fprintf(stdout,"Setting: FLOWPATHS %d\n",atoi(argv[1]));
00797     }
00798     else if (strncmp (argv[0], "FLOWCOVERS", (size_t)33)==0)
00799     {
00800       TESTG((rval=(argc!=2)), CLEANUP, "FLOWCOVERS has not 2 tokens");
00801       rval = CPXsetintparam (env, CPX_PARAM_FLOWCOVERS, atoi(argv[1]));
00802       MTcplexCHECKRVALG(env,rval,CLEANUP);
00803       fprintf(stdout,"Setting: FLOWCOVERS %d\n",atoi(argv[1]));
00804     }
00805     else if (strncmp (argv[0], "FRACCUTS", (size_t)33)==0)
00806     {
00807       TESTG((rval=(argc!=2)), CLEANUP, "FRACCUTS has not 2 tokens");
00808       rval = CPXsetintparam (env, CPX_PARAM_FRACCUTS, atoi(argv[1]));
00809       MTcplexCHECKRVALG(env,rval,CLEANUP);
00810       fprintf(stdout,"Setting: FRACCUTS %d\n",atoi(argv[1]));
00811     }
00812     else if (strncmp (argv[0], "INTTOL", (size_t)33)==0)
00813     {
00814       TESTG((rval=(argc!=2)), CLEANUP, "INTTOL has not 2 tokens");
00815       rval = CPXsetdblparam (env, CPX_PARAM_EPINT, strtod(argv[1],0));
00816       MTcplexCHECKRVALG(env,rval,CLEANUP);
00817       fprintf(stdout,"Setting: INTTOL %.6lg\n",strtod(argv[1],0));
00818     }
00819     else if (strncmp (argv[0], "PROBE", (size_t)33) == 0)
00820     {
00821       TESTG((rval=(argc!=2)), CLEANUP, "PROBE has not 2 tokens");
00822       rval = CPXsetintparam (env, CPX_PARAM_PROBE, atoi(argv[1]));
00823       MTcplexCHECKRVALG(env,rval,CLEANUP);
00824       fprintf(stdout,"Setting: PROBE %d\n",atoi(argv[1]));
00825     }
00826     else if (strncmp (argv[0], "VARSEL", (size_t)33) == 0)
00827     {
00828       TESTG((rval=(argc!=2)), CLEANUP, "VARSEL has not 2 tokens");
00829       rval = CPXsetintparam (env, CPX_PARAM_VARSEL, atoi(argv[1]));
00830       MTcplexCHECKRVALG(env,rval,CLEANUP);
00831       fprintf(stdout,"Setting: VARSEL %d\n",atoi(argv[1]));
00832     }
00833     else if (strncmp (argv[0], "MIPDISPLAY", (size_t)33)==0)
00834     {
00835       TESTG((rval=(argc!=2)), CLEANUP, "MIPDISPLAY has not 2 tokens");
00836       rval = CPXsetintparam (env, CPX_PARAM_MIPDISPLAY, atoi(argv[1]));
00837       MTcplexCHECKRVALG(env,rval,CLEANUP);
00838       fprintf(stdout,"Setting: MIPDISPLAY %d\n",atoi(argv[1]));
00839     }
00840     else if (strncmp (argv[0], "MIPINTERVAL", (size_t)33)==0)
00841     {
00842       TESTG((rval=(argc!=2)), CLEANUP, "MIPINTERVAL has not 2 tokens");
00843       rval = CPXsetintparam (env, CPX_PARAM_MIPINTERVAL, atoi(argv[1]));
00844       MTcplexCHECKRVALG(env,rval,CLEANUP);
00845       fprintf(stdout,"Setting: MIPINTERVAL %d\n",atoi(argv[1]));
00846     }
00847     else if (strncmp (argv[0], "RINSHEUR", (size_t)33) == 0)
00848     {
00849       TESTG((rval=(argc!=2)), CLEANUP, "RINSHEUR has not 2 tokens");
00850       rval = CPXsetintparam (env, CPX_PARAM_RINSHEUR, atoi(argv[1]));
00851       MTcplexCHECKRVALG(env,rval,CLEANUP);
00852       fprintf(stdout,"Setting: RINSHEUR %d\n",atoi(argv[1]));
00853     }
00854     else if (strncmp (argv[0], "HEURFREQ", (size_t)33) == 0)
00855     {
00856       TESTG((rval=(argc!=2)), CLEANUP, "HEURFREQ has not 2 tokens");
00857       rval = CPXsetintparam (env, CPX_PARAM_HEURFREQ, atoi(argv[1]));
00858       MTcplexCHECKRVALG(env,rval,CLEANUP);
00859       fprintf(stdout,"Setting: HEURFREQ %d\n",atoi(argv[1]));
00860     }
00861     else if (strncmp (argv[0], "NODELIM", (size_t)33) == 0)
00862     {
00863       TESTG((rval=(argc!=2)), CLEANUP, "NODELIM has not 2 tokens");
00864       rval = CPXsetintparam (env, CPX_PARAM_NODELIM, atoi(argv[1]));
00865       MTcplexCHECKRVALG(env,rval,CLEANUP);
00866       fprintf(stdout,"Setting: NODELIM %d\n",atoi(argv[1]));
00867     }
00868     else if (strncmp (argv[0], "LBHEUR", (size_t)33) == 0)
00869     {
00870       TESTG((rval=(argc!=2)), CLEANUP, "LBHEUR has not 2 tokens");
00871       if(strncmp(argv[1], "ON", (size_t)33)==0)
00872       {
00873         rval = CPXsetintparam (env, CPX_PARAM_LBHEUR, CPX_ON);
00874         MTcplexCHECKRVALG(env,rval,CLEANUP);
00875       }
00876       else if(strncmp(argv[1], "OFF", (size_t)33)==0)
00877       {
00878         rval = CPXsetintparam (env, CPX_PARAM_LBHEUR, CPX_OFF);
00879         MTcplexCHECKRVALG(env,rval,CLEANUP);
00880       }
00881       else
00882       {
00883         TESTG((rval=1),CLEANUP,"LBHEUR should be set to either ON or OFF");
00884       }
00885       fprintf(stdout,"Setting: LBHEUR %s\n",argv[1]);
00886     }
00887     else if (strncmp (argv[0], "SCREEN", (size_t)33) == 0)
00888     {
00889       TESTG((rval=(argc!=2)), CLEANUP, "SCREEN has not 2 tokens");
00890       if(strncmp(argv[1], "ON", (size_t)33)==0)
00891       {
00892         rval = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);
00893         MTcplexCHECKRVALG(env,rval,CLEANUP);
00894       }
00895       else if(strncmp(argv[1], "OFF", (size_t)33)==0)
00896       {
00897         rval = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_OFF);
00898         MTcplexCHECKRVALG(env,rval,CLEANUP);
00899       }
00900       else
00901       {
00902         TESTG((rval=1),CLEANUP,"SCREEN should be set to either ON or OFF");
00903       }
00904       fprintf(stdout,"Setting: SCREEN %s\n",argv[1]);
00905     }
00906     else if (strncmp (argv[0], "THREADS",  (size_t) 33) == 0)
00907     {
00908       TESTG((rval=(argc!=2)), CLEANUP, "THREADS has not 2 tokens");
00909       TESTG((rval=(atoi(argv[1])<1)), CLEANUP, "THREADS must be > 0");
00910       rval = CPXsetintparam(env,CPX_PARAM_THREADS,atoi(argv[1]));
00911       MTcplexCHECKRVALG(env,rval,CLEANUP);
00912       MTcplexCHECKRVALG(env,rval,CLEANUP);
00913       fprintf(stdout,"Setting: THREADS %d\n",atoi(argv[1]));
00914     }
00915     else if (strncmp (argv[0], "PRESOLVE",  (size_t) 33) == 0)
00916     {
00917       TESTG((rval=(argc!=2)), CLEANUP, "PRESOLVE has not 2 tokens");
00918       if (strncmp (argv[1], "ON", (size_t) 33) == 0)
00919       {
00920         rval = CPXsetintparam(env,CPX_PARAM_PREIND,CPX_ON);
00921         MTcplexCHECKRVALG(env,rval,CLEANUP);
00922         fprintf(stdout,"Setting: PRESOLVE ON\n");
00923       }else if  (strncmp (argv[1], "OFF", (size_t) 33) == 0)
00924       {
00925         rval = CPXsetintparam(env,CPX_PARAM_PREIND,CPX_OFF);
00926         MTcplexCHECKRVALG(env,rval,CLEANUP);
00927         fprintf(stdout,"Setting: PRESOLVE OFF\n");
00928       }
00929       else
00930       {
00931         rval = 1;
00932         fprintf(stdout, "PRESOLVE must be either ON or OFF\n");
00933         goto CLEANUP;
00934       }
00935     }
00936     else if (strncmp (argv[0], "MIPEMPHASIS", (size_t)33) == 0)
00937     {
00938       TESTG((rval=(argc!=2)), CLEANUP, "MIPEMPHASIS has not 2 tokens");
00939       itmp = atoi(argv[1]);
00940       TESTG((rval=(itmp<0)), CLEANUP, "MIPEMPHASIS must be >= 0");
00941       TESTG((rval=(itmp>4)), CLEANUP, "MIPEMPHASIS must be <= 4");
00942       rval=CPXsetintparam(env , CPX_PARAM_MIPEMPHASIS , itmp);
00943       MTcplexCHECKRVALG(env,rval,CLEANUP);
00944       fprintf(stdout,"Setting: MIPEMPHASIS to %d\n",itmp);
00945     }
00946     else if (strncmp (argv[0], "TIMELIMIT", (size_t) 33) == 0)
00947     {
00948       TESTG((rval=(argc!=2)), CLEANUP, "TIMELIMIT has not 2 tokens");
00949       TESTG((rval=(strtod(argv[1],0)<1)), CLEANUP, "TIMELIMIT must be > 0");
00950       rval=CPXsetintparam(env , CPX_PARAM_CLOCKTYPE , 1);
00951       MTcplexCHECKRVALG(env,rval,CLEANUP);
00952       rval=CPXsetdblparam(env , CPX_PARAM_TILIM , strtod(argv[1],0));
00953       MTcplexCHECKRVALG(env,rval,CLEANUP);
00954       fprintf(stdout,"Setting: TIMELIMIT %.2lf seconds (CPU time)\n",strtod(argv[1],0));
00955     }
00956     else if (strncmp (argv[0], "CUTSFACTOR", (size_t) 33) == 0)
00957     {
00958       TESTG((rval=(argc!=2)), CLEANUP, "CUTSFACTOR has not 2 tokens");
00959       tmpd = strtod (argv[1], 0);
00960       TESTG((rval=(tmpd<0.0)), CLEANUP, "CUTSFACTOR must be > 0,0");
00961       rval=CPXsetdblparam(env , CPX_PARAM_CUTSFACTOR , tmpd);
00962       MTcplexCHECKRVALG(env,rval,CLEANUP);
00963       fprintf(stdout,"Setting: CUTSFACTOR %f\n",tmpd);
00964     }
00965     else if (strncmp (argv[0], "GAP", (size_t) 33) == 0)
00966     {
00967       TESTG((rval=(argc!=2)), CLEANUP, "GAP has not 2 tokens");
00968       tmpd = strtod (argv[1], 0);
00969       TESTG((rval=(tmpd<0.0)), CLEANUP, "GAP must be > 0,0");
00970       rval=CPXsetdblparam(env , CPX_PARAM_EPGAP , tmpd);
00971       MTcplexCHECKRVALG(env,rval,CLEANUP);
00972       fprintf(stdout,"Setting: GAP %f\n",tmpd);
00973     }
00974     else if (strncmp (argv[0], "OBJ_UPLIM", (size_t) 33) == 0)
00975     {
00976       TESTG((rval=(argc!=2)), CLEANUP, "OBJ_UPLIM has not 2 tokens");
00977       tmpd = strtod (argv[1], 0);
00978       rval=CPXsetdblparam(env , CPX_PARAM_CUTUP , tmpd);
00979       MTcplexCHECKRVALG(env,rval,CLEANUP);
00980       fprintf(stdout,"Setting: OBJ_UPLIM %f\n",tmpd);
00981     }
00982     else if (strncmp (argv[0], "OBJ_LOLIM", (size_t) 33) == 0)
00983     {
00984       TESTG((rval=(argc!=2)), CLEANUP, "OBJ_LOLIM has not 2 tokens");
00985       tmpd = strtod (argv[1], 0);
00986       rval=CPXsetdblparam(env , CPX_PARAM_CUTLO , tmpd);
00987       MTcplexCHECKRVALG(env,rval,CLEANUP);
00988       fprintf(stdout,"Setting: OBJ_LOLIM %f\n",tmpd);
00989     }
00990     else if (strncmp (argv[0], "ROOT_ALG",  (size_t) 33) == 0)
00991     {
00992       TESTG((rval=(argc!=2)), CLEANUP, "ROOT_ALG has not 2 tokens");
00993       if (strncmp (argv[1], "PRIMAL", (size_t) 33) == 0)
00994       {
00995         rval=CPXsetintparam(env,CPX_PARAM_STARTALG,CPX_ALG_PRIMAL);
00996         MTcplexCHECKRVALG(env,rval,CLEANUP);
00997         fprintf(stdout,"Setting: ROOT_ALG PRIMAL\n");
00998       }
00999       else if (strncmp (argv[1], "DUAL", (size_t) 33) == 0)
01000       {
01001         rval=CPXsetintparam(env,CPX_PARAM_STARTALG,CPX_ALG_DUAL);
01002         MTcplexCHECKRVALG(env,rval,CLEANUP);
01003         fprintf(stdout,"Setting: ROOT_ALG DUAL\n");
01004       }
01005       else if (strncmp (argv[1], "BARRIER", (size_t) 33) == 0)
01006       {
01007         rval=CPXsetintparam(env,CPX_PARAM_STARTALG,CPX_ALG_BARRIER);
01008         MTcplexCHECKRVALG(env,rval,CLEANUP);
01009         fprintf(stdout,"Setting: ROOT_ALG BARRIER\n");
01010       }
01011       else 
01012       {
01013         rval = 1;
01014         fprintf(stdout, "ROOT_ALG must be PRIMAL, DUAL or BARRIER\n");
01015         goto CLEANUP;
01016       }
01017     }
01018     else if (strncmp (argv[0], "SUB_ALG",  (size_t) 33) == 0)
01019     {
01020       TESTG((rval=(argc!=2)), CLEANUP, "SUB_ALG has not 2 tokens");
01021       if (strncmp (argv[1], "PRIMAL", (size_t) 33) == 0)
01022       {
01023         rval=CPXsetintparam(env,CPX_PARAM_SUBALG,CPX_ALG_PRIMAL);
01024         MTcplexCHECKRVALG(env,rval,CLEANUP);
01025         fprintf(stdout,"Setting: SUB_ALG PRIMAL\n");
01026       }
01027       else if (strncmp (argv[1], "DUAL", (size_t) 33) == 0)
01028       {
01029         rval=CPXsetintparam(env,CPX_PARAM_SUBALG,CPX_ALG_DUAL);
01030         MTcplexCHECKRVALG(env,rval,CLEANUP);
01031         fprintf(stdout,"Setting: SUB_ALG DUAL\n");
01032       }
01033       else if (strncmp (argv[1], "BARRIER", (size_t) 33) == 0)
01034       {
01035         rval=CPXsetintparam(env,CPX_PARAM_SUBALG,CPX_ALG_BARRIER);
01036         MTcplexCHECKRVALG(env,rval,CLEANUP);
01037         fprintf(stdout,"Setting: SUB_ALG BARRIER\n");
01038       }
01039       else 
01040       {
01041         rval = 1;
01042         fprintf(stdout, "SUB_ALG must be PRIMAL, DUAL or BARRIER\n");
01043         goto CLEANUP;
01044       }
01045     }
01046     else if (strncmp (argv[0], "USE_TREE",  (size_t) 33) == 0)
01047     {
01048       //TESTG((rval=(argc!=2)), CLEANUP, "USE_TREE has not 2 tokens");
01049       //rval = CPXsetintparam(env,CPX_PARAM_ADVIND,CPX_ON);
01050       //MTcplexCHECKRVALG(env,rval,CLEANUP);
01051       //fprintf(stdout,"About to read tree\n");
01052       //rval = CPXreadcopytree(env,prob,argv[1]);
01053       //MTcplexCHECKRVALG(env,rval,CLEANUP);
01054       //fprintf(stdout,"Setting: USE_TREE %s\n",argv[1]);
01055       ;
01056     }
01057     else if (strncmp (argv[0], "EOF", (size_t) 33) == 0)
01058     {
01059       todo = 0;
01060     }
01061     /* default case */
01062     else
01063     {
01064       rval = 1;
01065       printf ("Unknown token %s\n", argv[0]);
01066       goto CLEANUP;
01067     }
01068     /* end loop */
01069   }
01070   CLEANUP:
01071   return rval;
01072 }
01073 /* ========================================================================= */
01074 void MTsol_rsz( MTsol_t*const MTsol,
01075     const int cnz,
01076     const int rnz)
01077 {
01078   if(MTsol->cspace < cnz)
01079   {
01080     MTsol->cspace = cnz;
01081     MTsol->x = EGrealloc(MTsol->x,sizeof(double)*((size_t)cnz));
01082     MTsol->pseudocost = EGrealloc(MTsol->pseudocost,sizeof(double)*((size_t)cnz));
01083   }
01084   if(MTsol->rspace < rnz)
01085   {
01086     MTsol->rspace = rnz;
01087     MTsol->slack = EGrealloc(MTsol->slack,sizeof(double)*((size_t)rnz));
01088   }
01089   return;
01090 }
01091 /* ========================================================================= */
01092 void MTrowmatrix_rsz(MTrowmatrix_t*const MTrm,
01093     const int rsz,
01094     const int tsz)
01095 {
01096   if(MTrm->rspace < rsz)
01097   {
01098     MTrm->rspace = rsz;
01099     MTrm->rowbeg = EGrealloc(MTrm->rowbeg, sizeof(int)*((size_t)(rsz+1)));
01100   }
01101   if(MTrm->nzspace < tsz)
01102   {
01103     MTrm->nzspace = tsz;
01104     MTrm->rowval = EGrealloc(MTrm->rowval,sizeof(double)*((size_t)tsz));
01105     MTrm->rowind = EGrealloc(MTrm->rowind,sizeof(int)*((size_t)tsz));
01106   }
01107   return;
01108 }
01109 /* ========================================================================= */
01110 void MTccbk_info_display(const MTccbk_info_t*const info,FILE*out)
01111 {
01112   fprintf(out,"Root LP Value : %.6lf\n",info->root_val);
01113   fprintf(out,"Added cuts    : %zd\n",info->added_cuts);
01114 }
01115 /* ========================================================================= */
01116 int MTccbk_info_process(MTccbk_info_t*const info,
01117     CPXCENVptr env,
01118     void*cbdata,
01119     int wherefrom,
01120     int*const action)
01121 {
01122   int rval = 0;
01123   int node_cnt;
01124   *action = 0;
01125   rval = CPXgetcallbackinfo(env, cbdata, wherefrom, 
01126       CPX_CALLBACK_INFO_NODE_COUNT, &node_cnt);
01127   MTcplexCHECKRVALG(env,rval,CLEANUP);
01128   if(((size_t)node_cnt) != info->node_cnt) info->node_cuts = 0;
01129   info->node_cnt = (size_t)node_cnt;
01130   if(node_cnt == 0)
01131   {
01132     rval = CPXgetcallbacknodeobjval(env, cbdata, wherefrom, &(info->root_val));
01133     MTcplexCHECKRVALG(env,rval,CLEANUP);
01134   }
01135   if(!info->use_cuts)
01136   {
01137     MESSAGE(MT_VERB_LVL,"Cuts disabled from beggining");
01138     *action = 1;
01139     goto CLEANUP;
01140   }
01141   if(info->root_only && node_cnt)
01142   {
01143     MESSAGE(MT_VERB_LVL,"node count limit reached");
01144     *action = 1;
01145     goto CLEANUP;
01146   }
01147   if(info->node_cuts >= info->max_cuts)
01148   {
01149     MESSAGE(MT_VERB_LVL,"Cut limit reached");
01150     *action = 1;
01151     goto CLEANUP;
01152   }
01153   CLEANUP:
01154   return rval;
01155 }
01156 /* ========================================================================= */
01157 void MTsol_clear(MTsol_t*const MTsol)
01158 {
01159   EGfree(MTsol->x);
01160   EGfree(MTsol->slack);
01161   EGfree(MTsol->pseudocost);
01162   memset(MTsol,0,sizeof(MTsol_t));
01163   return;
01164 }
01165 /* ========================================================================= */
01166 void MTrowmatrix_clear(MTrowmatrix_t*const MTrm)
01167 {
01168   EGfree(MTrm->rowbeg);
01169   EGfree(MTrm->rowind);
01170   EGfree(MTrm->rowval);
01171   memset(MTrm,0,sizeof(MTrowmatrix_t));
01172   return;
01173 }
01174 /* ========================================================================= */
01175 void MTlp_clear(MTlp_t*const MTlp)
01176 {
01177   EGfree(MTlp->sense);
01178   EGfree(MTlp->rhs);
01179   EGfree(MTlp->matbeg);
01180   EGfree(MTlp->rtype);
01181   EGfree(MTlp->rstat);
01182   EGfree(MTlp->ctype);
01183   EGfree(MTlp->cstat);
01184   EGfree(MTlp->matind);
01185   EGfree(MTlp->matval);
01186   EGfree(MTlp->ub);
01187   EGfree(MTlp->lb);
01188   #if MT_CCBK_USE_NAMES
01189   EGfree(MTlp->colnames);
01190   EGfree(MTlp->rownames);
01191   EGfree(MTlp->ccspace);
01192   EGfree(MTlp->rcspace);
01193   MTlp->ccsz = MTlp->rcsz = 0;
01194   #endif
01195   MTlp->cspace = MTlp->rspace = MTlp->nzspace = 0;
01196   return;
01197 }
01198 /* ========================================================================= */
01199 void MTcplex_display_tableau( FILE* file,
01200     const MTlp_t*const lp, 
01201     const double*const tableau)
01202 {
01203   const double*const slacks = tableau + lp->ncols;
01204   register int i;
01205   for( i = 0 ; i < lp->ncols ; i++)
01206   {
01207     if(EGabs(tableau[i]) < MT_ZERO_TOL) continue;
01208     fprintf(file, "%s%.6lg ",tableau[i] > 0 ? "+":"-", EGabs(tableau[i]));
01209     fprintf(file,"[%c]", lp->cstat[i]+'A');
01210     #if MT_CCBK_USE_NAMES
01211     fprintf(file,"%s ", lp->colnames[i]);
01212     #else
01213     fprintf(file,"x_%d ", i);
01214     #endif
01215   }
01216   for( i = 0 ; i < lp->nrows ; i++)
01217   {
01218     if(EGabs(slacks[i]) < MT_ZERO_TOL) continue;
01219     fprintf(file, "%s%.6lg ", slacks[i] > 0 ? "+":"-", EGabs(slacks[i]));
01220     fprintf(file,"[%c]", lp->rstat[i]+'A');
01221     #if MT_CCBK_USE_NAMES
01222     fprintf(file,"%s ", lp->rownames[i]);
01223     #else
01224     fprintf(file,"s_%d ", i);
01225     #endif
01226   }
01227   return;
01228 }
01229 /* ========================================================================= */
01230 void MTcplex_display_compress_tableau(FILE* file,
01231     const MTlp_t*const lp,
01232     const double*const matval,
01233     const int*const matind,
01234     const int nz)
01235 {
01236   register int i;
01237   int curind=0;
01238   for( i = 0 ; i < nz ; i++)
01239   {
01240     curind = matind[i];
01241     if(curind < lp->ncols)
01242     {
01243       fprintf(file, "%s%.5lf %s", matval[i]>0 ? "+":"-", EGabs(matval[i]),
01244           lp->cstat[curind] == CPX_AT_LOWER ? "": "~");
01245       #if MT_CCBK_USE_NAMES
01246       fprintf(file,"%s ",lp->colnames[curind]);
01247       #else
01248       fprintf(file,"x_%d ",curind);
01249       #endif
01250     }
01251     else
01252     {
01253       curind -= lp->ncols;
01254       fprintf(file, "%s%.5lf %s", matval[i]>0 ? "+":"-", EGabs(matval[i]),
01255           lp->rstat[curind] == CPX_AT_LOWER ? "": "~");
01256       #if MT_CCBK_USE_NAMES
01257       fprintf(file,"%s ",lp->rownames[curind]);
01258       #else
01259       fprintf(file,"s_%d ",curind);
01260       #endif
01261     }
01262   }
01263 }
01264 /* ========================================================================= */
01265 int MTget_best_k_integer_tableau_rows(const MTlp_t*const lp,
01266     const MTsol_t*const sol,
01267     const int max_tableau,
01268     CPXCENVptr env,
01269     CPXLPptr CPXlp,
01270     MTrowmatrix_t*const tb,
01271     double*const f)
01272 {
01273   const double*const x = sol->x;
01274   const double*const slack = sol->slack;
01275   const double*const pseudocost = sol->pseudocost;
01276   const int nrows = lp->nrows;
01277   const int ncols = lp->ncols;
01278   const int*const cstat = lp->cstat;
01279   register int i = 0;
01280   double dtmp,*binvrow = EGsMalloc(double,nrows),*tableau = 0, minratio=DBL_MAX, ratio;
01281   int rval = 0, nfrac = 0, brow = 0, tnz = 0, nfrac2 = 0, col = 0, status;
01282   int* rowmap = EGsMalloc(int,nrows), *tabind = 0;
01283   MESSAGE(MT_VERB_LVL,"Entering");
01284   tb->nz = tb->nrows = 0;
01285   tb->nrows = 0;
01286   tb->ncols = ncols+nrows;
01287   MTrowmatrix_rsz(tb,nrows,0);
01288   tb->rowbeg[0] = 0;
01289   /* compute all fractional variables */
01290   rval = MTcompute_integer_vars(lp, x, rowmap, &nfrac2, MT_CCBK_MIN_FRAC);
01291   CHECKRVALG(rval,CLEANUP);
01292   /* sort integer variables according to pseudocosts */
01293   dbl_EGutilPermSort2(((size_t)nfrac2),rowmap, pseudocost);
01294   /* now we start looping on all variables and check that they are basic
01295    * and that their asociated tableau is correct */
01296   for( i = 0 ; i < nfrac2 && nfrac < max_tableau ; i++ )
01297   {
01298     col = rowmap[i];
01299     /* check that the variable is basic */
01300     if(cstat[col] != CPX_BASIC) continue;
01301     /* get asociated basic row */
01302     rval = CPXgetijrow(env, CPXlp, -1, col, &brow);
01303     MTcplexCHECKRVALG(env,rval,CLEANUP);
01304     rval = CPXbinvrow(env, CPXlp, brow, binvrow);
01305     MTcplexCHECKRVALG(env,rval,CLEANUP);
01306     MTrowmatrix_rsz(tb,nrows,tb->nz+tb->ncols);
01307     tabind = tb->rowind+tb->nz;
01308     tableau = tb->rowval+tb->nz;
01309     rval = MTcplex_binv_to_tableau(lp, binvrow, tableau);
01310     CHECKRVALG(rval,CLEANUP);
01311     if(MT_VERB_LVL+10 < DEBUG) 
01312       MTcplex_display_tableau(stderr, lp, tableau);
01313     /* check that the asociated tableau is right */
01314     rval = MTcheckTabrow(lp,tableau,&status,&ratio,&tnz);
01315     CHECKRVALG(rval,CLEANUP);
01316     /* discard bad tableau rows */
01317     if(status) continue;
01318     if(minratio > ratio) minratio = ratio;
01319     /* store computed tableau */
01320     /* now we know that the tableau is good, and its length, so we save it */
01321     MTcompute_f(lp, tableau, x, slack, f+nfrac);
01322     tb->rowbeg[nfrac+1] = tb->nz+tnz;
01323     /* now we compress the tableau */
01324     rval = MTcompressTabRow(lp, sol, tableau, tabind, f+nfrac, &tnz);
01325     CHECKRVALG(rval,CLEANUP);
01326     MESSAGE(MT_VERB_LVL+5,"Tableau %d has %d non-zeros",nfrac,tnz);
01327     TESTG((rval=(tb->rowbeg[nfrac+1] !=tnz+tb->nz)), CLEANUP,"Wrong tnz count,"
01328         " found %d, should be %d", tnz + tb->nz, tb->rowbeg[nfrac+1]);
01329     /* up to now we have (but for integer factors) ax = f, and we want 
01330      * Z = f + Ax and f should be about zero */
01331     f[nfrac] = -f[nfrac];
01332     /* now we should check that |f[nfrac]|  approx 0.0 */
01333     dtmp = EGabs(f[nfrac]);
01334     if(dtmp > MT_CCBK_MIN_FRAC)
01335     {
01336       MTccbk_fail_tableau_frac++;
01337       MESSAGE(MT_VERB_LVL, "Discarding tableau %d, because f[%d] = %.4lf", 
01338           nfrac, nfrac, f[nfrac]);
01339       continue;
01340     }
01341     /* now we passed all test for the tableau, so we save it */
01342     nfrac++;
01343     tb->nrows++;
01344     tb->nz += tnz;
01345     /* done looking for fractional rows in the tableau */
01346   }
01347   CLEANUP:
01348   EGfree(rowmap);
01349   EGfree(binvrow);
01350   MESSAGE(MT_VERB_LVL,"Done, status %d, found %d integer tableau rows, minratio %.9lg",rval, nfrac,minratio);
01351   return rval;
01352 }
01353 /* ========================================================================= */
01354 int MTget_best_k_tableau_rows(const MTlp_t*const lp,
01355     const MTsol_t*const sol,
01356     const int max_tableau,
01357     CPXCENVptr env,
01358     CPXLPptr CPXlp,
01359     MTrowmatrix_t*const tb,
01360     double*const f)
01361 {
01362   const double*const x = sol->x;
01363   const double*const slack = sol->slack;
01364   const double*const pseudocost = sol->pseudocost;
01365   const int nrows = lp->nrows;
01366   const int ncols = lp->ncols;
01367   const int*const cstat = lp->cstat;
01368   register int i = 0;
01369   double dtmp,*binvrow = EGsMalloc(double,nrows),*tableau = 0, minratio=DBL_MAX, ratio;
01370   int rval = 0, nfrac = 0, brow = 0, tnz = 0, nfrac2 = 0, col = 0, status;
01371   int* rowmap = EGsMalloc(int,nrows), *tabind = 0;
01372   MESSAGE(MT_VERB_LVL,"Entering");
01373   tb->nz = tb->nrows = 0;
01374   tb->ncols = ncols+nrows;
01375   MTrowmatrix_rsz(tb,nrows,0);
01376   tb->rowbeg[0] = 0;
01377   /* compute all fractional variables */
01378   rval = MTcompute_fractional_vars(lp, x, rowmap, &nfrac2, MT_CCBK_MIN_FRAC);
01379   CHECKRVALG(rval,CLEANUP);
01380   /* sort fractional variables according to pseudocosts */
01381   dbl_EGutilPermSort2(((size_t)nfrac2),rowmap, pseudocost);
01382   /* now we start looping on all variables and check that they are fractional
01383    * and that their asociated tableau is correct */
01384   for( i = 0 ; i < nfrac2 && nfrac < max_tableau ; i++ )
01385   {
01386     col = rowmap[i];
01387     /* check that the variable is basic */
01388     TESTG((rval=(cstat[col] != CPX_BASIC)),CLEANUP, "Variable[%d] = %lf, status"
01389         " %d is not basic (%d)", col, x[col], cstat[col], CPX_BASIC); 
01390     /* get asociated basic row */
01391     rval = CPXgetijrow(env, CPXlp, -1, col, &brow);
01392     MTcplexCHECKRVALG(env,rval,CLEANUP);
01393     rval = CPXbinvrow(env, CPXlp, brow, binvrow);
01394     MTcplexCHECKRVALG(env,rval,CLEANUP);
01395     MTrowmatrix_rsz(tb,nrows,tb->nz+tb->ncols);
01396     tableau = tb->rowval+tb->nz;
01397     tabind = tb->rowind+tb->nz;
01398     rval = MTcplex_binv_to_tableau(lp, binvrow, tableau);
01399     CHECKRVALG(rval,CLEANUP);
01400     if(MT_VERB_LVL+10 < DEBUG) 
01401       MTcplex_display_tableau(stderr, lp, tableau);
01402     /* check that the asociated tableau is right */
01403     rval = MTcheckTabrow(lp,tableau,&status,&ratio,&tnz);
01404     CHECKRVALG(rval,CLEANUP);
01405     /* discard bad tableau rows */
01406     if(status) continue;
01407     if(minratio > ratio) minratio = ratio;
01408     /* store computed tableau */
01409     /* now we know that the tableau is good, and its length, so we save it */
01410     MTcompute_f(lp, tableau, x, slack, f+nfrac);
01411     tb->rowbeg[nfrac+1] = tb->nz+tnz;
01412     /* now we compress the tableau */
01413     rval = MTcompressTabRow(lp, sol, tableau, tabind, f+nfrac, &tnz);
01414     CHECKRVALG(rval,CLEANUP);
01415     MESSAGE(MT_VERB_LVL+5,"Tableau %d has %d non-zeros",nfrac,tnz);
01416     /* up to now we have (but for integer factors) ax = f, and we want 
01417      * Z = f + 0.5 * 1I + Ax */
01418     f[nfrac] = -f[nfrac];
01419     f[nfrac] -= (floor(f[nfrac])+0x1p-1);
01420     /* now we should check that -0.5 < f[nfrac] < 0.5 */
01421     dtmp = EGabs(f[nfrac])-0x1p-1;
01422     dtmp = EGabs(dtmp);
01423     if(dtmp < MT_ZERO_TOL)
01424     {
01425       MTccbk_fail_tableau_frac++;
01426       MESSAGE(MT_VERB_LVL, "Discarding tableau %d, because f[%d] = %.4lf", 
01427           nfrac, nfrac, f[nfrac]);
01428       continue;
01429     }
01430     /* now we passed all test for the tableau, so we save it */
01431     nfrac++;
01432     tb->nrows++;
01433     tb->nz += tnz;
01434     /* done looking for fractional rows in the tableau */
01435   }
01436   CLEANUP:
01437   EGfree(rowmap);
01438   EGfree(binvrow);
01439   MESSAGE(MT_VERB_LVL,"Done, status %d, found %d fractional tableau rows, minratio %.9lg",rval, nfrac,minratio);
01440   return rval;
01441 }
01442 /* ========================================================================= */
01443 int MTraw_to_uncomplemented_cut(const MTlp_t*const lp,
01444     const double*const extended,
01445     const int nactive,
01446     const int*const active,
01447     MTcut_t*const cut, 
01448     int*const status)
01449 {
01450   const int nvars = lp->ncols + lp->nrows;
01451   double*cutval = 0;
01452   int*cutind = 0;
01453   int nz = 0,pos;
01454   double min_abs = DBL_MAX, max_abs = 0.0, absval;
01455   register int i;
01456   MESSAGE(MT_VERB_LVL+3,"Entering");
01457   /* resize if needed */
01458   MTcut_rsz(cut,nvars);
01459   cutind = cut->cutind;
01460   cutval = cut->cutval;
01461   /* now we process the cut */
01462   if(!active)
01463   {
01464     i = nvars;
01465     while(i--)
01466     {
01467       absval = EGabs(extended[i]);
01468       if(absval < MT_ZERO_TOL) continue;
01469       cutind[nz] = i;
01470       cutval[nz++] = extended[i];
01471       if(absval < min_abs) min_abs = absval;
01472       if(absval > max_abs) max_abs = absval;
01473     }
01474   }
01475   else
01476   {
01477     i = nactive;
01478     while(i--)
01479     {
01480       pos = active[i];
01481       absval = EGabs(extended[pos]);
01482       if(absval < MT_ZERO_TOL) continue;
01483       cutind[nz] = pos;
01484       cutval[nz++] = extended[pos];
01485       if(absval < min_abs) min_abs = absval;
01486       if(absval > max_abs) max_abs = absval;
01487     }
01488   }
01489   cut->score = 1.0/max_abs;
01490   cut->nz = nz;
01491   cut->rhs = 1.0;
01492   cut->sense = MT_CPLEX_SENSE;
01493   cut->violation = 1.0;
01494   cut->min_abs = min_abs;
01495   cut->max_abs = max_abs;
01496   cut->form = 1;
01497   *status = 0;
01498   /*  if ratio is too small, discard the cut */
01499   absval = max_abs / min_abs;
01500   if(absval > MT_CCBK_MAX_RATIO)
01501   {
01502     MESSAGE(MT_VERB_LVL+1, "Discarded by ratio (%lf<= %lf)", 
01503             absval, MT_CCBK_MAX_RATIO);
01504     MTccbk_fail_ratio++;
01505     if(MTccbk_discarded_ratio > absval) 
01506       MTccbk_discarded_ratio = absval;
01507     *status = 1;
01508     goto CLEANUP;
01509   }
01510   /* if violation is too small, discard it */
01511   if(cut->violation < MT_CCBK_MIN_VIO)
01512   {
01513     if(MTccbk_discarded_violation < cut->violation) 
01514       MTccbk_discarded_violation = cut->violation;
01515     MTccbk_fail_violation++;
01516     MESSAGE(MT_VERB_LVL+1, "Discarded by violation (%lf<= %lf)", 
01517             cut->violation, MT_CCBK_MIN_VIO);
01518     *status = 1;
01519   }
01520   /* display cut */
01521   if(MT_VERB_LVL+9 <= DEBUG)
01522   {
01523     fprintf(stderr,"MTcut (extended): ");
01524     MTcplex_display_tableau(stderr,lp,extended);
01525     fprintf(stderr,">=1\n");
01526   }
01527   CLEANUP:
01528   MESSAGE(MT_VERB_LVL+3,"Done");
01529   return 0;
01530 }
01531 /* ========================================================================= */
01532 int MTuncomplemented_to_complemented_cut( const MTlp_t*const lp,
01533     const MTsol_t*const sol,
01534     MTcut_t*const cut)
01535 {
01536   const int ncols = lp->ncols;
01537   const int nrows = lp->nrows;
01538   double*extended = EGsMalloc(double,ncols+nrows);
01539   int rval = 0;
01540   register int i;
01541   /* build the extended cut */
01542   memset(extended,0,sizeof(double)*((size_t)(ncols+nrows)));
01543   for( i = cut->nz ; i-- ; ) extended[cut->cutind[i]] = cut->cutval[i];
01544   /* call raw to complemented cut */
01545   rval = MTraw_to_complemented_cut(lp,sol,extended,cut); 
01546   CHECKRVALG(rval,CLEANUP);
01547   CLEANUP:
01548   EGfree(extended);
01549   MESSAGE(MT_VERB_LVL+3,"Done");
01550   return rval;
01551 }
01552 /* ========================================================================= */
01553 int MTraw_to_complemented_cut( const MTlp_t*const lp,
01554     const MTsol_t*const sol,
01555     double*const extended,
01556     MTcut_t*const cut)
01557 {
01558   const int ncols = lp->ncols;
01559   const int nrows = lp->nrows;
01560   const int*const cstat = lp->cstat;
01561   const int*const rstat = lp->rstat;
01562   const char*const sense = lp->sense;
01563   const double*const rhs_arr = lp->rhs;
01564   const int*const matbeg = lp->matbeg;
01565   const int*const matind = lp->matind;
01566   const double*const matval = lp->matval;
01567   const double*const x = sol->x;
01568   const double*const slack = sol->slack;
01569   double*cutval = 0;
01570   int*cutind = 0;
01571   int rval = 0, nz = 0;
01572   double rhs = 1.0, dtmp, violation = 0.0, min_abs = DBL_MAX, max_abs = 0.0;
01573   register int i,j;
01574   /* resize if needed */
01575   MTcut_rsz(cut,lp->ncols);
01576   cut->form = 0;
01577   cutind = cut->cutind;
01578   cutval = cut->cutval;
01579   /* display cut */
01580   MESSAGE(MT_VERB_LVL+3,"Entering");
01581   if(MT_VERB_LVL+3 <= DEBUG)
01582   {
01583     fprintf(stderr,"MTcut (before compl): ");
01584     MTcplex_display_tableau(stderr,lp,extended);
01585     fprintf(stderr,">=1\n");
01586   }
01587   /* compute raw score */
01588   for(dtmp = 0.0, i = ncols + nrows; i-- ; )
01589   {
01590     if(fabs(extended[i])>MT_ZERO_TOL) 
01591       dtmp += extended[i]*extended[i];
01592   }
01593   cut->score = 1.0/dtmp;
01594   /* un-complement structural variables */
01595   if(MT_VERB_LVL + 5 <= DEBUG) fprintf(stderr,"rhs = %.4lf\n",rhs);
01596   for( i = ncols ; i-- ; )
01597   {
01598     if(EGabs(extended[i]) < MT_ZERO_TOL) continue;
01599     TESTG((rval=(cstat[i] == CPX_BASIC)),CLEANUP,"Imposible, basic variable with cutvalue %.9lf",extended[i]);
01600     if(cstat[i] == CPX_AT_UPPER) extended[i] *= -1.0;
01601     rhs += extended[i]*x[i];
01602     if(MT_VERB_LVL + 5 <= DEBUG)
01603     {
01604       fprintf(stderr,"rhs = %.3lf ",rhs);
01605       #if MT_CCBK_USE_NAMES
01606       fprintf(stderr,"%s\n",lp->colnames[i]);
01607       #else
01608       fprintf(stderr,"x_%d\n",i);
01609       #endif
01610     }
01611   }
01612   /*  un-complement logical variables */
01613   for( i = nrows; i-- ; )
01614   {
01615     if(EGabs(extended[i+ncols]) < MT_ZERO_TOL) continue;
01616     TESTG((rval=(rstat[i]==CPX_BASIC)), CLEANUP, "Imposible");
01617     if(rstat[i] == CPX_AT_UPPER) extended[i+ncols] *= -1.0;
01618     rhs += extended[i+ncols]*slack[i];
01619     if(MT_VERB_LVL + 5 <= DEBUG)
01620     {
01621       fprintf(stderr,"rhs = %.3lf ",rhs);
01622       #if MT_CCBK_USE_NAMES
01623       fprintf(stderr,"%s\n",lp->rownames[i]);
01624       #else
01625       fprintf(stderr,"s_%d\n",i);
01626       #endif
01627     }
01628   }
01629   /* now we substitute back each logical variable */
01630   for( i = nrows ; i-- ; )
01631   {
01632     if(EGabs(extended[i+ncols]) < MT_ZERO_TOL) continue;
01633     switch(sense[i])
01634     {
01635       case 'L':
01636       case 'R':
01637       case 'l':
01638       case 'r':
01639         dtmp = 1.0;
01640         break;
01641       case 'G':
01642       case 'g':
01643         dtmp = -1.0;
01644         break;
01645       default:
01646         FTESTG((rval = 1),CLEANUP,"Imposible!");
01647         break;
01648     }
01649     rhs -= dtmp*extended[i+ncols]*rhs_arr[i];
01650     for( j = matbeg[i] ; j < matbeg[i+1] ; j++)
01651       extended[matind[j]] -= dtmp*extended[i+ncols]*matval[j];
01652     extended[i+ncols] = 0;
01653     if(MT_VERB_LVL + 5 <= DEBUG)
01654     {
01655       fprintf(stderr,"rhs = %.3lf ",rhs);
01656       #if MT_CCBK_USE_NAMES
01657       fprintf(stderr,"complement %s\n",lp->rownames[i]);
01658       #else
01659       fprintf(stderr,"complement s_%d\n",i);
01660       #endif
01661     }
01662   }
01663   /* compute violation, minimum/maximum absolute values */
01664   violation = rhs;
01665   for( i = ncols ; i-- ; )
01666   {
01667     dtmp = EGabs(extended[i]);
01668     if(dtmp < MT_ZERO_TOL) continue;
01669     violation -= x[i]*extended[i];
01670     if(dtmp < min_abs) min_abs = dtmp;
01671     if(dtmp > max_abs) max_abs = dtmp;
01672   }
01673   /* re-scale so that min_abs = 1 */
01674   if(min_abs > 1.0)
01675   {
01676     rhs /= min_abs;
01677     violation /= min_abs;
01678     max_abs /= min_abs;
01679     for( i = ncols ; i-- ; ) extended[i] /= min_abs;
01680     min_abs = 1.0;
01681   }
01682   /* display the complemented cut */
01683   if(MT_VERB_LVL <= DEBUG)
01684   {
01685     fprintf(stderr,"min_abs %.9lf max_abs %.9lf MTcut (complemented): ",
01686         min_abs, max_abs);
01687     MTcplex_display_tableau(stderr,lp,extended);
01688     fprintf(stderr," >= %.9lf\n", rhs);
01689   }
01690   /* compress the coefficients */
01691   for( i = 0 ; i < ncols ; i++)
01692   {
01693     if(EGabs(extended[i]) < MT_ZERO_TOL) continue;
01694     cutval[nz] = extended[i];
01695     cutind[nz++] = i;
01696   }
01697   cut->nz = nz;
01698   cut->rhs = rhs;
01699   cut->sense = MT_CPLEX_SENSE;
01700   cut->violation = violation;
01701   cut->min_abs = min_abs;
01702   cut->max_abs = max_abs;
01703   if(MT_VERB_LVL+5 <= DEBUG)
01704   {
01705     fprintf(stderr,"MTcut (complemented+compresed): ");
01706     MTcplex_display_compress_tableau(stderr, lp, cutval, cutind, nz);
01707     fprintf(stderr," >= %.4lf\n", rhs);
01708   }
01709   CLEANUP:
01710   MESSAGE(MT_VERB_LVL+3,"Done");
01711   return rval;
01712 }
01713 /* ========================================================================= */
01714 int MTget_sol(const MTlp_t*const lp,
01715     CPXCENVptr env,
01716     CPXLPptr CPXlp,
01717     void*cbdata,
01718     int wherefrom,
01719     MTsol_t*const sol)
01720 {
01721   const int ncols = lp->ncols;
01722   const int nrows = lp->nrows;
01723   register int i;
01724   int rval = EGmax(ncols,nrows);
01725   double dtmp;
01726   MTsol_rsz(sol,ncols,rval);
01727   sol->ncols = ncols;
01728   sol->nrows = nrows;
01729   /* get the x solution */
01730   rval = CPXgetcallbacknodex(env,cbdata,wherefrom,sol->x,0,ncols-1);
01731   MTcplexCHECKRVALG(env,rval,CLEANUP);
01732   /* get pseudocost information */
01733   rval = CPXgetcallbackpseudocosts(env, cbdata, wherefrom, sol->pseudocost, sol->slack, 0, ncols-1);
01734   MTcplexCHECKRVALG(env,rval,CLEANUP);
01735   for( i = ncols ; i-- ; ) 
01736   {
01737     dtmp = sol->x[i] - floor(sol->x[i]);
01738     sol->pseudocost[i] = sol->pseudocost[i]*(1-dtmp)+dtmp*sol->slack[i];
01739     WARNING(fabs(sol->pseudocost[i])<MT_ZERO_TOL, "pseudocost %d is zero",i);
01740   }
01741   /* get the slack */
01742   rval = CPXslackfromx(env,CPXlp,sol->x,sol->slack);
01743   MTcplexCHECKRVALG(env,rval,CLEANUP);
01744   CLEANUP:
01745   return rval;
01746 }
01747 /* ========================================================================= */
01748 int MTcplex_test_cut( const MTcut_t*const cut,
01749     const MTlp_t*const lp)
01750 {
01751   const int ncols = lp->ncols;
01752   const int nrows = lp->nrows;
01753   const int*const cutind = cut->cutind;
01754   const double*const cutval = cut->cutval;
01755   int rval = 0;
01756   CPXENVptr env = 0;
01757   CPXLPptr curip = 0;
01758   double*obj = EGsMalloc(double,ncols);
01759   double curip_obj = 0.0;
01760   double curip_bound = 0.0,dtmp;
01761   int curip_stat = 0;
01762   register int i;
01763   MESSAGE(MT_VERB_LVL+3,"Entering");
01764   TESTG((rval = cut->form),CLEANUP,"Trying to test cut in complemented form!");
01765   for( i = ncols ; i-- ; ) obj[0] = 0;
01766   for( i = cut->nz ; i-- ; ) obj[cutind[i]] = cutval[i]; 
01767   env = CPXopenCPLEX(&rval);
01768   MTcplexCHECKRVALG(env,rval,CLEANUP);
01769   rval = CPXsetintparam(env,CPX_PARAM_PREIND,CPX_OFF);
01770   MTcplexCHECKRVALG(env,rval,CLEANUP);
01771   rval = CPXsetdblparam (env, CPX_PARAM_EPINT, MT_ZERO_TOL);
01772   MTcplexCHECKRVALG(env,rval,CLEANUP);
01773   rval = CPXsetdblparam(env , CPX_PARAM_EPAGAP , MT_ZERO_TOL);
01774   MTcplexCHECKRVALG(env,rval,CLEANUP);
01775   rval = CPXsetdblparam(env , CPX_PARAM_EPGAP , MT_ZERO_TOL);
01776   MTcplexCHECKRVALG(env,rval,CLEANUP);
01777   curip = CPXcreateprob (env, &rval, "node_ip");
01778   MTcplexCHECKRVALG(env,rval,CLEANUP);
01779   #if MT_CCBK_USE_NAMES
01780   rval = CPXnewcols(env, curip, lp->ncols, obj, lp->lb, 
01781       lp->ub, lp->ctype, lp->colnames);
01782   #else
01783   rval = CPXnewcols(env, curip, lp->ncols, obj, lp->lb, 
01784       lp->ub, lp->ctype, 0);
01785   #endif
01786   MTcplexCHECKRVALG(env,rval,CLEANUP);
01787   #if MT_CCBK_USE_NAMES
01788   rval = CPXaddrows(env, curip, 0, nrows, 
01789       lp->matbeg[nrows], lp->rhs, 
01790       lp->sense, lp->matbeg, lp->matind, 
01791       lp->matval, 0, lp->rownames);
01792   #else
01793   rval = CPXaddrows(env, curip, 0, nrows, 
01794       lp->matbeg[nrows], lp->rhs, 
01795       lp->sense, lp->matbeg, lp->matind, 
01796       lp->matval, 0, 0);
01797   #endif
01798   MTcplexCHECKRVALG(env,rval,CLEANUP);
01799   rval = CPXwriteprob(env,curip, "node_ip.lp.gz","lp");
01800   MTcplexCHECKRVALG(env,rval,CLEANUP);
01801   /* Optimize the problem and obtain solution */
01802   rval = CPXmipopt (env, curip);
01803   MTcplexCHECKRVALG(env,rval,CLEANUP);
01804   curip_stat = CPXgetstat (env, curip);
01805   switch(curip_stat)
01806   {
01807     case CPXMIP_OPTIMAL:
01808     case CPXMIP_OPTIMAL_TOL:
01809       rval = CPXgetobjval (env, curip, &curip_obj);
01810       MTcplexCHECKRVALG(env,rval,CLEANUP);
01811       rval = CPXgetbestobjval (env, curip, &curip_bound);
01812       MTcplexCHECKRVALG(env,rval,CLEANUP);
01813       /* now test that the value is right! */
01814       rval = (curip_obj+MT_ZERO_TOL  < cut->rhs);
01815       dtmp = fabs(cut->rhs);
01816       if(dtmp < 1.0) dtmp=1.0;
01817       FTESTG(rval, CLEANUP,
01818           "Found an infeasible cut! rhs = %.6lf, actual minimum "
01819           "%.6lf, difference %.6lg (%.3lf%%)", cut->rhs, curip_obj, 
01820           cut->rhs - curip_obj, 100*(curip_bound-cut->rhs)/dtmp);
01821       fprintf(stderr,"\t\tdiff %.6lg (%.3lg%%)\n ",curip_bound-cut->rhs, 
01822           100*(curip_bound-cut->rhs)/dtmp);
01823       break;
01824     case CPXMIP_INFEASIBLE:
01825       break;
01826     case CPXMIP_INForUNBD:
01827       break;
01828     case CPXMIP_UNBOUNDED:
01829       break;
01830     case CPXMIP_OPTIMAL_INFEAS:
01831       WARNING(1,"Optimal with unscaled infeasibilites");
01832       break;
01833     default: 
01834       rval = 1;
01835       MESSAGE(0, "Unexpected status %d for node_ip", curip_stat);
01836       goto CLEANUP;
01837   }
01838   /* close and free local information */
01839   rval = CPXcloseCPLEX(&env);
01840   MTcplexCHECKRVALG(env,rval,CLEANUP);
01841   CLEANUP:
01842   EGfree(obj);
01843   MESSAGE(MT_VERB_LVL+3,"Done, with status %d",rval);
01844   return rval;
01845 }
01846 /* ========================================================================= */
01847 #define __MT_NS_LVL 0
01848 int MTnode_solve( CPXCENVptr env, 
01849     void*cbdata,
01850     int wherefrom, 
01851     void*cbhandle, 
01852     int*useraction_p)
01853 {
01854   MTns_ccbk_t*const ns_data = (MTns_ccbk_t*)cbhandle;
01855   int nfixed=0, prevfixed=0, loop=0, ncols=0,nrows=0, rval=0, i=0, objsen=0, 
01856       nextcol=INT_MAX, *rstat=0, *cstat=0, *cind=0, lp_stat=-1, stat0=-1, nz=0, 
01857       surplus=0, matbeg=0, *matind=0;
01858   char*cfix=0, *rfix=0, *sense=0, ctmp='E';
01859   double*dj=0,*pi=0, *oobj=0, *lb=0, *ub=0, *x0=0, *x1=0, dtmp=0, *matval=0, 
01860     llim=0, ulim=0;
01861   CPXLPptr lp = 0;
01862   /* basic set-up */
01863   *useraction_p = CPX_CALLBACK_DEFAULT;
01864   ns_data->ncalls++;
01865   /* get LP and minimal data/memory */
01866   rval = CPXgetcallbacknodelp(env,cbdata,wherefrom,&lp);
01867   MTcplexCHECKRVALG(env,rval,CLEANUP);
01868   ncols = CPXgetnumcols(env,lp);
01869   nrows = CPXgetnumrows(env,lp);
01870   objsen = CPXgetobjsen(env,lp);
01871   if(objsen == CPX_MAX) objsen = 1;
01872   else objsen = -1;
01873   MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Entering, ncols %d, nrows %d, sense %d",
01874       ncols, nrows, objsen);
01875   MTns_ccbk_rsz(ns_data,ncols,nrows);
01876   rval = CPXgetdblparam(env, CPX_PARAM_OBJLLIM, &llim);
01877   MTcplexCHECKRVALG(env,rval,CLEANUP);
01878   rval = CPXgetdblparam(env, CPX_PARAM_OBJULIM, &ulim);
01879   MTcplexCHECKRVALG(env,rval,CLEANUP);
01880   rval = CPXsetdblparam(env, CPX_PARAM_OBJLLIM, -DBL_MAX);
01881   MTcplexCHECKRVALG(env,rval,CLEANUP);
01882   rval = CPXsetdblparam(env, CPX_PARAM_OBJULIM, DBL_MAX);
01883   MTcplexCHECKRVALG(env,rval,CLEANUP);
01884   cfix = ns_data->cfix;
01885   rfix = ns_data->rfix;
01886   memset(cfix,'F',sizeof(char)*((size_t)ncols));
01887   memset(rfix,'F',sizeof(char)*((size_t)nrows));
01888   dj = ns_data->dj;
01889   pi = ns_data->pi;
01890   lb = ns_data->lb;
01891   ub = ns_data->ub;
01892   x0 = ns_data->x0;
01893   x1 = ns_data->x1;
01894   matind = ns_data->matind;
01895   matval = ns_data->matval;
01896   memset(x0,0,sizeof(double)*((size_t)ncols));
01897   memset(x1,0,sizeof(double)*((size_t)ncols));
01898   sense = ns_data->sense;
01899   oobj = ns_data->oobj;
01900   rstat = ns_data->rstat;
01901   cstat = ns_data->cstat;
01902   nfixed = EGmax(ncols,nrows);
01903   cind = EGsMalloc(int,nfixed);
01904   for( i = nfixed ; i-- ; ) cind[i] = i;
01905   rval = CPXgetobj(env, lp, oobj, 0, ncols-1);
01906   MTcplexCHECKRVALG(env,rval,CLEANUP);
01907   rval = CPXgetub(env, lp, ub, 0, ncols-1);
01908   MTcplexCHECKRVALG(env,rval,CLEANUP);
01909   rval = CPXgetlb(env, lp, lb, 0, ncols-1);
01910   MTcplexCHECKRVALG(env,rval,CLEANUP);
01911   rval = CPXgetsense(env, lp, sense, 0, nrows-1);
01912   MTcplexCHECKRVALG(env,rval,CLEANUP);
01913   /* now we start optimization loop */
01914   do
01915   {
01916     /* solve and get duals and basis solution */
01917     rval = CPXdualopt(env,lp);
01918     MTcplexCHECKRVALG(env,rval,CLEANUP);
01919     lp_stat = CPXgetstat(env, lp);  
01920     MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Optimizing... lp_stat %d",lp_stat);
01921     rval = CPXgetbase(env, lp, cstat, rstat);
01922     MTcplexCHECKRVALG(env,rval,CLEANUP);
01923     /* check what status we receive */
01924     if(loop==0)
01925     {
01926       stat0 = lp_stat;
01927       if(stat0 != CPX_STAT_OPTIMAL) break;
01928       rval = CPXgetx(env,lp,x0,0,ncols-1);
01929       MTcplexCHECKRVALG(env,rval,CLEANUP);
01930     }
01931     /* depending on the status we proceed */
01932     TESTG((rval=(lp_stat!= CPX_STAT_OPTIMAL)), CLEANUP, "Must handle this case!");
01933     /* if optimal get reduced costs etc. */
01934     rval = CPXgetdj(env,lp,dj,0,ncols-1);
01935     MTcplexCHECKRVALG(env,rval,CLEANUP);
01936     rval = CPXgetpi(env,lp,pi,0,nrows-1);
01937     MTcplexCHECKRVALG(env,rval,CLEANUP);
01938     rval = CPXgetx(env,lp,x1,0,ncols-1);
01939     MTcplexCHECKRVALG(env,rval,CLEANUP);
01940     /* check the status and look for a lexicographicly minimal LP solution */
01941     nfixed = 0;
01942   SOFT_REDO:
01943     loop++;
01944     nextcol = INT_MAX;
01945     /* process structural variables */
01946     for( i = ncols ; i-- ; )
01947     {
01948       if(cfix[i] != 'F') continue;
01949       switch(cstat[i])
01950       {
01951         case CPX_AT_LOWER:
01952           if(dj[i]*objsen < -MT_ZERO_TOL)
01953           {
01954             nfixed++;
01955             cfix[i] = 'L';
01956             ctmp = 'U';
01957             rval = CPXchgbds(env, lp, 1, &i, &ctmp, lb+i);
01958             MTcplexCHECKRVALG(env,rval,CLEANUP);
01959             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing %d to %.4lf",i,lb[i]);
01960           }
01961           else
01962           {
01963             nextcol = EGmin(i, nextcol);
01964             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Nextcol %d rd = %.7lg L",i,dj[i]);
01965           }
01966           break;
01967         case CPX_AT_UPPER:
01968           if(dj[i]*objsen > MT_ZERO_TOL)
01969           {
01970             nfixed++;
01971             cfix[i] = 'U';
01972             ctmp = 'L';
01973             rval = CPXchgbds(env, lp, 1, &i, &ctmp, ub+i);
01974             MTcplexCHECKRVALG(env,rval,CLEANUP);
01975             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing %d to %.4lf",i,ub[i]);
01976           }
01977           else
01978           {
01979             nextcol = EGmin(i, nextcol);
01980             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Nextcol %d rd = %.7lg U",i,dj[i]);
01981           }
01982           break;
01983         case CPX_FREE_SUPER:
01984           nfixed++;
01985           cfix[i] = 'S';
01986           break;
01987         case CPX_BASIC:
01988           break;
01989         default:
01990           TESTG((rval=1),CLEANUP,"Imposible!");
01991           break;
01992       }
01993     }
01994     /* process logical variables */
01995     for( i = nrows ; i-- ; )
01996     {
01997       if(rfix[i] != 'F') continue;
01998       switch(rstat[i])
01999       {
02000         case CPX_AT_LOWER:
02001           if(fabs(pi[i]) > MT_ZERO_TOL)
02002           {
02003             rfix[i] = 'L';
02004             nfixed++;
02005             ctmp = 'E';
02006             rval = CPXchgsense(env, lp, 1, &i, &ctmp);
02007             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing row %d to equality",i);
02008             MTcplexCHECKRVALG(env,rval,CLEANUP);
02009           }
02010           else
02011           {
02012             nextcol = EGmin(nextcol, i+ncols);
02013             MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Nextcol %d (row) pi %.7lg",i,pi[i]);
02014           }
02015           break;
02016         case CPX_BASIC:
02017           break;
02018         default:
02019           TESTG((rval=1),CLEANUP,"Imposible! slack variable with upper bound!");
02020           break;
02021       }
02022     }
02023     /* update statistics */
02024     if(prevfixed == 0)
02025     {
02026       if(ncols > nfixed)
02027         ns_data->nzrd += ((size_t)(ncols - nfixed));
02028     }
02029     prevfixed += nfixed;
02030     if(prevfixed == ncols)
02031     {
02032       MESSAGE(MT_VERB_LVL-__MT_NS_LVL, "Nothing to do!, nextcol %d, prevfixed"
02033           " %d, ncols %d", nextcol, prevfixed, ncols);
02034       break;
02035     }
02036     TESTG((rval=(nextcol >= ncols+nrows)), CLEANUP, "Imposible! ncols %d "
02037         "nrows %d prevfixed %d loop %d nextcol %d", ncols, nrows, prevfixed, 
02038         loop, nextcol);
02039     MESSAGE(MT_VERB_LVL-__MT_NS_LVL, "Nextcol %d (%c), loop %d, to fix %d", 
02040         nextcol, (nextcol < ncols) ? ((cstat[nextcol]== CPX_AT_LOWER) ? 
02041           'L': 'U'): 'R', loop, ncols - prevfixed);
02042     /* now we know we must continue iterating, but we may not need to
02043      * re-optimize */
02044     if(nextcol < ncols)
02045     {
02046       /* if the next free variable was at its upper bound, we 'soft-fix it' */
02047       if(cstat[nextcol] == CPX_AT_UPPER)
02048       {
02049         ctmp = 'L';
02050         rval = CPXchgbds(env, lp, 1, &nextcol, &ctmp, ub+nextcol);
02051         MTcplexCHECKRVALG(env,rval,CLEANUP);
02052         cfix[nextcol] = 'U';
02053         nfixed=1;
02054         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Soft redo");
02055         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing %d to %.4lf",nextcol,ub[nextcol]);
02056         goto SOFT_REDO;
02057       }
02058       /* if it is at lower, but the space to upper bound is zero, we just fix
02059        * it */
02060       else if( fabs(ub[nextcol] - lb[nextcol]) < MT_ZERO_TOL)
02061       {
02062         cfix[nextcol] = 'L';
02063         nfixed = 1;
02064         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Soft redo");
02065         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing %d to %.4lf",nextcol,lb[nextcol]);
02066         goto SOFT_REDO;
02067       }
02068       /* otherwise we must optimize */
02069       else
02070       {
02071         memset(dj,0,sizeof(double)*((size_t)ncols));
02072         dj[nextcol] = 1.0*objsen;
02073         rval = CPXchgobj(env, lp, ncols, cind, dj);
02074         MTcplexCHECKRVALG(env,rval,CLEANUP);
02075         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Hard redo, prevfixed %d ncols %d", 
02076             prevfixed, ncols);
02077       }
02078     }
02079     /* if the next free is a slack variable ... we are done! */
02080     else 
02081     {
02082       nextcol -= ncols;
02083       /* soft-reoptimize */
02084       if(sense[nextcol] == 'E' || sense[nextcol] == 'e')
02085       {
02086         rfix[nextcol] = 'L';
02087         nfixed = 1;
02088         ctmp = 'E';
02089         rval = CPXchgsense(env, lp, 1, &nextcol, &ctmp);
02090         MTcplexCHECKRVALG(env,rval,CLEANUP);
02091         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Soft redo");
02092         MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Fixing row %d to equality",nextcol);
02093         goto SOFT_REDO;
02094       }
02095       /* hard re-optimization */
02096       rval = CPXgetrows(env, lp, &nz, &matbeg, matind, matval, ncols, &surplus, 
02097           nextcol, nextcol);
02098       MTcplexCHECKRVALG(env,rval,CLEANUP);
02099       if(sense[nextcol] == 'L') surplus = 1;
02100       else surplus = -1;
02101       memset(dj,0,sizeof(double)*((size_t)ncols));
02102       for( i = nz ; i-- ; ) dj[matind[i]] = matval[i]*objsen*surplus;
02103       rval = CPXchgobj(env, lp, ncols, cind, dj);
02104       MTcplexCHECKRVALG(env,rval,CLEANUP);
02105       MESSAGE(MT_VERB_LVL-__MT_NS_LVL,"Hard redo");
02106     }
02107   } while(1);
02108   /* now we copy the final basis to the original LP */
02109   /* get current basis (if any) in the ndoe LP */
02110   switch(lp_stat)
02111   {
02112     case CPX_STAT_OPTIMAL:
02113       dtmp = 0;
02114       for( i = ncols ; i-- ; ) dtmp += (x0[i]-x1[i])*(x0[i]-x1[i]);
02115       dtmp = sqrt(dtmp);
02116       ns_data->xdiff += dtmp;
02117       for( i = ncols ; i-- ; )
02118       {
02119         switch(cfix[i])
02120         {
02121           case 'F':
02122             cstat[i] = CPX_BASIC;
02123             break;
02124           case 'L':
02125             cstat[i] = CPX_AT_LOWER;
02126             break;
02127           case 'U':
02128             cstat[i] = CPX_AT_UPPER;
02129             break;
02130           case 'S':
02131             cstat[i] = CPX_FREE_SUPER;
02132             break;
02133           default:
02134             TESTG((rval=1), CLEANUP, "Imposible!");
02135             break;
02136         }
02137       }
02138       for( i = nrows ; i-- ; )
02139       {
02140         switch(rfix[i])
02141         {
02142           case 'L':
02143             rstat[i] = CPX_AT_LOWER;
02144             break;
02145           case 'F':
02146             rstat[i] = CPX_BASIC;
02147             break;
02148           default:
02149             TESTG((rval=1), CLEANUP, "Imposible!");
02150             break;
02151         }
02152       }
02153       rval = CPXcopybase(env, lp, cstat, rstat);
02154       MTcplexCHECKRVALG(env,rval,CLEANUP);
02155     case CPX_STAT_ABORT_OBJ_LIM:
02156     case CPX_STAT_UNBOUNDED:
02157     case CPX_STAT_INFEASIBLE:
02158       break;
02159     default:
02160       CPXgetstatstring(env, lp_stat, __mt_cplex_errbuf);
02161       TESTG((rval=1), CLEANUP, "Unknown LP status %d! LP0 stat %d, %s", 
02162           lp_stat,stat0, __mt_cplex_errbuf);
02163       break;
02164   }
02165   /* undo all bound, sense, and objective changes and objective change */
02166   rval = CPXchgobj(env, lp, ncols, cind, oobj);
02167   MTcplexCHECKRVALG(env,rval,CLEANUP);
02168   memset(cfix,'B',((size_t)ncols)*sizeof(char));
02169   rval = CPXchgbds(env, lp, ncols, cind, cfix, lb);
02170   MTcplexCHECKRVALG(env,rval,CLEANUP);
02171   memset(cfix,'U',((size_t)ncols)*sizeof(char));
02172   rval = CPXchgbds(env, lp, ncols, cind, cfix, ub);
02173   MTcplexCHECKRVALG(env,rval,CLEANUP);
02174   rval = CPXchgsense(env, lp, nrows, cind, sense);
02175   MTcplexCHECKRVALG(env,rval,CLEANUP);
02176   /* optimize with optimal basis */
02177   rval = CPXdualopt(env,lp);
02178   MTcplexCHECKRVALG(env,rval,CLEANUP);
02179   if(MT_DEBUG_LVL <= DEBUG && lp_stat == CPX_STAT_OPTIMAL)
02180   {
02181     rval = CPXgetx(env,lp,x0,0,ncols-1);
02182     MTcplexCHECKRVALG(env,rval,CLEANUP);
02183     for( i = 0 ; i<ncols ; i++ )
02184     {
02185       dtmp = fabs(x0[i]-x1[i]);
02186       TESTG((rval=(dtmp > MT_ZERO_TOL)), CLEANUP, "Imposible!, new solution is"
02187           " different! in column %d x_lex = %.5lf x_new = %.5lf, "
02188           "difference %.7le", i, x1[i], x0[i], dtmp);
02189     }
02190   }
02191   /* ending */
02192   rval = CPXsetdblparam(env, CPX_PARAM_OBJLLIM, llim);
02193   MTcplexCHECKRVALG(env,rval,CLEANUP);
02194   rval = CPXsetdblparam(env, CPX_PARAM_OBJULIM, ulim);
02195   MTcplexCHECKRVALG(env,rval,CLEANUP);
02196   CLEANUP:
02197   MESSAGE((rval ? 0 : MT_VERB_LVL-__MT_NS_LVL),"Done, code %d, loops %d, prevfixed %d, ncols %d"
02198       " LP0_stat %d", rval, loop, prevfixed, ncols, stat0 );
02199   ns_data->loops+=(size_t)loop;
02200   EGfree(cind);
02201   return rval;
02202 }
02203 /* ========================================================================= */
02204 /** @} */
02205 /* end of mt_cplex_cbk.c */
Generated on Mon Oct 26 09:16:29 2009 for MTgomory by  doxygen 1.4.6