exact.c

Go to the documentation of this file.
00001 /* ========================================================================= */
00002 /* ESolver "Exact Mixed Integer Linear Solver" provides some basic structures 
00003  * and algorithms commons in solving MIP's
00004  *
00005  * Copyright (C) 2005 Daniel Espinoza.
00006  * 
00007  * This library is free software; you can redistribute it and/or modify it
00008  * under the terms of the GNU Lesser General Public License as published by the
00009  * Free Software Foundation; either version 2.1 of the License, or (at your
00010  * option) any later version.
00011  *
00012  * This library is distributed in the hope that it will be useful, but 
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
00014  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
00015  * License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public License
00018  * along with this library; if not, write to the Free Software Foundation,
00019  * Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA 
00020  * */
00021 /* ========================================================================= */
00022 /** @file
00023  * @ingroup Esolver */
00024 /** @addtogroup Esolver */
00025 /** @{ */
00026 #include "exact.h"
00027 #include "util.h"
00028 extern mpq_t mpq_ILL_MAXDOUBLE;
00029 extern mpq_t mpq_ILL_MINDOUBLE;
00030 extern void mpq_ILLfct_set_variable_type (mpq_lpinfo * lp);
00031 extern void mpq_ILLfct_compute_piz (mpq_lpinfo * lp);
00032 extern void mpq_ILLfct_compute_dz (mpq_lpinfo * lp);
00033 extern void mpq_ILLfct_compute_xbz (mpq_lpinfo * lp);
00034 extern void mpq_ILLfct_compute_dobj ( mpq_lpinfo * lp);
00035 extern void mpq_ILLfct_check_pfeasible (mpq_lpinfo * lp,
00036                                         mpq_feas_info * fs,
00037                                         const mpq_t ftol);
00038 extern void mpq_ILLfct_check_dfeasible (mpq_lpinfo * lp,
00039                                         mpq_feas_info * fs,
00040                                         const mpq_t ftol);
00041 extern void mpq_ILLfct_set_status_values (mpq_lpinfo * lp,
00042                                           int pstatus,
00043                                           int dstatus,
00044                                           int ptype,
00045                                           int dtype);
00046 extern int mpq_grab_cache (mpq_QSdata * p,
00047                            int status);
00048 extern void mpq_ILLfct_compute_phaseI_piz (mpq_lpinfo * lp);
00049 
00050 /* ========================================================================= */
00051 int QSexact_print_sol (mpq_QSdata * p,
00052                        EGioFile_t * out_f)
00053 {
00054   int rval = 0,
00055     status;
00056   const int ncols = mpq_QSget_colcount (p);
00057   const int nrows = mpq_QSget_rowcount (p);
00058   mpq_t *x = mpq_EGlpNumAllocArray (ncols);
00059   mpq_t *rc = mpq_EGlpNumAllocArray (ncols);
00060   mpq_t *slack = mpq_EGlpNumAllocArray (nrows);
00061   mpq_t *pi = mpq_EGlpNumAllocArray (nrows);
00062   mpq_t value;
00063   register int i;
00064   char *str1 = 0;
00065   mpq_init (value);
00066   EGcallD(mpq_QSget_status (p, &status));
00067   if(mpq_QSget_x_array (p, x)) mpq_EGlpNumFreeArray (x);
00068   if(mpq_QSget_slack_array (p, slack)) mpq_EGlpNumFreeArray (slack);
00069   if(mpq_QSget_pi_array (p, pi)) mpq_EGlpNumFreeArray (pi);
00070   if(mpq_QSget_rc_array (p, rc)) mpq_EGlpNumFreeArray (rc);
00071 
00072   switch (status)
00073   {
00074   case QS_LP_OPTIMAL:
00075     EGcallD(mpq_QSget_objval (p, &value));
00076     str1 = mpq_EGlpNumGetStr (value);
00077     EGioPrintf (out_f, "status OPTIMAL\n\tValue = %s\n", str1);
00078     free (str1);
00079     str1 = 0;
00080     break;
00081   case QS_LP_INFEASIBLE:
00082     EGioPrintf (out_f, "status INFEASIBLE\n");
00083     break;
00084   case QS_LP_UNBOUNDED:
00085     EGioPrintf (out_f, "status UNBOUNDED\n");
00086     break;
00087   case QS_LP_ITER_LIMIT:
00088   case QS_LP_TIME_LIMIT:
00089   case QS_LP_UNSOLVED:
00090   case QS_LP_ABORTED:
00091   case QS_LP_MODIFIED:
00092     EGioPrintf (out_f, "status NOT_SOLVED\n");
00093     break;
00094   }
00095   if (x)
00096   {
00097     EGioPrintf (out_f, "VARS:\n");
00098     for (i = 0; i < ncols; i++)
00099       if (!mpq_equal (x[i], __zeroLpNum_mpq__))
00100       {
00101         str1 = mpq_EGlpNumGetStr (x[i]);
00102         EGioPrintf (out_f, "%s = %s\n", p->qslp->colnames[i], str1);
00103         free (str1);
00104       }
00105   }
00106   if (rc)
00107   {
00108     EGioPrintf (out_f, "REDUCED COST:\n");
00109     for (i = 0; i < ncols; i++)
00110       if (!mpq_equal (rc[i], __zeroLpNum_mpq__))
00111       {
00112         str1 = mpq_EGlpNumGetStr (rc[i]);
00113         EGioPrintf (out_f, "%s = %s\n", p->qslp->colnames[i], str1);
00114         free (str1);
00115       }
00116   }
00117   if (pi)
00118   {
00119     EGioPrintf (out_f, "PI:\n");
00120     for (i = 0; i < nrows; i++)
00121       if (!mpq_equal (pi[i], __zeroLpNum_mpq__))
00122       {
00123         str1 = mpq_EGlpNumGetStr (pi[i]);
00124         EGioPrintf (out_f, "%s = %s\n", p->qslp->rownames[i], str1);
00125         free (str1);
00126       }
00127   }
00128   if (slack)
00129   {
00130     EGioPrintf (out_f, "SLACK:\n");
00131     for (i = 0; i < nrows; i++)
00132       if (!mpq_equal (slack[i], __zeroLpNum_mpq__))
00133       {
00134         str1 = mpq_EGlpNumGetStr (slack[i]);
00135         EGioPrintf (out_f, "%s = %s\n", p->qslp->rownames[i], str1);
00136         free (str1);
00137       }
00138   }
00139 
00140   /* ending */
00141 CLEANUP:
00142   if (x)
00143     mpq_EGlpNumFreeArray (x);
00144   if (pi)
00145     mpq_EGlpNumFreeArray (pi);
00146   if (rc)
00147     mpq_EGlpNumFreeArray (rc);
00148   if (slack)
00149     mpq_EGlpNumFreeArray (slack);
00150   mpq_clear (value);
00151   return rval;
00152 }
00153 
00154 /* ========================================================================= */
00155 dbl_QSdata *QScopy_prob_mpq_dbl (mpq_QSdata * p,
00156                                  const char *newname)
00157 {
00158   const int ncol = mpq_QSget_colcount(p);
00159   const int nrow = mpq_QSget_rowcount(p);
00160   char*sense=0;
00161   int*rowcnt=0;
00162   int*rowbeg=0;
00163   int*rowind=0;
00164   int objsense;
00165   mpq_t*mpq_lb=0;
00166   mpq_t*mpq_ub=0;
00167   mpq_t*mpq_obj=0;
00168   mpq_t*mpq_range=0;
00169   mpq_t*mpq_rhs=0;
00170   mpq_t*mpq_rowval=0;
00171   double*dbl_lb=0;
00172   double*dbl_ub=0;
00173   double*dbl_obj=0;
00174   double*dbl_range=0;
00175   double*dbl_rhs=0;
00176   double*dbl_rowval=0;
00177   dbl_QSdata *p2 = 0;
00178   int rval = 0;
00179   register int i;
00180   mpq_t mpq_val;
00181   double dbl_val;
00182   mpq_init(mpq_val);
00183   /* get all information */
00184   EGcallD(mpq_QSget_objsense(p,&objsense));
00185   mpq_lb = mpq_EGlpNumAllocArray(ncol);
00186   mpq_ub = mpq_EGlpNumAllocArray(ncol);
00187   EGcallD(mpq_QSget_bounds(p,mpq_lb,mpq_ub));
00188   dbl_lb = QScopy_array_mpq_dbl(mpq_lb);
00189   dbl_ub = QScopy_array_mpq_dbl(mpq_ub);
00190   mpq_EGlpNumFreeArray(mpq_ub);
00191   mpq_obj = mpq_lb;
00192   mpq_lb = 0;
00193   EGcallD(mpq_QSget_obj(p, mpq_obj));
00194   dbl_obj = QScopy_array_mpq_dbl(mpq_obj);
00195   mpq_EGlpNumFreeArray(mpq_obj);
00196   EGcallD(mpq_QSget_ranged_rows(p, &rowcnt, &rowbeg, &rowind, &mpq_rowval,
00197                                 &mpq_rhs, &sense, &mpq_range, 0));
00198   dbl_rowval = QScopy_array_mpq_dbl(mpq_rowval);
00199   mpq_EGlpNumFreeArray(mpq_rowval);
00200   dbl_range = QScopy_array_mpq_dbl(mpq_range);
00201   mpq_EGlpNumFreeArray(mpq_range);
00202   dbl_rhs = QScopy_array_mpq_dbl(mpq_rhs);
00203   mpq_EGlpNumFreeArray(mpq_rhs);
00204   /* create copy */
00205   p2 = dbl_QScreate_prob (newname, objsense);
00206   if (!p2) goto CLEANUP;
00207   for( i = 0 ; i < ncol; i++)
00208   {
00209     EGcallD(dbl_QSnew_col(p2, dbl_obj[i], dbl_lb[i], dbl_ub[i], 0));
00210   }
00211   dbl_EGlpNumFreeArray(dbl_lb);
00212   dbl_EGlpNumFreeArray(dbl_ub);
00213   dbl_EGlpNumFreeArray(dbl_obj);
00214   EGcallD(dbl_QSadd_ranged_rows(p2, nrow, rowcnt, rowbeg, rowind, 
00215                                 dbl_rowval, dbl_rhs, sense, dbl_range, 0));
00216   /* set parameters */
00217   EGcallD(mpq_QSget_param(p, QS_PARAM_PRIMAL_PRICING, &objsense));
00218   EGcallD(dbl_QSset_param(p2, QS_PARAM_PRIMAL_PRICING, objsense));
00219   EGcallD(mpq_QSget_param(p, QS_PARAM_DUAL_PRICING, &objsense));
00220   EGcallD(dbl_QSset_param(p2, QS_PARAM_DUAL_PRICING, objsense));
00221   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_DISPLAY, &objsense));
00222   EGcallD(dbl_QSset_param(p2, QS_PARAM_SIMPLEX_DISPLAY, objsense));
00223   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_MAX_ITERATIONS, &objsense));
00224   EGcallD(dbl_QSset_param(p2, QS_PARAM_SIMPLEX_MAX_ITERATIONS, objsense));
00225   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_SCALING, &objsense));
00226   EGcallD(dbl_QSset_param(p2, QS_PARAM_SIMPLEX_SCALING, objsense));
00227   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_SIMPLEX_MAX_TIME, &mpq_val));
00228   dbl_val = mpq_get_d(mpq_val);
00229   EGcallD(dbl_QSset_param_EGlpNum(p2, QS_PARAM_SIMPLEX_MAX_TIME, dbl_val));
00230   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_OBJULIM, &mpq_val));
00231   dbl_val = mpq_get_d(mpq_val);
00232   EGcallD(dbl_QSset_param_EGlpNum(p2, QS_PARAM_OBJULIM, dbl_val));
00233   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_OBJLLIM, &mpq_val));
00234   dbl_val = mpq_get_d(mpq_val);
00235   EGcallD(dbl_QSset_param_EGlpNum(p2, QS_PARAM_OBJLLIM, dbl_val));
00236   /* ending */
00237   CLEANUP:
00238   mpq_clear(mpq_val);
00239   dbl_EGlpNumFreeArray(dbl_rowval);
00240   dbl_EGlpNumFreeArray(dbl_range);
00241   dbl_EGlpNumFreeArray(dbl_rhs);
00242   dbl_EGlpNumFreeArray(dbl_lb);
00243   dbl_EGlpNumFreeArray(dbl_ub);
00244   dbl_EGlpNumFreeArray(dbl_obj);
00245   mpq_EGlpNumFreeArray(mpq_rowval);
00246   mpq_EGlpNumFreeArray(mpq_range);
00247   mpq_EGlpNumFreeArray(mpq_rhs);
00248   mpq_EGlpNumFreeArray(mpq_lb);
00249   mpq_EGlpNumFreeArray(mpq_ub);
00250   mpq_EGlpNumFreeArray(mpq_obj);
00251   EGfree(rowcnt);
00252   EGfree(rowbeg);
00253   EGfree(rowind);
00254   EGfree(sense);
00255   if (rval)
00256   {
00257     dbl_QSfree_prob (p2);
00258     p2 = 0;
00259   }
00260 #if QSEXACT_SAVE_INT
00261   else
00262   {
00263     dbl_QSwrite_prob (p2, "prob.dbl.lp", "LP");
00264   }
00265 #endif
00266   return p2;
00267 }
00268 
00269 /* ========================================================================= */
00270 mpf_QSdata *QScopy_prob_mpq_mpf (mpq_QSdata * p,
00271                                  const char *newname)
00272 {
00273   const int ncol = mpq_QSget_colcount(p);
00274   const int nrow = mpq_QSget_rowcount(p);
00275   char*sense=0;
00276   int*rowcnt=0;
00277   int*rowbeg=0;
00278   int*rowind=0;
00279   int objsense;
00280   mpq_t*mpq_lb=0;
00281   mpq_t*mpq_ub=0;
00282   mpq_t*mpq_obj=0;
00283   mpq_t*mpq_range=0;
00284   mpq_t*mpq_rhs=0;
00285   mpq_t*mpq_rowval=0;
00286   mpf_t*mpf_lb=0;
00287   mpf_t*mpf_ub=0;
00288   mpf_t*mpf_obj=0;
00289   mpf_t*mpf_range=0;
00290   mpf_t*mpf_rhs=0;
00291   mpf_t*mpf_rowval=0;
00292   mpf_QSdata *p2 = 0;
00293   int rval = 0;
00294   mpq_t mpq_val;
00295   mpf_t mpf_val;
00296   register int i;
00297   mpq_init(mpq_val);
00298   mpf_init(mpf_val);
00299   /* get all information */
00300   EGcallD(mpq_QSget_objsense(p,&objsense));
00301   mpq_lb = mpq_EGlpNumAllocArray(ncol);
00302   mpq_ub = mpq_EGlpNumAllocArray(ncol);
00303   EGcallD(mpq_QSget_bounds(p,mpq_lb,mpq_ub));
00304   mpf_lb = QScopy_array_mpq_mpf(mpq_lb);
00305   mpf_ub = QScopy_array_mpq_mpf(mpq_ub);
00306   mpq_EGlpNumFreeArray(mpq_ub);
00307   mpq_obj = mpq_lb;
00308   mpq_lb = 0;
00309   EGcallD(mpq_QSget_obj(p, mpq_obj));
00310   mpf_obj = QScopy_array_mpq_mpf(mpq_obj);
00311   mpq_EGlpNumFreeArray(mpq_obj);
00312   EGcallD(mpq_QSget_ranged_rows(p, &rowcnt, &rowbeg, &rowind, &mpq_rowval, &mpq_rhs, &sense, &mpq_range, 0));
00313   mpf_rowval = QScopy_array_mpq_mpf(mpq_rowval);
00314   mpq_EGlpNumFreeArray(mpq_rowval);
00315   mpf_range = QScopy_array_mpq_mpf(mpq_range);
00316   mpq_EGlpNumFreeArray(mpq_range);
00317   mpf_rhs = QScopy_array_mpq_mpf(mpq_rhs);
00318   mpq_EGlpNumFreeArray(mpq_rhs);
00319   /* create copy */
00320   p2 = mpf_QScreate_prob (newname, objsense);
00321   if (!p2) goto CLEANUP;
00322   for( i = 0 ; i < ncol; i++)
00323   {
00324     EGcallD(mpf_QSnew_col(p2, mpf_obj[i], mpf_lb[i], mpf_ub[i], 0));
00325   }
00326   mpf_EGlpNumFreeArray(mpf_lb);
00327   mpf_EGlpNumFreeArray(mpf_ub);
00328   mpf_EGlpNumFreeArray(mpf_obj);
00329   EGcallD(mpf_QSadd_ranged_rows(p2, nrow, rowcnt, rowbeg, rowind, (const mpf_t*)mpf_rowval,(const mpf_t*) mpf_rhs, sense,(const mpf_t*) mpf_range, 0));
00330   /* set parameters */
00331   EGcallD(mpq_QSget_param(p, QS_PARAM_PRIMAL_PRICING, &objsense));
00332   EGcallD(mpf_QSset_param(p2, QS_PARAM_PRIMAL_PRICING, objsense));
00333   EGcallD(mpq_QSget_param(p, QS_PARAM_DUAL_PRICING, &objsense));
00334   EGcallD(mpf_QSset_param(p2, QS_PARAM_DUAL_PRICING, objsense));
00335   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_DISPLAY, &objsense));
00336   EGcallD(mpf_QSset_param(p2, QS_PARAM_SIMPLEX_DISPLAY, objsense));
00337   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_MAX_ITERATIONS, &objsense));
00338   EGcallD(mpf_QSset_param(p2, QS_PARAM_SIMPLEX_MAX_ITERATIONS, objsense));
00339   EGcallD(mpq_QSget_param(p, QS_PARAM_SIMPLEX_SCALING, &objsense));
00340   EGcallD(mpf_QSset_param(p2, QS_PARAM_SIMPLEX_SCALING, objsense));
00341   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_SIMPLEX_MAX_TIME, &mpq_val));
00342   mpf_set_q(mpf_val,mpq_val);
00343   EGcallD(mpf_QSset_param_EGlpNum(p2, QS_PARAM_SIMPLEX_MAX_TIME, mpf_val));
00344   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_OBJULIM, &mpq_val));
00345   mpf_set_q(mpf_val,mpq_val);
00346   EGcallD(mpf_QSset_param_EGlpNum(p2, QS_PARAM_OBJULIM, mpf_val));
00347   EGcallD(mpq_QSget_param_EGlpNum(p, QS_PARAM_OBJLLIM, &mpq_val));
00348   mpf_set_q(mpf_val,mpq_val);
00349   EGcallD(mpf_QSset_param_EGlpNum(p2, QS_PARAM_OBJLLIM, mpf_val));
00350   /* ending */
00351   CLEANUP:
00352   mpq_clear(mpq_val);
00353   mpf_clear(mpf_val);
00354   mpf_EGlpNumFreeArray(mpf_rowval);
00355   mpf_EGlpNumFreeArray(mpf_range);
00356   mpf_EGlpNumFreeArray(mpf_rhs);
00357   mpf_EGlpNumFreeArray(mpf_lb);
00358   mpf_EGlpNumFreeArray(mpf_ub);
00359   mpf_EGlpNumFreeArray(mpf_obj);
00360   mpq_EGlpNumFreeArray(mpq_rowval);
00361   mpq_EGlpNumFreeArray(mpq_range);
00362   mpq_EGlpNumFreeArray(mpq_rhs);
00363   mpq_EGlpNumFreeArray(mpq_lb);
00364   mpq_EGlpNumFreeArray(mpq_ub);
00365   mpq_EGlpNumFreeArray(mpq_obj);
00366   EGfree(rowcnt);
00367   EGfree(rowbeg);
00368   EGfree(rowind);
00369   EGfree(sense);
00370   if (rval)
00371   {
00372     mpf_QSfree_prob (p2);
00373     p2 = 0;
00374   }
00375 #if QSEXACT_SAVE_INT
00376   else
00377   {
00378     mpf_QSwrite_prob (p2, "prob.mpf.lp", "LP");
00379   }
00380 #endif
00381   return p2;
00382 }
00383 
00384 #if QSEXACT_SAVE_OPTIMAL
00385 /* ========================================================================= */
00386 /** @brief used to enumerate the generated optimal tests */
00387 static int QSEXACT_SAVE_OPTIMAL_IND = 0;
00388 #endif
00389 
00390 /* ========================================================================= */
00391 int QSexact_optimal_test (mpq_QSdata * p,
00392                           mpq_t * p_sol,
00393                           mpq_t * d_sol,
00394                           QSbasis * basis)
00395 {
00396   /* local variables */
00397   register int i,
00398     j;
00399   mpq_ILLlpdata *qslp = p->lp->O;
00400   int *iarr1 = 0,
00401    *rowmap = qslp->rowmap,
00402    *structmap = qslp->structmap,
00403     col;
00404   mpq_t *arr1 = 0,
00405    *arr2 = 0,
00406    *arr3 = 0,
00407    *arr4 = 0,
00408    *rhs_copy = 0;
00409   mpq_t *dz = 0;
00410   int objsense = (qslp->objsense == QS_MIN) ? 1 : -1;
00411   int const msg_lvl = __QS_SB_VERB <= DEBUG ? 0 : 100000 * (1 - p->simplex_display);
00412   int rval = 1;                 /* store whether or not the solution is optimal, we start 
00413                                  * assuming it is. */
00414   mpq_t num1,
00415     num2,
00416     num3,
00417     p_obj,
00418     d_obj;
00419   mpq_init (num1);
00420   mpq_init (num2);
00421   mpq_init (num3);
00422   mpq_init (p_obj);
00423   mpq_init (d_obj);
00424   mpq_set_ui (p_obj, 0UL, 1UL);
00425   mpq_set_ui (d_obj, 0UL, 1UL);
00426 
00427   /* now check if the given basis is the optimal basis */
00428   arr3 = qslp->lower;
00429   arr4 = qslp->upper;
00430   if (mpq_QSload_basis (p, basis))
00431   {
00432     rval = 0;
00433     MESSAGE (msg_lvl, "QSload_basis failed");
00434     goto CLEANUP;
00435   }
00436   for (i = basis->nstruct; i--;)
00437   {
00438     /* check that the upper and lower bound define a non-empty space */
00439     if (mpq_cmp (arr3[structmap[i]], arr4[structmap[i]]) > 0)
00440     {
00441       rval = 0;
00442       if(!msg_lvl)
00443       {
00444         MESSAGE(0, "variable %s has empty feasible range [%lg,%lg]",
00445                  qslp->colnames[i], mpq_EGlpNumToLf(arr3[structmap[i]]), 
00446                  mpq_EGlpNumToLf(arr4[structmap[i]]));
00447       }
00448       goto CLEANUP;
00449     }
00450     /* set the variable to its apropiate values, depending its status */
00451     switch (basis->cstat[i])
00452     {
00453     case QS_COL_BSTAT_FREE:
00454     case QS_COL_BSTAT_BASIC:
00455       if (mpq_cmp (p_sol[i], arr4[structmap[i]]) > 0)
00456         mpq_set (p_sol[i], arr4[structmap[i]]);
00457       else if (mpq_cmp (p_sol[i], arr3[structmap[i]]) < 0)
00458         mpq_set (p_sol[i], arr3[structmap[i]]);
00459       break;
00460     case QS_COL_BSTAT_UPPER:
00461       mpq_set (p_sol[i], arr4[structmap[i]]);
00462       break;
00463     case QS_COL_BSTAT_LOWER:
00464       mpq_set (p_sol[i], arr3[structmap[i]]);
00465       break;
00466     default:
00467       rval = 0;
00468       MESSAGE (msg_lvl, "Unknown Variable basic status %d, for variable "
00469                "(%s,%d)", basis->cstat[i], qslp->colnames[i], i);
00470       goto CLEANUP;
00471       break;
00472     }
00473   }
00474   for (i = basis->nrows; i--;)
00475   {
00476     /* check that the upper and lower bound define a non-empty space */
00477     if (mpq_cmp (arr3[rowmap[i]], arr4[rowmap[i]]) > 0)
00478     {
00479       rval = 0;
00480       if(!msg_lvl)
00481       {
00482         MESSAGE(0, "constraint %s logical has empty feasible range "
00483                  "[%lg,%lg]", qslp->rownames[i], 
00484                  mpq_EGlpNumToLf(arr3[rowmap[i]]), 
00485                  mpq_EGlpNumToLf(arr4[rowmap[i]]));
00486       }
00487       goto CLEANUP;
00488     }
00489     /* set the variable to its apropiate values, depending its status */
00490     switch (basis->rstat[i])
00491     {
00492     case QS_ROW_BSTAT_BASIC:
00493       if (mpq_cmp (p_sol[i + basis->nstruct], arr4[rowmap[i]]) > 0)
00494         mpq_set (p_sol[i + basis->nstruct], arr4[rowmap[i]]);
00495       else if (mpq_cmp (p_sol[i + basis->nstruct], arr3[rowmap[i]]) < 0)
00496         mpq_set (p_sol[i + basis->nstruct], arr3[rowmap[i]]);
00497       break;
00498     case QS_ROW_BSTAT_UPPER:
00499       mpq_set (p_sol[i + basis->nstruct], arr4[rowmap[i]]);
00500       break;
00501     case QS_ROW_BSTAT_LOWER:
00502       mpq_set (p_sol[i + basis->nstruct], arr3[rowmap[i]]);
00503       break;
00504     default:
00505       rval = 0;
00506       MESSAGE (msg_lvl, "Unknown Variable basic status %d, for constraint "
00507                "(%s,%d)", basis->cstat[i], qslp->rownames[i], i);
00508       goto CLEANUP;
00509       break;
00510     }
00511   }
00512 
00513   /* compute the actual RHS */
00514   rhs_copy = mpq_EGlpNumAllocArray (qslp->nrows);
00515   for (i = qslp->nstruct; i--;)
00516   {
00517     if (!mpq_equal (p_sol[i], mpq_zeroLpNum))
00518     {
00519       arr1 = qslp->A.matval + qslp->A.matbeg[structmap[i]];
00520       iarr1 = qslp->A.matind + qslp->A.matbeg[structmap[i]];
00521       for (j = qslp->A.matcnt[structmap[i]]; j--;)
00522       {
00523         mpq_mul (num1, arr1[j], p_sol[i]);
00524         mpq_add (rhs_copy[iarr1[j]], rhs_copy[iarr1[j]], num1);
00525       }
00526     }
00527   }
00528 
00529   /* now check if both rhs and copy_rhs are equal */
00530   arr4 = qslp->upper;
00531   arr1 = qslp->rhs;
00532   arr2 = qslp->lower;
00533   for (i = qslp->nrows; i--;)
00534   {
00535     mpq_mul (num1, arr1[i], d_sol[i]);
00536     mpq_add (d_obj, d_obj, num1);
00537     mpq_sub (num2, arr1[i], rhs_copy[i]);
00538     EXIT (qslp->A.matcnt[rowmap[i]] != 1, "Imposible!");
00539     if (basis->rstat[i] == QS_ROW_BSTAT_BASIC)
00540       mpq_div (p_sol[qslp->nstruct + i], num2,
00541                qslp->A.matval[qslp->A.matbeg[rowmap[i]]]);
00542     else
00543     {
00544       mpq_mul (num1, p_sol[qslp->nstruct + i],
00545                qslp->A.matval[qslp->A.matbeg[rowmap[i]]]);
00546       if (!mpq_equal (num1, num2))
00547       {
00548         rval = 0;
00549         if(!msg_lvl)
00550         {
00551           MESSAGE(0, "solution is infeasible for constraint %s, violation"
00552                  " %lg", qslp->rownames[i],
00553                  mpq_get_d (num1) - mpq_get_d (num2));
00554         }
00555         goto CLEANUP;
00556       }
00557     }
00558     mpq_set (num2, p_sol[qslp->nstruct + i]);
00559     /* now we check the bounds on the logical variables */
00560     if (mpq_cmp (num2, arr2[rowmap[i]]) < 0)
00561     {
00562       rval = 0;
00563       if(!msg_lvl)
00564       {
00565         MESSAGE(0, "constraint %s artificial (%lg) bellow lower"
00566                  " bound (%lg), actual LHS (%lg), actual RHS (%lg)",
00567                  qslp->rownames[i], mpq_get_d (num2), 
00568                  mpq_get_d (arr2[rowmap[i]]), mpq_get_d (rhs_copy[i]), 
00569                  mpq_get_d (arr1[i]));
00570       }
00571       goto CLEANUP;
00572     }
00573     else if (mpq_cmp (num2, arr4[rowmap[i]]) > 0)
00574     {
00575       rval = 0;
00576       if(!msg_lvl)
00577       {
00578         MESSAGE(0, "constraint %s artificial (%lg) bellow lower bound"
00579                  " (%lg)", qslp->rownames[i], mpq_get_d (num2),
00580                  mpq_get_d (arr4[rowmap[i]]));
00581       }
00582       goto CLEANUP;
00583     }
00584   }
00585 
00586   /* compute the upper and lower bound dual variables, note that dl is the dual
00587    * of the lower bounds, and du the dual of the upper bound, dl >= 0 and du <=
00588    * 0 and A^t y + Idl + Idu = c, and the dual objective value is 
00589    * max y*b + l*dl + u*du, we colapse both vector dl and du into dz, note that
00590    * if we are maximizing, then dl <= 0 and du >=0 */
00591   dz = mpq_EGlpNumAllocArray (qslp->ncols);
00592   arr2 = qslp->obj;
00593   arr3 = qslp->lower;
00594   arr4 = qslp->upper;
00595   for (i = qslp->nstruct; i--;)
00596   {
00597     col = structmap[i];
00598     mpq_mul (num1, arr2[col], p_sol[i]);
00599     mpq_add (p_obj, p_obj, num1);
00600     arr1 = qslp->A.matval + qslp->A.matbeg[col];
00601     iarr1 = qslp->A.matind + qslp->A.matbeg[col];
00602     mpq_set (num1, arr2[col]);
00603     for (j = qslp->A.matcnt[col]; j--;)
00604     {
00605       mpq_mul (num2, arr1[j], d_sol[iarr1[j]]);
00606       mpq_sub (num1, num1, num2);
00607     }
00608     mpq_set (dz[col], num1);
00609     /* objective update */
00610     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) > 0)
00611     {
00612       mpq_mul (num3, dz[col], arr3[col]);
00613       mpq_add (d_obj, d_obj, num3);
00614     }
00615     else
00616     {
00617       mpq_mul (num3, dz[col], arr4[col]);
00618       mpq_add (d_obj, d_obj, num3);
00619     }
00620     /* now we check that only when the logical is tight then the dual
00621      * variable may be non-zero, also check for primal feasibility with respect
00622      * to lower/upper bounds. */
00623     mpq_set_ui (num2, 0UL, 1UL);
00624     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) > 0)
00625     {
00626       mpq_sub (num1, p_sol[i], arr3[col]);
00627       mpq_mul (num2, num1, dz[col]);
00628     }
00629     if (mpq_cmp_ui (num2, 0UL, 1UL) != 0)
00630     {
00631       rval = 0;
00632       if(!msg_lvl)
00633       {
00634         MESSAGE(0, "lower bound (%s,%d) slack (%lg) and dual variable (%lg)"
00635                  " don't satisfy complementary slacknes %s",
00636                  qslp->colnames[i], i, mpq_get_d(num1), mpq_get_d(dz[col]),
00637                  "(real)");
00638       }
00639       goto CLEANUP;
00640     }
00641     mpq_set_ui (num2, 0UL, 1UL);
00642     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) < 0)
00643     {
00644       mpq_sub (num1, p_sol[i], arr4[col]);
00645       mpq_mul (num2, num1, dz[col]);
00646     }
00647     if (mpq_cmp_ui (num2, 0UL, 1UL) != 0)
00648     {
00649       rval = 0;
00650       if(!msg_lvl)
00651       {
00652         MESSAGE(0, "upper bound (%lg) variable (%lg) and dual variable"
00653                   " (%lg) don't satisfy complementary slacknes for variable"
00654                   " (%s,%d) %s", mpq_get_d(arr4[col]), mpq_get_d(p_sol[i]),
00655                   mpq_get_d(dz[col]), qslp->colnames[i], i, "(real)");
00656       }
00657       goto CLEANUP;
00658     }
00659   }
00660   /* complenetary slackness checked, now update the same for the logical
00661    * variables */
00662   for (i = qslp->nrows; i--;)
00663   {
00664     col = rowmap[i];
00665     mpq_mul (num1, arr2[col], p_sol[i + qslp->nstruct]);
00666     WARNING (mpq_cmp (arr2[col], mpq_zeroLpNum), "logical variable %s with "
00667              "non-zero objective function %lf", qslp->rownames[i],
00668              mpq_get_d (arr2[col]));
00669     mpq_add (p_obj, p_obj, num1);
00670     arr1 = qslp->A.matval + qslp->A.matbeg[col];
00671     iarr1 = qslp->A.matind + qslp->A.matbeg[col];
00672     mpq_set (num1, arr2[col]);
00673     for (j = qslp->A.matcnt[col]; j--;)
00674     {
00675       mpq_mul (num2, arr1[j], d_sol[iarr1[j]]);
00676       mpq_sub (num1, num1, num2);
00677     }
00678     mpq_set (dz[col], num1);
00679     /* objective update */
00680     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) > 0)
00681     {
00682       mpq_mul (num3, dz[col], arr3[col]);
00683       mpq_add (d_obj, d_obj, num3);
00684     }
00685     else
00686     {
00687       mpq_mul (num3, dz[col], arr4[col]);
00688       mpq_add (d_obj, d_obj, num3);
00689     }
00690     /* now we check that only when the primal variable is tight then the dual
00691      * variable may be non-zero, also check for primal feasibility with respect
00692      * to lower/upper bounds. */
00693     mpq_set_ui (num2, 0UL, 1UL);
00694     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) > 0)
00695     {
00696       mpq_sub (num1, p_sol[i + qslp->nstruct], arr3[col]);
00697       mpq_mul (num2, num1, dz[col]);
00698     }
00699     if (mpq_cmp_ui (num2, 0UL, 1UL) != 0)
00700     {
00701       rval = 0;
00702       if(!msg_lvl)
00703       {
00704         MESSAGE(0, "lower bound (%s,%d) slack (%lg) and dual variable (%lg)"
00705                  " don't satisfy complementary slacknes %s", 
00706                  qslp->colnames[col], i, mpq_get_d(num1), mpq_get_d(dz[col]), 
00707                  "(real)");
00708       }
00709       goto CLEANUP;
00710     }
00711     mpq_set_ui (num2, 0UL, 1UL);
00712     if (objsense * mpq_cmp_ui (dz[col], 0UL, 1UL) < 0)
00713     {
00714       mpq_sub (num1, p_sol[i + qslp->nstruct], arr4[col]);
00715       mpq_mul (num2, num1, dz[col]);
00716     }
00717     if (mpq_cmp_ui (num2, 0UL, 1UL) != 0)
00718     {
00719       rval = 0;
00720       if(!msg_lvl)
00721       {
00722         MESSAGE(0, "upper bound (%lg) variable (%lg) and dual variable"
00723                 " (%lg) don't satisfy complementary slacknes for variable "
00724                 "(%s,%d) %s", mpq_get_d(arr4[col]), 
00725                 mpq_get_d(p_sol[i+qslp->nstruct]), mpq_get_d(dz[col]), qslp->colnames[col], i,
00726                 "(real)");
00727       }
00728       goto CLEANUP;
00729     }
00730   }
00731 
00732   /* now check the objective values */
00733   if (mpq_cmp (p_obj, d_obj) != 0)
00734   {
00735     rval = 0;
00736     if(!msg_lvl)
00737     {
00738       MESSAGE(0, "primal and dual objective value differ %lg %lg",
00739                mpq_get_d(p_obj), mpq_get_d(d_obj));
00740     }
00741     goto CLEANUP;
00742   }
00743   /* now we report optimality */
00744   if(!msg_lvl)
00745   {
00746     MESSAGE(0, "Problem solved to optimality, LP value %lg", mpq_get_d(p_obj));
00747   }
00748   /* now we load into cache the solution */
00749   if (!p->cache)
00750   {
00751     p->cache = EGsMalloc (mpq_ILLlp_cache, 1);
00752     mpq_EGlpNumInitVar (p->cache->val);
00753     mpq_ILLlp_cache_init (p->cache);
00754   }
00755   if (qslp->nrows != p->cache->nrows || qslp->nstruct != p->cache->nstruct)
00756   {
00757     mpq_ILLlp_cache_free (p->cache);
00758     EGcallD(mpq_ILLlp_cache_alloc (p->cache, qslp->nstruct, qslp->nrows));
00759   }
00760   p->cache->status = QS_LP_OPTIMAL;
00761   p->qstatus = QS_LP_OPTIMAL;
00762   p->lp->basisstat.optimal = 1;
00763   mpq_set (p->cache->val, p_obj);
00764   for (i = qslp->nstruct; i--;)
00765   {
00766     mpq_set (p->cache->x[i], p_sol[i]);
00767     mpq_set (p->cache->rc[i], dz[structmap[i]]);
00768   }
00769   for (i = qslp->nrows; i--;)
00770   {
00771     mpq_set (p->cache->slack[i], p_sol[i + qslp->nstruct]);
00772     mpq_set (p->cache->pi[i], d_sol[i]);
00773   }
00774 
00775   /* save the problem and solution if enablred */
00776 #if QSEXACT_SAVE_OPTIMAL
00777   {
00778     char stmp[1024];
00779     EGioFile_t *out_f = 0;
00780     snprintf (stmp, 1023, "%s-opt%03d.lp", p->name ? p->name : "UNNAMED",
00781               QSEXACT_SAVE_OPTIMAL_IND);
00782     if (mpq_QSwrite_prob (p, stmp, "LP"))
00783     {
00784       rval = 0;
00785       MESSAGE (0, "Couldn't write output problem %s", stmp);
00786       goto CLEANUP;
00787     }
00788     snprintf (stmp, 1023, "%s-opt%03d.sol.gz", p->name ? p->name : "UNNAMED",
00789               QSEXACT_SAVE_OPTIMAL_IND);
00790     if (!(out_f = EGioOpen (stmp, "w+")))
00791     {
00792       rval = 0;
00793       MESSAGE (0, "Couldn't open solution file %s", stmp);
00794       goto CLEANUP;
00795     }
00796     if (QSexact_print_sol (p, out_f))
00797     {
00798       rval = 0;
00799       MESSAGE (0, "Couldn't write output solution %s", stmp);
00800       goto CLEANUP;
00801     }
00802     EGioClose (out_f);
00803     QSEXACT_SAVE_OPTIMAL_IND++;
00804   }
00805 #endif
00806   rval = 1;
00807 
00808   /* ending */
00809 CLEANUP:
00810   mpq_EGlpNumFreeArray (dz);
00811   mpq_EGlpNumFreeArray (rhs_copy);
00812   mpq_clear (num1);
00813   mpq_clear (num2);
00814   mpq_clear (num3);
00815   mpq_clear (p_obj);
00816   mpq_clear (d_obj);
00817   return rval;
00818 }
00819 
00820 /* ========================================================================= */
00821 int QSexact_infeasible_test (mpq_QSdata * p,
00822                              mpq_t * d_sol)
00823 {
00824   /* local variables */
00825   register int i,
00826     j;
00827   int *iarr1;
00828   mpq_ILLlpdata *qslp = p->lp->O;
00829   mpq_t *arr1,
00830    *arr2,
00831    *arr3,
00832    *arr4;
00833   mpq_t *dl = 0,
00834    *du = 0;
00835   int const msg_lvl = __QS_SB_VERB <= DEBUG ? 0 : 100000 * (1 - p->simplex_display);
00836   int rval = 1;                 /* store whether or not the solution is optimal, we start 
00837                                  * assuming it is. */
00838   mpq_t num1,
00839     num2,
00840     num3,
00841     d_obj;
00842   mpq_init (num1);
00843   mpq_init (num2);
00844   mpq_init (num3);
00845   mpq_init (d_obj);
00846   mpq_set_ui (d_obj, 0UL, 1UL);
00847 
00848   /* compute the dual objective value */
00849   arr2 = qslp->rhs;
00850   for (i = qslp->nrows; i--;)
00851   {
00852     mpq_mul (num1, arr2[i], d_sol[i]);
00853     mpq_add (d_obj, d_obj, num1);
00854   }
00855 
00856   /* compute the upper and lower bound dual variables, note that dl is the dual
00857    * of the lower bounds, and du the dual of the upper bound, dl <= 0 and du >=
00858    * 0 and A^t y + Idl + Idu = c, and the dual objective value is 
00859    * max y*b + l*dl + u*du */
00860   du = mpq_EGlpNumAllocArray (qslp->ncols);
00861   dl = mpq_EGlpNumAllocArray (qslp->ncols);
00862   arr3 = qslp->lower;
00863   arr4 = qslp->upper;
00864   for (i = qslp->ncols; i--;)
00865   {
00866     arr1 = qslp->A.matval + qslp->A.matbeg[i];
00867     iarr1 = qslp->A.matind + qslp->A.matbeg[i];
00868     mpq_set_ui (num1, 0UL, 1UL);
00869     mpq_set_ui (du[i], 0UL, 1UL);
00870     mpq_set_ui (dl[i], 0UL, 1UL);
00871     for (j = qslp->A.matcnt[i]; j--;)
00872     {
00873       mpq_mul (num2, arr1[j], d_sol[iarr1[j]]);
00874       mpq_sub (num1, num1, num2);
00875     }
00876     if (mpq_cmp_ui (num1, 0UL, 1UL) < 0)
00877       mpq_set (du[i], num1);
00878     else
00879       mpq_set (dl[i], num1);
00880     if (mpq_equal (arr4[i], mpq_ILL_MAXDOUBLE) &&
00881         mpq_cmp_ui (du[i], 0UL, 1UL) != 0)
00882     {
00883       rval = 0;
00884       if(!msg_lvl)
00885       {
00886         MESSAGE(0, "upper bound of variable is INFTY, and it's dual is "
00887                  "non-zero %lg", mpq_get_d(du[i]));
00888       }
00889       goto CLEANUP;
00890     }
00891     if (mpq_equal (arr3[i], mpq_ILL_MINDOUBLE) &&
00892         mpq_cmp_ui (dl[i], 0UL, 1UL) != 0)
00893     {
00894       rval = 0;
00895       if(!msg_lvl)
00896       {
00897         MESSAGE(0, "lower bound of variable is -INFTY, and it's dual is "
00898                  "non-zero %lg", mpq_get_d(dl[i]));
00899       }
00900       goto CLEANUP;
00901     }
00902     mpq_mul (num3, dl[i], arr3[i]);
00903     mpq_add (d_obj, d_obj, num3);
00904     mpq_mul (num3, du[i], arr4[i]);
00905     mpq_add (d_obj, d_obj, num3);
00906   }
00907   /* now check the objective values */
00908   if (mpq_cmp_ui (d_obj, 0UL, 1UL) <= 0)
00909   {
00910     rval = 0;
00911     if(!msg_lvl)
00912     {
00913       MESSAGE(0, "dual ray is feasible, but objective is non "
00914               "positive %lg", mpq_get_d(d_obj));
00915     }
00916     goto CLEANUP;
00917   }
00918   p->qstatus = QS_LP_INFEASIBLE;
00919 
00920   /* ending */
00921 CLEANUP:
00922   mpq_EGlpNumFreeArray (dl);
00923   mpq_EGlpNumFreeArray (du);
00924   mpq_clear (num1);
00925   mpq_clear (num2);
00926   mpq_clear (num3);
00927   mpq_clear (d_obj);
00928   return rval;
00929 }
00930 
00931 /* ========================================================================= */
00932 /** @brief Used as separator while printing output to the screen (controled by
00933  * enabling simplex_display in the mpq_QSdata */
00934 /* ========================================================================= */
00935 static const char __sp[81] =
00936   "================================================================================";
00937 
00938 /* ========================================================================= */
00939 /** @brief print into screen (if enable) a message indicating that we have
00940  * successfully prove infeasibility, and save (if y is non
00941  * NULL ) the dual ray solution provided in y_mpq.
00942  * @param p_mpq the problem data.
00943  * @param y where to store the optimal dual solution (if not null).
00944  * @param y_mpq  the optimal dual solution.
00945  * */
00946 /* ========================================================================= */
00947 static void infeasible_output (mpq_QSdata * p_mpq,
00948                                mpq_t * const y,
00949                                mpq_t * y_mpq)
00950 {
00951   if (p_mpq->simplex_display)
00952   {
00953     fprintf (stdout, "%s\n\tProblem Is Infeasible\n%s\n", __sp, __sp);
00954     fflush(stdout);
00955   }
00956   if (y)
00957   {
00958     unsigned sz = __EGlpNumArraySize (y_mpq);
00959     while (sz--)
00960       mpq_set (y[sz], y_mpq[sz]);
00961   }
00962 }
00963 
00964 /* ========================================================================= */
00965 /** @brief print into screen (if enable) a message indicating that we have
00966  * successfully solved the problem at optimality, and save (if x and y are non
00967  * NULL respectivelly) the optimal primal/dual solution provided in x_mpq and
00968  * y_mpq. 
00969  * @param p_mpq the problem data.
00970  * @param x where to store the optimal primal solution (if not null).
00971  * @param y where to store the optimal dual solution (if not null).
00972  * @param x_mpq  the optimal primal solution.
00973  * @param y_mpq  the optimal dual solution.
00974  * */
00975 /* ========================================================================= */
00976 static void optimal_output (mpq_QSdata * p_mpq,
00977                             mpq_t * const x,
00978                             mpq_t * const y,
00979                             mpq_t * x_mpq,
00980                             mpq_t * y_mpq)
00981 {
00982   if (p_mpq->simplex_display)
00983   {
00984     fprintf (stdout, "%s\n\tProblem Solved Exactly\n%s\n", __sp, __sp);
00985     fflush(stdout);
00986   }
00987   if (y)
00988   {
00989     unsigned sz = __EGlpNumArraySize (y_mpq);
00990     while (sz--)
00991       mpq_set (y[sz], y_mpq[sz]);
00992   }
00993   if (x)
00994   {
00995     unsigned sz = __EGlpNumArraySize (x_mpq);
00996     while (sz--)
00997       mpq_set (x[sz], x_mpq[sz]);
00998   }
00999 }
01000 
01001 /* ========================================================================= */
01002 /** @brief get the status for a given basis in rational arithmetic, it should
01003  * also leave everything set to get primal/dual solutions when needed.
01004  * */
01005 static int QSexact_basis_status (mpq_QSdata * p_mpq,
01006                                  int *status,
01007                                  QSbasis * const basis,
01008                                  const int msg_lvl,
01009                                  int *const simplexalgo)
01010 {
01011   int rval = 0,
01012   singular;
01013   mpq_feas_info fi;
01014   EGtimer_t local_timer;
01015   mpq_EGlpNumInitVar (fi.totinfeas);
01016   EGtimerReset (&local_timer);
01017   EGtimerStart (&local_timer);
01018   EGcallD(mpq_QSload_basis (p_mpq, basis));
01019   if (p_mpq->cache) 
01020   {
01021     mpq_ILLlp_cache_free (p_mpq->cache);
01022     mpq_clear (p_mpq->cache->val);
01023     ILL_IFFREE (p_mpq->cache, mpq_ILLlp_cache); 
01024   }
01025   p_mpq->qstatus = QS_LP_MODIFIED;
01026   if(p_mpq->qslp->sinfo) 
01027   {
01028     mpq_ILLlp_sinfo_free(p_mpq->qslp->sinfo);
01029     ILL_IFFREE(p_mpq->qslp->sinfo, mpq_ILLlp_sinfo); 
01030   }
01031   if(p_mpq->qslp->rA)
01032   {
01033     mpq_ILLlp_rows_clear (p_mpq->qslp->rA);
01034     ILL_IFFREE (p_mpq->qslp->rA, mpq_ILLlp_rows);
01035   }
01036   mpq_free_internal_lpinfo (p_mpq->lp);
01037   mpq_init_internal_lpinfo (p_mpq->lp);
01038   EGcallD(mpq_build_internal_lpinfo (p_mpq->lp));
01039   mpq_ILLfct_set_variable_type (p_mpq->lp);
01040   EGcallD(mpq_ILLbasis_load (p_mpq->lp, p_mpq->basis));
01041   EGcallD(mpq_ILLbasis_factor (p_mpq->lp, &singular));
01042   memset (&(p_mpq->lp->basisstat), 0, sizeof (mpq_lp_status_info));
01043   mpq_ILLfct_compute_piz (p_mpq->lp);
01044   mpq_ILLfct_compute_dz (p_mpq->lp);
01045   mpq_ILLfct_compute_xbz (p_mpq->lp);
01046   mpq_ILLfct_check_pfeasible (p_mpq->lp, &fi, mpq_zeroLpNum);
01047   mpq_ILLfct_check_dfeasible (p_mpq->lp, &fi, mpq_zeroLpNum);
01048   mpq_ILLfct_set_status_values (p_mpq->lp, fi.pstatus, fi.dstatus, PHASEII,
01049       PHASEII);
01050   if (p_mpq->lp->basisstat.optimal)
01051   {
01052     *status = QS_LP_OPTIMAL;
01053     EGcallD(mpq_grab_cache (p_mpq, QS_LP_OPTIMAL));
01054   }
01055   else if (p_mpq->lp->basisstat.primal_infeasible
01056       || p_mpq->lp->basisstat.dual_unbounded)
01057   {
01058     if (*status == QS_LP_INFEASIBLE)
01059       *simplexalgo = PRIMAL_SIMPLEX;
01060     *status = QS_LP_INFEASIBLE;
01061     p_mpq->lp->final_phase = PRIMAL_PHASEI;
01062     p_mpq->lp->pIpiz = mpq_EGlpNumAllocArray (p_mpq->lp->nrows);
01063     mpq_ILLfct_compute_phaseI_piz (p_mpq->lp);
01064   }
01065   else if (p_mpq->lp->basisstat.primal_unbounded)
01066     *status = QS_LP_UNBOUNDED;
01067   else
01068     *status = QS_LP_UNSOLVED;
01069   EGtimerStop (&local_timer);
01070   if(!msg_lvl)
01071   {
01072     MESSAGE(0, "Performing Rational Basic Solve on %s, %s, check"
01073         " done in %lg seconds, PS %s %lg, DS %s %lg", p_mpq->name, 
01074         (*status == QS_LP_OPTIMAL) ? "RAT_optimal" : 
01075         ((*status == QS_LP_INFEASIBLE) ?  "RAT_infeasible" : 
01076          ((*status == QS_LP_UNBOUNDED) ?  "RAT_unbounded" : "RAT_unsolved")),
01077         local_timer.time, p_mpq->lp->basisstat.primal_feasible ? 
01078         "F":(p_mpq->lp->basisstat.primal_infeasible ? "I" : "U"), 
01079         p_mpq->lp->basisstat.primal_feasible ?
01080         mpq_get_d(p_mpq->lp->objval) : 
01081         (p_mpq->lp->basisstat.primal_infeasible ?
01082          mpq_get_d(p_mpq->lp->pinfeas) : mpq_get_d(p_mpq->lp->objbound)), 
01083         p_mpq->lp->basisstat.dual_feasible ? 
01084         "F":(p_mpq->lp->basisstat.dual_infeasible ? "I" : "U"), 
01085         p_mpq->lp->basisstat.dual_feasible ? mpq_get_d(p_mpq->lp->dobjval) 
01086         :(p_mpq->lp->basisstat.dual_infeasible ? 
01087           mpq_get_d(p_mpq->lp->dinfeas) : mpq_get_d(p_mpq->lp->objbound)) );
01088   }
01089 CLEANUP:
01090   mpq_EGlpNumClearVar (fi.totinfeas);
01091   return rval;
01092 }
01093 
01094 /* ========================================================================= */
01095 /** @brief test whether given basis is primal and dual feasible in rational arithmetic. 
01096  * @param p_mpq   the problem data.
01097  * @param basis   basis to be tested.
01098  * @param result  where to store whether given basis is primal and dual feasible.
01099  * @param msg_lvl message level.
01100  */
01101 int QSexact_basis_optimalstatus(
01102    mpq_QSdata * p_mpq,
01103    QSbasis* basis,
01104    char* result,
01105    const int msg_lvl
01106    )
01107 {
01108    int rval = 0,
01109       singular;
01110    mpq_feas_info fi;
01111    EGtimer_t local_timer;
01112 
01113    /* test primal and dual feasibility of basic solution */
01114    mpq_EGlpNumInitVar (fi.totinfeas);
01115 
01116    EGtimerReset (&local_timer);
01117    EGtimerStart (&local_timer);
01118 
01119    EGcallD(mpq_QSload_basis (p_mpq, basis));
01120 
01121    if (p_mpq->cache) 
01122    {
01123       mpq_ILLlp_cache_free (p_mpq->cache);
01124       mpq_clear (p_mpq->cache->val);
01125       ILL_IFFREE (p_mpq->cache, mpq_ILLlp_cache); 
01126    }
01127    p_mpq->qstatus = QS_LP_MODIFIED;
01128 
01129    if(p_mpq->qslp->sinfo) 
01130    {
01131       mpq_ILLlp_sinfo_free(p_mpq->qslp->sinfo);
01132       ILL_IFFREE(p_mpq->qslp->sinfo, mpq_ILLlp_sinfo); 
01133    }
01134    if(p_mpq->qslp->rA)
01135    {
01136       mpq_ILLlp_rows_clear (p_mpq->qslp->rA);
01137       ILL_IFFREE (p_mpq->qslp->rA, mpq_ILLlp_rows);
01138    }
01139 
01140    mpq_free_internal_lpinfo (p_mpq->lp);
01141    mpq_init_internal_lpinfo (p_mpq->lp);
01142    EGcallD(mpq_build_internal_lpinfo (p_mpq->lp));
01143 
01144    mpq_ILLfct_set_variable_type (p_mpq->lp);
01145 
01146    EGcallD(mpq_ILLbasis_load (p_mpq->lp, p_mpq->basis));
01147    EGcallD(mpq_ILLbasis_factor (p_mpq->lp, &singular));
01148 
01149    memset (&(p_mpq->lp->basisstat), 0, sizeof (mpq_lp_status_info));
01150    mpq_ILLfct_compute_piz (p_mpq->lp); 
01151    mpq_ILLfct_compute_dz (p_mpq->lp);
01152    mpq_ILLfct_compute_xbz (p_mpq->lp);
01153    mpq_ILLfct_check_pfeasible (p_mpq->lp, &fi, mpq_zeroLpNum);
01154    mpq_ILLfct_check_dfeasible (p_mpq->lp, &fi, mpq_zeroLpNum);
01155    mpq_ILLfct_set_status_values (p_mpq->lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
01156    
01157    if( p_mpq->lp->basisstat.optimal )
01158    {
01159       *result = 1;
01160    }
01161    else
01162    {
01163       *result = 0;
01164    }
01165 
01166    EGtimerStop (&local_timer);
01167 
01168    if( !msg_lvl )
01169    {
01170       MESSAGE(0, "Performing rational solution check for accuratelp on %s, sucess=%s", 
01171          p_mpq->name, 
01172          *result ? "YES" : "NO");
01173    }
01174    
01175  CLEANUP:
01176    mpq_EGlpNumClearVar (fi.totinfeas);
01177    return rval;
01178 }
01179 
01180 /* ========================================================================= */
01181 /** @brief test whether given basis is dual feasible in rational arithmetic. 
01182  * @param p_mpq   the problem data.
01183  * @param basis   basis to be tested.
01184  * @param result  where to store whether given basis is dual feasible.
01185  * @param dobjval where to store dual solution value in case of dual feasibility (if not NULL).
01186  * @param msg_lvl message level.
01187  */
01188 int QSexact_basis_dualstatus(
01189    mpq_QSdata * p_mpq,
01190    QSbasis* basis,
01191    char* result,
01192    mpq_t* dobjval,
01193    const int msg_lvl
01194    )
01195 {
01196   int rval = 0,
01197   singular;
01198   mpq_feas_info fi;
01199   EGtimer_t local_timer;
01200 
01201   mpq_EGlpNumInitVar (fi.totinfeas);
01202   EGtimerReset (&local_timer);
01203   EGtimerStart (&local_timer);
01204   EGcallD(mpq_QSload_basis (p_mpq, basis));
01205   if (p_mpq->cache) 
01206   {
01207     mpq_ILLlp_cache_free (p_mpq->cache);
01208     mpq_clear (p_mpq->cache->val);
01209     ILL_IFFREE (p_mpq->cache, mpq_ILLlp_cache); 
01210   }
01211   p_mpq->qstatus = QS_LP_MODIFIED;
01212 
01213   if(p_mpq->qslp->sinfo) 
01214   {
01215     mpq_ILLlp_sinfo_free(p_mpq->qslp->sinfo);
01216     ILL_IFFREE(p_mpq->qslp->sinfo, mpq_ILLlp_sinfo); 
01217   }
01218 
01219   if(p_mpq->qslp->rA)
01220   {
01221     mpq_ILLlp_rows_clear (p_mpq->qslp->rA);
01222     ILL_IFFREE (p_mpq->qslp->rA, mpq_ILLlp_rows);
01223   }
01224 
01225   mpq_free_internal_lpinfo (p_mpq->lp);
01226   mpq_init_internal_lpinfo (p_mpq->lp);
01227   EGcallD(mpq_build_internal_lpinfo (p_mpq->lp));
01228   mpq_ILLfct_set_variable_type (p_mpq->lp);
01229   EGcallD(mpq_ILLbasis_load (p_mpq->lp, p_mpq->basis));
01230   EGcallD(mpq_ILLbasis_factor (p_mpq->lp, &singular));
01231 
01232   memset (&(p_mpq->lp->basisstat), 0, sizeof (mpq_lp_status_info));
01233   mpq_ILLfct_compute_piz (p_mpq->lp); 
01234   mpq_ILLfct_compute_dz (p_mpq->lp);
01235   mpq_ILLfct_compute_dobj(p_mpq->lp); 
01236   mpq_ILLfct_check_dfeasible (p_mpq->lp, &fi, mpq_zeroLpNum);
01237   mpq_ILLfct_set_status_values (p_mpq->lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
01238 
01239   if( p_mpq->lp->basisstat.dual_feasible )
01240   {
01241     *result = 1;
01242     if( dobjval )
01243     {
01244       mpq_EGlpNumCopy(*dobjval, p_mpq->lp->dobjval);
01245     }
01246   }
01247   else if( p_mpq->lp->basisstat.dual_infeasible )
01248   {
01249     *result = 0;
01250   }
01251   else 
01252   {
01253     TESTG((rval=!(p_mpq->lp->basisstat.dual_unbounded)), CLEANUP, "Internal BUG, problem should be dual unbounded but is not");
01254     *result = 1;
01255     if( dobjval )
01256     {
01257       mpq_EGlpNumCopy(*dobjval, p_mpq->lp->objbound);
01258     }
01259   }
01260 
01261   EGtimerStop (&local_timer);
01262 
01263   if(!msg_lvl)
01264   {
01265     MESSAGE(0, "Performing Rational Basic Test on %s, check done in %lg seconds, DS %s %lg", 
01266         p_mpq->name, local_timer.time, 
01267         p_mpq->lp->basisstat.dual_feasible ? "F": (p_mpq->lp->basisstat.dual_infeasible ? "I" : "U"), 
01268         p_mpq->lp->basisstat.dual_feasible ? mpq_get_d(p_mpq->lp->dobjval) : (p_mpq->lp->basisstat.dual_infeasible ? mpq_get_d(p_mpq->lp->dinfeas) : mpq_get_d(p_mpq->lp->objbound)) );
01269   }
01270 
01271 CLEANUP:
01272   mpq_EGlpNumClearVar (fi.totinfeas);
01273   return rval;
01274 }
01275 
01276 /* ========================================================================= */
01277 /** @brief test whether given basis is dual feasible in rational arithmetic. 
01278  * if wanted it will first directly test the corresponding approximate dual and primal solution 
01279  * (corrected via dual variables for bounds and primal variables for slacks if possible) for optimality
01280  * before performing the dual feasibility test on the more expensive exact basic solution. 
01281  * @param p_mpq   the problem data.
01282  * @param basis   basis to be tested.
01283  * @param useprestep whether to directly test approximate primal and dual solution first.
01284  * @param dbl_p_sol  approximate primal solution to use in prestep 
01285  *                   (NULL in order to compute it by dual simplex in double precision with given starting basis).
01286  * @param dbl_d_sol  approximate dual solution to use in prestep 
01287  *                   (NULL in order to compute it by dual simplex in double precision with given starting basis).
01288  * @param result  where to store whether given basis is dual feasible.
01289  * @param dobjval where to store dual solution value in case of dual feasibility (if not NULL).
01290  * @param msg_lvl message level.
01291  */
01292 int QSexact_verify (
01293    mpq_QSdata * p_mpq,
01294    QSbasis* basis,
01295    int useprestep,
01296    double* dbl_p_sol,
01297    double* dbl_d_sol,
01298    char* result,
01299    mpq_t* dobjval,
01300    const int msg_lvl
01301 )
01302 {
01303    int rval = 0;
01304 
01305    //assert(basis);
01306    //assert(basis->nstruct);
01307 
01308    *result = 0;
01309             
01310    if( useprestep )
01311    {
01312       mpq_t *x_mpq = 0;
01313       mpq_t *y_mpq = 0;
01314       int status = 0;
01315 
01316       if( dbl_p_sol == NULL || dbl_d_sol == NULL  )
01317       {
01318          dbl_QSdata *p_dbl = 0;
01319          double *x_dbl = 0;
01320          double *y_dbl = 0;
01321 
01322          /* create double problem, warmstart with given basis and solve it using double precision 
01323           * this is only done to get approximate primal and dual solution corresponding to the given basis 
01324           */
01325          p_dbl = QScopy_prob_mpq_dbl(p_mpq, "dbl_problem");
01326    
01327          dbl_QSload_basis(p_dbl, basis);
01328          rval = dbl_ILLeditor_solve(p_dbl, DUAL_SIMPLEX);
01329          CHECKRVALG(rval, CLEANUP);
01330       
01331          rval = dbl_QSget_status(p_dbl, &status);
01332          CHECKRVALG(rval, CLEANUP);
01333       
01334          if( status == QS_LP_OPTIMAL )
01335          {
01336             /* get continued fraction approximation of approximate solution */
01337             x_dbl = dbl_EGlpNumAllocArray(p_dbl->qslp->ncols);
01338             y_dbl = dbl_EGlpNumAllocArray(p_dbl->qslp->nrows);
01339 
01340             rval = dbl_QSget_x_array(p_dbl, x_dbl);
01341             CHECKRVALG(rval, CLEANUP);
01342             rval = dbl_QSget_pi_array(p_dbl, y_dbl);
01343             CHECKRVALG(rval, CLEANUP);
01344             x_mpq = QScopy_array_dbl_mpq(x_dbl);
01345             y_mpq = QScopy_array_dbl_mpq(y_dbl);
01346             
01347             /* test optimality of constructed solution */
01348             basis = dbl_QSget_basis(p_dbl);
01349             rval = QSexact_optimal_test(p_mpq, x_mpq, y_mpq, basis);
01350             if( rval )
01351             {
01352                *result = 1;
01353                if( dobjval )
01354                {
01355                   rval = mpq_QSget_objval(p_mpq, dobjval);
01356                   if( rval )
01357                      *result = 0;
01358                }         
01359             }
01360             if( !msg_lvl )
01361             {
01362                MESSAGE(0, "Performing approximated solution check on %s, sucess=%s dobjval=%lg", 
01363                   p_mpq->name, 
01364                   *result ? "YES" : "NO",
01365                   *result ? mpq_get_d(*dobjval) : mpq_get_d(*dobjval));
01366             }
01367          }
01368       CLEANUP:
01369          dbl_EGlpNumFreeArray(x_dbl);
01370          dbl_EGlpNumFreeArray(y_dbl);
01371          mpq_EGlpNumFreeArray(x_mpq);
01372          mpq_EGlpNumFreeArray(y_mpq);
01373          dbl_QSfree_prob(p_dbl);
01374          rval = 0;
01375       }
01376       else
01377       {
01378          dbl_QSdata *p_dbl = 0;
01379          int i;
01380 
01381          /* for some reason, this help to avoid fails in QSexact_basis_dualstatus() after 
01382           * the test here fails, i.e., if we would not perform the test here, than QSexact_basis_dualstatus() would normally not fail
01383           * something happens with the basis... if we do not set up the dbl-prob (?) ????????????????????????
01384           */
01385          p_dbl = QScopy_prob_mpq_dbl(p_mpq, "dbl_problem");
01386          dbl_QSload_basis(p_dbl, basis);
01387 
01388          x_mpq = mpq_EGlpNumAllocArray(p_mpq->qslp->ncols);
01389          y_mpq = mpq_EGlpNumAllocArray(p_mpq->qslp->nrows);
01390 
01391          /* get continued fraction approximation of approximate solution */
01392          for( i = 0; i < p_mpq->qslp->ncols; ++i )
01393             mpq_EGlpNumSet(x_mpq[i], dbl_p_sol[i]);
01394 
01395          for( i = 0; i < p_mpq->qslp->nrows; ++i )
01396             mpq_EGlpNumSet(y_mpq[i], dbl_d_sol[i]);
01397             
01398          /* test optimality of constructed solution */
01399          basis = dbl_QSget_basis(p_dbl);
01400          rval = QSexact_optimal_test(p_mpq, x_mpq, y_mpq, basis);
01401          if( rval )
01402          {
01403             *result = 1;
01404             if( dobjval )
01405             {
01406                rval = mpq_QSget_objval(p_mpq, dobjval);
01407                if( rval )
01408                   *result = 0;
01409             }         
01410          }
01411          if( !msg_lvl )
01412          {
01413             MESSAGE(0, "Performing approximated solution check on %s, sucess=%s dobjval=%lg", 
01414                p_mpq->name, 
01415                *result ? "YES" : "NO",
01416                *result ? mpq_get_d(*dobjval) : mpq_get_d(*dobjval));
01417          }
01418          mpq_EGlpNumFreeArray(x_mpq);
01419          mpq_EGlpNumFreeArray(y_mpq);
01420          dbl_QSfree_prob(p_dbl);
01421          rval = 0;
01422       }
01423    }
01424 
01425    if( !(*result) )
01426    {
01427       rval = QSexact_basis_dualstatus(p_mpq, basis, result, dobjval, msg_lvl);
01428       if( !msg_lvl )
01429       {
01430          MESSAGE(0, "Performing rational solution check on %s, sucess=%s dobjval=%lg", 
01431             p_mpq->name, 
01432             *result ? "YES" : "NO",
01433             *result ? mpq_get_d(*dobjval) : mpq_get_d(*dobjval));
01434       }
01435    }
01436 
01437    return rval;
01438 }
01439 
01440 /* ========================================================================= */
01441 int QSexact_solver (mpq_QSdata * p_mpq,
01442                     mpq_t * const x,
01443                     mpq_t * const y,
01444                     QSbasis * const ebasis,
01445                     int simplexalgo,
01446                     int *status)
01447 {
01448   /* local variables */
01449   int last_status = 0, last_iter = 0;
01450   QSbasis *basis = 0;
01451   unsigned precision = EGLPNUM_PRECISION;
01452   int rval = 0,
01453     it = QS_EXACT_MAX_ITER;
01454   dbl_QSdata *p_dbl = 0;
01455   mpf_QSdata *p_mpf = 0;
01456   double *x_dbl = 0,
01457    *y_dbl = 0;
01458   mpq_t *x_mpq = 0,
01459    *y_mpq = 0;
01460   mpf_t *x_mpf = 0,
01461    *y_mpf = 0;
01462   int const msg_lvl = __QS_SB_VERB <= DEBUG ? 0: (1 - p_mpq->simplex_display) * 10000;
01463   *status = 0;
01464   /* save the problem if we are really debugging */
01465   if(DEBUG >= __QS_SB_VERB)
01466   {
01467     EGcallD(mpq_QSwrite_prob(p_mpq, "qsxprob.lp","LP"));
01468   }
01469   /* try first with doubles */
01470   if (p_mpq->simplex_display || DEBUG >= __QS_SB_VERB)
01471   {
01472     fprintf (stdout, "%s\n\tTrying double precision\n%s\n", __sp, __sp);
01473     fflush(stdout);
01474   }
01475   p_dbl = QScopy_prob_mpq_dbl (p_mpq, "dbl_problem");
01476   if(__QS_SB_VERB <= DEBUG) p_dbl->simplex_display = 1;
01477   if (ebasis && ebasis->nstruct)
01478     dbl_QSload_basis (p_dbl, ebasis);
01479   if (dbl_ILLeditor_solve (p_dbl, simplexalgo))
01480   {
01481     MESSAGE(p_mpq->simplex_display ? 0: __QS_SB_VERB, 
01482             "double approximation failed, code %d, "
01483             "continuing in extended precision", rval);
01484     goto MPF_PRECISION;
01485   }
01486   EGcallD(dbl_QSget_status (p_dbl, status));
01487   if ((*status == QS_LP_INFEASIBLE) &&
01488       (p_dbl->lp->final_phase != PRIMAL_PHASEI) &&
01489       (p_dbl->lp->final_phase != DUAL_PHASEII))
01490     dbl_QSopt_primal (p_dbl, status);
01491   EGcallD(dbl_QSget_status (p_dbl, status));
01492   last_status = *status;
01493   EGcallD(dbl_QSget_itcnt(p_dbl, 0, 0, 0, 0, &last_iter));
01494   /* deal with the problem depending on what status we got from our optimizer */
01495   switch (*status)
01496   {
01497   case QS_LP_OPTIMAL:
01498     x_dbl = dbl_EGlpNumAllocArray (p_dbl->qslp->ncols);
01499     y_dbl = dbl_EGlpNumAllocArray (p_dbl->qslp->nrows);
01500     EGcallD(dbl_QSget_x_array (p_dbl, x_dbl));
01501     EGcallD(dbl_QSget_pi_array (p_dbl, y_dbl));
01502     x_mpq = QScopy_array_dbl_mpq (x_dbl);
01503     y_mpq = QScopy_array_dbl_mpq (y_dbl);
01504     dbl_EGlpNumFreeArray (x_dbl);
01505     dbl_EGlpNumFreeArray (y_dbl);
01506     basis = dbl_QSget_basis (p_dbl);
01507     if (QSexact_optimal_test (p_mpq, x_mpq, y_mpq, basis))
01508     {
01509       optimal_output (p_mpq, x, y, x_mpq, y_mpq);
01510       goto CLEANUP;
01511     }
01512     else
01513     {
01514       EGcallD(QSexact_basis_status (p_mpq, status, basis, msg_lvl, &simplexalgo));
01515       if (*status == QS_LP_OPTIMAL)
01516       {
01517         if(!msg_lvl)
01518         {
01519           MESSAGE(0,"Retesting solution");
01520         }
01521         EGcallD(mpq_QSget_x_array (p_mpq, x_mpq));
01522         EGcallD(mpq_QSget_pi_array (p_mpq, y_mpq));
01523         if (QSexact_optimal_test (p_mpq, x_mpq, y_mpq, basis))
01524         {
01525           optimal_output (p_mpq, x, y, x_mpq, y_mpq);
01526           goto CLEANUP;
01527         }
01528         else
01529         {
01530           last_status = *status = QS_LP_UNSOLVED;
01531         }
01532       }
01533       else
01534       {
01535         if(!msg_lvl)
01536         {
01537           MESSAGE(0,"Status is not optimal, but %d", *status);
01538         }
01539       }
01540     }
01541     mpq_EGlpNumFreeArray (x_mpq);
01542     mpq_EGlpNumFreeArray (y_mpq);
01543     break;
01544   case QS_LP_INFEASIBLE:
01545     y_dbl = dbl_EGlpNumAllocArray (p_dbl->qslp->nrows);
01546     if (dbl_QSget_infeas_array (p_dbl, y_dbl))
01547     {
01548       MESSAGE(p_mpq->simplex_display ? 0 : __QS_SB_VERB, "double approximation"
01549               " failed, code %d, continuing in extended precision\n", rval);
01550       goto MPF_PRECISION;
01551     }
01552     y_mpq = QScopy_array_dbl_mpq (y_dbl);
01553     dbl_EGlpNumFreeArray (y_dbl);
01554     if (QSexact_infeasible_test (p_mpq, y_mpq))
01555     {
01556       infeasible_output (p_mpq, y, y_mpq);
01557       goto CLEANUP;
01558     }
01559     else
01560     {
01561       MESSAGE (msg_lvl, "Retesting solution in exact arithmetic");
01562       basis = dbl_QSget_basis (p_dbl);
01563       EGcallD(QSexact_basis_status (p_mpq, status, basis, msg_lvl, &simplexalgo));
01564       #if 0
01565       mpq_QSset_param (p_mpq, QS_PARAM_SIMPLEX_MAX_ITERATIONS, 1);
01566       mpq_QSload_basis (p_mpq, basis);
01567       mpq_QSfree_basis (basis);
01568       EGcallD(mpq_ILLeditor_solve (p_mpq, simplexalgo));
01569       EGcallD(mpq_QSget_status (p_mpq, status));
01570       #endif
01571       if (*status == QS_LP_INFEASIBLE)
01572       {
01573         mpq_EGlpNumFreeArray (y_mpq);
01574         y_mpq = mpq_EGlpNumAllocArray (p_mpq->qslp->nrows);
01575         EGcallD(mpq_QSget_infeas_array (p_mpq, y_mpq));
01576         if (QSexact_infeasible_test (p_mpq, y_mpq))
01577         {
01578           infeasible_output (p_mpq, y, y_mpq);
01579           goto CLEANUP;
01580         }
01581         else
01582         {
01583           last_status = *status = QS_LP_UNSOLVED;
01584         }
01585       }
01586     }
01587     mpq_EGlpNumFreeArray (y_mpq);
01588     break;
01589   case QS_LP_UNBOUNDED:
01590     MESSAGE(p_mpq->simplex_display ? 0 : __QS_SB_VERB, "%s\n\tUnbounded "
01591             "Problem found, not implemented to deal with this\n%s\n",__sp,__sp);
01592     break;
01593   case QS_LP_OBJ_LIMIT:
01594     rval=1;
01595     IFMESSAGE(p_mpq->simplex_display,"Objective limit reached (in floating point) ending now");
01596     goto CLEANUP;
01597     break;
01598   default:
01599     IFMESSAGE(p_mpq->simplex_display,"Re-trying inextended precision");
01600     break;
01601   }
01602   /* if we reach this point, then we have to keep going, we use the previous
01603    * basis ONLY if the previous precision think that it has the optimal
01604    * solution, otherwise we start from scratch. */
01605   precision = 128;
01606   MPF_PRECISION:
01607   dbl_QSfree_prob (p_dbl);
01608   p_dbl = 0;
01609   /* try with multiple precision floating points */
01610   for (; it--; precision = (unsigned) (precision * 1.5))
01611   {
01612     QSexact_set_precision (precision);
01613     if (p_mpq->simplex_display || DEBUG >= __QS_SB_VERB)
01614     {
01615       fprintf (stdout, "%s\n\tTrying mpf with %u bits\n%s\n", __sp, precision,
01616                __sp);
01617       fflush(stdout);
01618     }
01619     p_mpf = QScopy_prob_mpq_mpf (p_mpq, "mpf_problem");
01620     if(DEBUG >= __QS_SB_VERB)
01621     {
01622       EGcallD(mpf_QSwrite_prob(p_mpf, "qsxprob.mpf.lp","LP"));
01623     }
01624     if(__QS_SB_VERB <= DEBUG) p_mpf->simplex_display = 1;
01625     simplexalgo = PRIMAL_SIMPLEX;
01626     if(!last_iter) last_status = QS_LP_UNSOLVED;
01627     if(last_status == QS_LP_OPTIMAL || last_status == QS_LP_INFEASIBLE)
01628     {
01629       if (p_mpq->simplex_display || DEBUG >= __QS_SB_VERB)
01630       {
01631         fprintf(stdout,"Re-using previous basis\n");
01632         fflush(stdout);
01633       }
01634       if (basis)
01635       {
01636         EGcallD(mpf_QSload_basis (p_mpf, basis));
01637         mpf_QSfree_basis (basis);
01638         simplexalgo = DUAL_SIMPLEX;
01639         basis = 0;
01640       }
01641       else if (ebasis && ebasis->nstruct)
01642       {
01643         mpf_QSload_basis (p_mpf, ebasis);
01644         simplexalgo = DUAL_SIMPLEX;
01645       }
01646     }
01647     else
01648     {
01649       if(p_mpf->basis)
01650       {
01651         mpf_ILLlp_basis_free(p_mpf->basis);
01652         p_mpf->lp->basisid = -1;
01653         p_mpf->factorok = 0;
01654       }
01655       if (p_mpq->simplex_display || DEBUG >= __QS_SB_VERB)
01656       {
01657         fprintf(stdout,"Not-using previous basis\n");
01658         fflush(stdout);
01659       }
01660     }
01661     if (mpf_ILLeditor_solve (p_mpf, simplexalgo))
01662     {
01663       if (p_mpq->simplex_display || DEBUG >= __QS_SB_VERB)
01664       {
01665         fprintf (stdout,
01666                  "mpf_%u precision falied, error code %d, continuing with "
01667                  "next precision", precision, rval);
01668         fflush(stdout);
01669        }
01670       goto NEXT_PRECISION;
01671     }
01672     EGcallD(mpf_QSget_status (p_mpf, status));
01673     if ((*status == QS_LP_INFEASIBLE) &&
01674         (p_mpf->lp->final_phase != PRIMAL_PHASEI) &&
01675         (p_mpf->lp->final_phase != DUAL_PHASEII))
01676       mpf_QSopt_primal (p_mpf, status);
01677     EGcallD(mpf_QSget_status (p_mpf, status));
01678     last_status = *status;
01679     EGcallD(mpf_QSget_itcnt(p_mpf, 0, 0, 0, 0, &last_iter));
01680     /* deal with the problem depending on status we got from our optimizer */
01681     switch (*status)
01682     {
01683     case QS_LP_OPTIMAL:
01684       basis = mpf_QSget_basis (p_mpf);
01685       x_mpf = mpf_EGlpNumAllocArray (p_mpf->qslp->ncols);
01686       y_mpf = mpf_EGlpNumAllocArray (p_mpf->qslp->nrows);
01687       EGcallD(mpf_QSget_x_array (p_mpf, x_mpf));
01688       EGcallD(mpf_QSget_pi_array (p_mpf, y_mpf));
01689       x_mpq = QScopy_array_mpf_mpq (x_mpf);
01690       y_mpq = QScopy_array_mpf_mpq (y_mpf);
01691       mpf_EGlpNumFreeArray (x_mpf);
01692       mpf_EGlpNumFreeArray (y_mpf);
01693       if (QSexact_optimal_test (p_mpq, x_mpq, y_mpq, basis))
01694       {
01695         optimal_output (p_mpq, x, y, x_mpq, y_mpq);
01696         goto CLEANUP;
01697       }
01698       else
01699       {
01700         EGcallD(QSexact_basis_status (p_mpq, status, basis, msg_lvl, &simplexalgo));
01701         if (*status == QS_LP_OPTIMAL)
01702         {
01703           MESSAGE (msg_lvl, "Retesting solution");
01704           EGcallD(mpq_QSget_x_array (p_mpq, x_mpq));
01705           EGcallD(mpq_QSget_pi_array (p_mpq, y_mpq));
01706           if (QSexact_optimal_test (p_mpq, x_mpq, y_mpq, basis))
01707           {
01708             optimal_output (p_mpq, x, y, x_mpq, y_mpq);
01709             goto CLEANUP;
01710           }
01711           else
01712           {
01713             last_status = *status = QS_LP_UNSOLVED;
01714           }
01715         }
01716         else
01717           MESSAGE (msg_lvl, "Status is not optimal, but %d", *status);
01718       }
01719       mpq_EGlpNumFreeArray (x_mpq);
01720       mpq_EGlpNumFreeArray (y_mpq);
01721       break;
01722     case QS_LP_INFEASIBLE:
01723       y_mpf = mpf_EGlpNumAllocArray (p_mpf->qslp->nrows);
01724       EGcallD(mpf_QSget_infeas_array (p_mpf, y_mpf));
01725       y_mpq = QScopy_array_mpf_mpq (y_mpf);
01726       mpf_EGlpNumFreeArray (y_mpf);
01727       if (QSexact_infeasible_test (p_mpq, y_mpq))
01728       {
01729         infeasible_output (p_mpq, y, y_mpq);
01730         goto CLEANUP;
01731       }
01732       else
01733       {
01734         MESSAGE (msg_lvl, "Retesting solution in exact arithmetic");
01735         basis = mpf_QSget_basis (p_mpf);
01736         EGcallD(QSexact_basis_status (p_mpq, status, basis, msg_lvl, &simplexalgo));
01737 #if 0
01738         mpq_QSset_param (p_mpq, QS_PARAM_SIMPLEX_MAX_ITERATIONS, 1);
01739         mpq_QSload_basis (p_mpq, basis);
01740         mpq_QSfree_basis (basis);
01741         EGcallD(mpq_ILLeditor_solve (p_mpq, simplexalgo));
01742         EGcallD(mpq_QSget_status (p_mpq, status));
01743 #endif
01744         if (*status == QS_LP_INFEASIBLE)
01745         {
01746           mpq_EGlpNumFreeArray (y_mpq);
01747           y_mpq = mpq_EGlpNumAllocArray (p_mpq->qslp->nrows);
01748           EGcallD(mpq_QSget_infeas_array (p_mpq, y_mpq));
01749           if (QSexact_infeasible_test (p_mpq, y_mpq))
01750           {
01751             infeasible_output (p_mpq, y, y_mpq);
01752             goto CLEANUP;
01753           }
01754           else
01755           {
01756             last_status = *status = QS_LP_UNSOLVED;
01757           }
01758         }
01759       }
01760       mpq_EGlpNumFreeArray (y_mpq);
01761       break;
01762       break;
01763     case QS_LP_OBJ_LIMIT:
01764       rval=1;
01765       IFMESSAGE(p_mpq->simplex_display,"Objective limit reached (in floating point) ending now");
01766       goto CLEANUP;
01767       break;
01768     case QS_LP_UNBOUNDED:
01769     default:
01770       MESSAGE(__QS_SB_VERB,"Re-trying inextended precision");
01771       break;
01772     }
01773   NEXT_PRECISION:
01774     mpf_QSfree_prob (p_mpf);
01775     p_mpf = 0;
01776   }
01777   /* ending */
01778 CLEANUP:
01779   dbl_EGlpNumFreeArray (x_dbl);
01780   dbl_EGlpNumFreeArray (y_dbl);
01781   mpq_EGlpNumFreeArray (x_mpq);
01782   mpq_EGlpNumFreeArray (y_mpq);
01783   mpf_EGlpNumFreeArray (x_mpf);
01784   mpf_EGlpNumFreeArray (y_mpf);
01785   if (ebasis && basis)
01786   {
01787     ILL_IFFREE (ebasis->cstat, char);
01788     ILL_IFFREE (ebasis->rstat, char);
01789     ebasis->nstruct = basis->nstruct;
01790     ebasis->nrows = basis->nrows;
01791     ebasis->cstat = basis->cstat;
01792     ebasis->rstat = basis->rstat;
01793     basis->cstat = basis->rstat = 0;
01794   }
01795   mpq_QSfree_basis (basis);
01796   dbl_QSfree_prob (p_dbl);
01797   mpf_QSfree_prob (p_mpf);
01798   return rval;
01799 }
01800 
01801 /* ========================================================================= */
01802 int __QSexact_setup = 0;
01803 /* ========================================================================= */
01804 void QSexactStart(void)
01805 {
01806   /* if we have been initialized before, do nothing */
01807   if(__QSexact_setup) return;
01808   /* we should call EGlpNumStart() */
01809   EGlpNumStart();
01810   /* now we call all setups */
01811   EXutilDoInit();
01812   dbl_ILLstart();
01813   mpf_ILLstart();
01814   mpq_ILLstart();
01815   fp20_ILLstart();
01816   #ifdef HAVE_SOFTFLOAT
01817   float128_ILLstart();
01818   #endif
01819   #if ENABLE_LONG_DOUBLE
01820   ldbl_ILLstart();
01821   #endif
01822   /* ending */
01823   __QSexact_setup = 1;
01824 }
01825 /* ========================================================================= */
01826 void QSexactClear(void)
01827 {
01828   if(!__QSexact_setup) return;
01829   /* now we call all ends */
01830   #ifdef HAVE_SOFTFLOAT
01831   float128_ILLend();
01832   #endif
01833   #if ENABLE_LONG_DOUBLE
01834   ldbl_ILLend();
01835   #endif
01836   fp20_ILLend();
01837   dbl_ILLend();
01838   mpf_ILLend();
01839   mpq_ILLend();
01840   EXutilDoClear();
01841   /* ending */
01842   EGlpNumClear();
01843   __QSexact_setup = 0;
01844 }
01845 /* ========================================================================= */
01846 /** @} */
01847 /* end of exact.c */
01848 

Generated on Thu Mar 29 09:32:29 2012 for QSopt_ex by  doxygen 1.4.7