mpf_simplex.c

Go to the documentation of this file.
00001 /****************************************************************************/
00002 /*                                                                          */
00003 /*  This file is part of QSopt_ex.                                          */
00004 /*                                                                          */
00005 /*  (c) Copyright 2006 by David Applegate, William Cook, Sanjeeb Dash,      */
00006 /*  and Daniel Espinoza                                                     */
00007 /*                                                                          */
00008 /*  Sanjeeb Dash ownership of copyright in QSopt_ex is derived from his     */
00009 /*  copyright in QSopt.                                                     */
00010 /*                                                                          */
00011 /*  This code may be used under the terms of the GNU General Public License */
00012 /*  (Version 2.1 or later) as published by the Free Software Foundation.    */
00013 /*                                                                          */
00014 /*  Alternatively, use is granted for research purposes only.               */
00015 /*                                                                          */
00016 /*  It is your choice of which of these two licenses you are operating      */
00017 /*  under.                                                                  */
00018 /*                                                                          */
00019 /*  We make no guarantees about the correctness or usefulness of this code. */
00020 /*                                                                          */
00021 /****************************************************************************/
00022 
00023 /* RCS_INFO = "$RCSfile: mpf_simplex.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 static int TRACE = 0;
00025 
00026 #define mpf_QSOPT_CURRENT_PRECICION
00027 #include "basicdefs.h"
00028 #include "config.h"
00029 #include "mpf_iqsutil.h"
00030 #include "mpf_lpdata.h"
00031 #include "mpf_lpdefs.h"
00032 
00033 #include "stddefs.h"
00034 #include "mpf_fct.h"
00035 #include "mpf_ratio.h"
00036 #include "mpf_price.h"
00037 #include "mpf_basis.h"
00038 #include "mpf_simplex.h"
00039 #include "mpf_dstruct.h"
00040 #include "mpf_qstruct.h"
00041 #include "mpf_qsopt.h"
00042 #include "mpf_lib.h"                /* for mpf_ILLlib_writebasis */
00043 #include "mpf_lp.h"                 /* for mpf_ILLwrite_lp */
00044 
00045 #ifdef USEDMALLOC
00046 #include "dmalloc.h"
00047 #endif
00048 
00049 static void mpf_init_lp_status_info (
00050   mpf_lp_status_info * ls),
00051   mpf_init_simplex_tols (
00052   mpf_lpinfo * lp),
00053   mpf_monitor_iter (
00054   mpf_lpinfo * lp,
00055   mpf_price_info * p,
00056   mpf_iter_info * it,
00057   int cphase),
00058   mpf_get_current_stat (
00059   mpf_lp_status_info * p,
00060   int algorithm,
00061   int *bstat);
00062 
00063 static int mpf_terminate_simplex (
00064   mpf_lpinfo * lp,
00065   int phase,
00066   mpf_iter_info * it),
00067   mpf_primal_phaseI_step (
00068   mpf_lpinfo * lp,
00069   mpf_price_info * pinf,
00070   mpf_svector * updz,
00071   mpf_svector * wz,
00072   mpf_iter_info * it),
00073   mpf_primal_phaseII_step (
00074   mpf_lpinfo * lp,
00075   mpf_price_info * pinf,
00076   mpf_svector * updz,
00077   mpf_svector * wz,
00078   mpf_iter_info * it),
00079   mpf_dual_phaseI_step (
00080   mpf_lpinfo * lp,
00081   mpf_price_info * pinf,
00082   mpf_svector * updz,
00083   mpf_svector * wz,
00084   mpf_iter_info * it),
00085   mpf_dual_phaseII_step (
00086   mpf_lpinfo * lp,
00087   mpf_price_info * pinf,
00088   mpf_svector * updz,
00089   mpf_svector * wz,
00090   mpf_iter_info * it),
00091   mpf_report_value (
00092   mpf_lpinfo * lp,
00093   mpf_iter_info * it,
00094   const char *value_name,
00095   mpf_t value);
00096 
00097 
00098 void mpf_ILLsimplex_init_lpinfo (
00099   mpf_lpinfo * lp)
00100 {
00101   mpf_ILLbasis_init_basisinfo (lp);
00102   mpf_init_internal_lpinfo (lp);
00103 }
00104 
00105 void mpf_ILLsimplex_free_lpinfo (
00106   mpf_lpinfo * lp)
00107 {
00108   if (lp)
00109   {
00110     mpf_EGlpNumFreeArray (lp->lz);
00111     mpf_EGlpNumFreeArray (lp->uz);
00112     mpf_EGlpNumFreeArray (lp->cz);
00113     mpf_ILLbasis_free_basisinfo (lp);
00114     mpf_free_internal_lpinfo (lp);
00115   }
00116 }
00117 
00118 void mpf_ILLsimplex_load_lpinfo (
00119   mpf_ILLlpdata * qslp,
00120   mpf_lpinfo * lp)
00121 {
00122   lp->basisid = -1;
00123   lp->maxiter = 500000;
00124   lp->maxtime = 300000;
00125   //lp->iterskip = 10;
00126   lp->iterskip = 100;
00127   mpf_EGlpNumCopy (lp->objbound, mpf_INFTY);
00128   lp->O = qslp;
00129 }
00130 
00131 void mpf_ILLsimplex_set_bound (
00132   mpf_lpinfo * lp,
00133   const mpf_t * objbound,
00134   int sense)
00135 {
00136   mpf_EGlpNumCopy (lp->objbound, *objbound);
00137   if (sense == mpf_ILL_MAX)
00138     mpf_EGlpNumSign (lp->objbound);
00139 }
00140 
00141 static void mpf_init_lp_status_info (
00142   mpf_lp_status_info * ls)
00143 {
00144   ls->optimal = 0;
00145   ls->primal_feasible = 0;
00146   ls->primal_infeasible = 0;
00147   ls->primal_unbounded = 0;
00148   ls->dual_feasible = 0;
00149   ls->dual_infeasible = 0;
00150   ls->dual_unbounded = 0;
00151 }
00152 
00153 static void mpf_init_simplex_tols (
00154   mpf_lpinfo * lp)
00155 {
00156   mpf_EGlpNumCopy (lp->tol->pfeas_tol, mpf_PFEAS_TOLER);
00157   mpf_EGlpNumCopy (lp->tol->dfeas_tol, mpf_DFEAS_TOLER);
00158   mpf_EGlpNumCopy (lp->tol->pivot_tol, mpf_PIVOT_TOLER);
00159   mpf_EGlpNumCopy (lp->tol->szero_tol, mpf_SZERO_TOLER);
00160   mpf_EGlpNumCopy (lp->tol->ip_tol, lp->tol->pfeas_tol);
00161   mpf_EGlpNumCopy (lp->tol->id_tol, lp->tol->dfeas_tol);
00162   if (mpf_EGlpNumIsNeqqZero (lp->tol->ip_tol))
00163   {
00164 #if VERBOSE_LEVEL <= DEBUG
00165     MESSAGE (VERBOSE_LEVEL, "ip_tol %lg", mpf_EGlpNumToLf (lp->tol->ip_tol));
00166     MESSAGE (VERBOSE_LEVEL, "eps %lg", mpf_EGlpNumToLf (mpf_epsLpNum));
00167     MESSAGE (VERBOSE_LEVEL, "mpf_PFEAS_TOLER %lg", mpf_EGlpNumToLf (mpf_PFEAS_TOLER));
00168 #endif
00169     mpf_EGlpNumDivUiTo (lp->tol->ip_tol, 2UL);
00170   }
00171   if (mpf_EGlpNumIsNeqqZero (lp->tol->id_tol))
00172   {
00173 #if VERBOSE_LEVEL <= DEBUG
00174     MESSAGE (VERBOSE_LEVEL, "id_tol %lg", mpf_EGlpNumToLf (lp->tol->id_tol));
00175 #endif
00176     mpf_EGlpNumDivUiTo (lp->tol->id_tol, 2UL);
00177   }
00178 }
00179 
00180 void mpf_init_internal_lpinfo (
00181   mpf_lpinfo * lp)
00182 {
00183   int rval = 0;
00184 
00185   lp->nrows = 0;
00186   lp->nnbasic = 0;
00187   lp->localrows = 0;
00188   lp->rowcnt = 0;
00189   lp->rowbeg = 0;
00190   lp->rowind = 0;
00191   lp->rowval = 0;
00192   lp->cz = 0;
00193   lp->lz = 0;
00194   lp->uz = 0;
00195   lp->xbz = 0;
00196   lp->piz = 0;
00197   lp->dz = 0;
00198   lp->pIxbz = 0;
00199   lp->pIpiz = 0;
00200   lp->pIdz = 0;
00201   lp->vtype = 0;
00202   lp->vclass = 0;
00203   lp->iwork = 0;
00204   lp->upd.perm = 0;
00205   lp->upd.ix = 0;
00206   lp->upd.t = 0;
00207   lp->bfeas = 0;
00208   lp->dfeas = 0;
00209   lp->tol = 0;
00210   lp->cnts = 0;
00211   lp->bchanges = 0;
00212   lp->cchanges = 0;
00213   mpf_ILLsvector_init (&(lp->zz));
00214   mpf_ILLsvector_init (&(lp->yjz));
00215   mpf_ILLsvector_init (&(lp->zA));
00216   mpf_ILLsvector_init (&(lp->work));
00217   mpf_ILLsvector_init (&(lp->srhs));
00218   mpf_ILLsvector_init (&(lp->ssoln));
00219   ILL_SAFE_MALLOC (lp->tol, 1, mpf_tol_struct);
00220   mpf_EGlpNumInitVar (lp->tol->pfeas_tol);
00221   mpf_EGlpNumInitVar (lp->tol->dfeas_tol);
00222   mpf_EGlpNumInitVar (lp->tol->pivot_tol);
00223   mpf_EGlpNumInitVar (lp->tol->szero_tol);
00224   mpf_EGlpNumInitVar (lp->tol->ip_tol);
00225   mpf_EGlpNumInitVar (lp->tol->id_tol);
00226   ILL_SAFE_MALLOC (lp->cnts, 1, mpf_count_struct);
00227   mpf_EGlpNumInitVar (lp->cnts->y_ravg);
00228   mpf_EGlpNumInitVar (lp->cnts->z_ravg);
00229   mpf_EGlpNumInitVar (lp->cnts->za_ravg);
00230 CLEANUP:
00231   if (rval)
00232   {
00233     fprintf (stderr, "\nno memory, in %s, exit\n", __func__);
00234     exit (1);
00235   }
00236 }
00237 
00238 void mpf_free_internal_lpinfo (
00239   mpf_lpinfo * lp)
00240 {
00241   mpf_bndinfo *binfo = 0;
00242   mpf_coefinfo *cinfo = 0;
00243 
00244   if (lp->localrows)
00245   {
00246     ILL_IFFREE (lp->rowcnt, int);
00247     ILL_IFFREE (lp->rowbeg, int);
00248     ILL_IFFREE (lp->rowind, int);
00249 
00250     mpf_EGlpNumFreeArray (lp->rowval);
00251     lp->localrows = 0;
00252   }
00253   mpf_EGlpNumFreeArray (lp->lz);
00254   mpf_EGlpNumFreeArray (lp->uz);
00255   mpf_EGlpNumFreeArray (lp->cz);
00256   mpf_EGlpNumFreeArray (lp->xbz);
00257   mpf_EGlpNumFreeArray (lp->piz);
00258   mpf_EGlpNumFreeArray (lp->pIpiz);
00259   mpf_EGlpNumFreeArray (lp->dz);
00260   mpf_EGlpNumFreeArray (lp->pIdz);
00261   mpf_EGlpNumFreeArray (lp->pIxbz);
00262 
00263   ILL_IFFREE (lp->vtype, int);
00264   ILL_IFFREE (lp->vclass, char);
00265 
00266   mpf_ILLsvector_free (&(lp->zz));
00267   mpf_ILLsvector_free (&(lp->yjz));
00268   mpf_ILLsvector_free (&(lp->zA));
00269   mpf_ILLsvector_free (&(lp->work));
00270   mpf_ILLsvector_free (&(lp->srhs));
00271   mpf_ILLsvector_free (&(lp->ssoln));
00272   ILL_IFFREE (lp->iwork, int);
00273   ILL_IFFREE (lp->upd.perm, int);
00274   ILL_IFFREE (lp->upd.ix, int);
00275 
00276   mpf_EGlpNumFreeArray (lp->upd.t);
00277 
00278   ILL_IFFREE (lp->bfeas, int);
00279   ILL_IFFREE (lp->dfeas, int);
00280 
00281   if (lp->tol)
00282   {
00283     mpf_EGlpNumClearVar (lp->tol->pfeas_tol);
00284     mpf_EGlpNumClearVar (lp->tol->dfeas_tol);
00285     mpf_EGlpNumClearVar (lp->tol->pivot_tol);
00286     mpf_EGlpNumClearVar (lp->tol->szero_tol);
00287     mpf_EGlpNumClearVar (lp->tol->ip_tol);
00288     mpf_EGlpNumClearVar (lp->tol->id_tol);
00289     ILL_IFFREE (lp->tol, mpf_tol_struct);
00290   }
00291   if (lp->cnts)
00292   {
00293     mpf_EGlpNumClearVar (lp->cnts->y_ravg);
00294     mpf_EGlpNumClearVar (lp->cnts->z_ravg);
00295     mpf_EGlpNumClearVar (lp->cnts->za_ravg);
00296     ILL_IFFREE (lp->cnts, mpf_count_struct);
00297   }
00298 
00299   while (lp->bchanges)
00300   {
00301     binfo = lp->bchanges;
00302     mpf_EGlpNumClearVar (binfo->pbound);
00303     mpf_EGlpNumClearVar (binfo->cbound);
00304     lp->bchanges = binfo->next;
00305     ILL_IFFREE (binfo, mpf_bndinfo);
00306   }
00307 
00308   while (lp->cchanges)
00309   {
00310     cinfo = lp->cchanges;
00311     mpf_EGlpNumClearVar (cinfo->pcoef);
00312     mpf_EGlpNumClearVar (cinfo->ccoef);
00313     lp->cchanges = cinfo->next;
00314     ILL_IFFREE (cinfo, mpf_coefinfo);
00315   }
00316 }
00317 
00318 int mpf_build_internal_lpinfo (
00319   mpf_lpinfo * lp)
00320 {
00321   int rval = 0;
00322   int i, n;
00323   mpf_ILLlpdata *qslp = lp->O;
00324   mpf_ILLlp_sinfo *S = lp->O->sinfo;
00325   mpf_t *lower, *upper, *obj;
00326   mpf_ILLlp_rows lprows;
00327   mpf_ILLmatrix *A;
00328 
00329   mpf_init_lp_status_info (&(lp->probstat));
00330   mpf_init_lp_status_info (&(lp->basisstat));
00331 
00332   if (S != 0)
00333   {
00334     lp->nrows = S->nrows;
00335     lp->ncols = S->ncols;
00336     lp->bz = S->rhs;
00337     lower = S->lower;
00338     upper = S->upper;
00339     obj = S->obj;
00340     A = &(S->A);
00341   }
00342   else
00343   {
00344     lp->nrows = qslp->nrows;
00345     lp->ncols = qslp->ncols;
00346     lp->bz = qslp->rhs;
00347     lower = qslp->lower;
00348     upper = qslp->upper;
00349     obj = qslp->obj;
00350     A = &(qslp->A);
00351   }
00352 
00353   lp->matbeg = A->matbeg;
00354   lp->matcnt = A->matcnt;
00355   lp->matind = A->matind;
00356   lp->matval = A->matval;
00357 
00358   lp->nnbasic = lp->ncols - lp->nrows;
00359 
00360   lp->lz = mpf_EGlpNumAllocArray (lp->ncols);
00361   lp->cz = mpf_EGlpNumAllocArray (lp->ncols);
00362   lp->uz = mpf_EGlpNumAllocArray (lp->ncols);
00363   if (!lp->lz || !lp->uz || !lp->cz)
00364   {
00365     fprintf (stderr, "mpf_build_internal_lpinfo\n");
00366     rval = 1;
00367     goto CLEANUP;
00368   }
00369   for (i = 0; i < lp->ncols; i++)
00370   {
00371     mpf_EGlpNumCopy (lp->lz[i], lower[i]);
00372     mpf_EGlpNumCopy (lp->uz[i], upper[i]);
00373     mpf_EGlpNumCopy (lp->cz[i], obj[i]);
00374     if (qslp->objsense == mpf_ILL_MAX)
00375     {
00376       mpf_EGlpNumSign (lp->cz[i]);
00377     }
00378   }
00379 
00380   if (!lp->O->rA)
00381   {
00382     rval = mpf_ILLlp_rows_init (&lprows, lp->O, 1);
00383     CHECKRVALG (rval, CLEANUP);
00384     lp->rowbeg = lprows.rowbeg;
00385     lp->rowcnt = lprows.rowcnt;
00386     lp->rowind = lprows.rowind;
00387     lp->rowval = lprows.rowval;
00388     lp->localrows = 1;
00389   }
00390   else
00391   {
00392     /* row format exists, just use pointers */
00393     lp->rowbeg = lp->O->rA->rowbeg;
00394     lp->rowcnt = lp->O->rA->rowcnt;
00395     lp->rowind = lp->O->rA->rowind;
00396     lp->rowval = lp->O->rA->rowval;
00397     lp->localrows = 0;
00398   }
00399 
00400   lp->xbz = mpf_EGlpNumAllocArray (lp->nrows);
00401   lp->piz = mpf_EGlpNumAllocArray (lp->nrows);
00402   lp->dz = mpf_EGlpNumAllocArray (lp->nnbasic);
00403   lp->final_phase = -1;
00404   lp->infub_ix = -1;
00405 
00406   ILL_SAFE_MALLOC (lp->vtype, lp->ncols, int);
00407   ILL_SAFE_MALLOC (lp->vclass, lp->ncols, char);
00408 
00409   rval = mpf_ILLsvector_alloc (&(lp->zz), lp->nrows);
00410   CHECKRVALG (rval, CLEANUP);
00411   rval = mpf_ILLsvector_alloc (&(lp->yjz), lp->nrows);
00412   CHECKRVALG (rval, CLEANUP);
00413   rval = mpf_ILLsvector_alloc (&(lp->zA), lp->nnbasic);
00414   CHECKRVALG (rval, CLEANUP);
00415   rval = mpf_ILLsvector_alloc (&(lp->work), lp->ncols);
00416   CHECKRVALG (rval, CLEANUP);
00417   rval = mpf_ILLsvector_alloc (&(lp->srhs), lp->nrows);
00418   CHECKRVALG (rval, CLEANUP);
00419   rval = mpf_ILLsvector_alloc (&(lp->ssoln), lp->nrows);
00420   CHECKRVALG (rval, CLEANUP);
00421   ILL_SAFE_MALLOC (lp->iwork, lp->ncols, int);
00422 
00423   for (i = 0; i < lp->ncols; i++)
00424   {
00425     lp->work.indx[i] = 0;
00426     mpf_EGlpNumZero (lp->work.coef[i]);
00427     lp->iwork[i] = 0;
00428   }
00429   n = lp->nrows > lp->ncols ? 2 * (lp->nrows) + 1 : 2 * (lp->ncols) + 1;
00430   lp->upd.t = mpf_EGlpNumAllocArray (n);
00431   ILL_SAFE_MALLOC (lp->upd.perm, n, int);
00432   ILL_SAFE_MALLOC (lp->upd.ix, n, int);
00433 
00434 
00435   ILL_SAFE_MALLOC (lp->bfeas, lp->nrows, int);
00436   ILL_SAFE_MALLOC (lp->dfeas, lp->nnbasic, int);
00437 
00438   mpf_init_simplex_tols (lp);
00439   mpf_ILLfct_init_counts (lp);
00440 
00441   lp->nbchange = 0;
00442   lp->ncchange = 0;
00443 
00444   lp->pIratio = RATIOTEST_HARRIS;
00445   lp->pIIratio = RATIOTEST_HARRIS;
00446   lp->dIratio = RATIOTEST_HARRIS;
00447   lp->dIIratio = RATIOTEST_HARRIS;
00448   lp->starttime = ILLutil_zeit ();
00449   ILLutil_sprand (1, &(lp->rstate));
00450 
00451 CLEANUP:
00452   if (rval)
00453     mpf_free_internal_lpinfo (lp);
00454   EG_RETURN (rval);
00455 }
00456 
00457 int mpf_ILLsimplex_retest_psolution (
00458   mpf_lpinfo * lp,
00459   mpf_price_info * p,
00460   int phase,
00461   mpf_feas_info * fi)
00462 {
00463   int rval = 0;
00464   int fbid = lp->fbasisid;
00465   int bid = lp->basisid;
00466   mpf_t *ptol = &(lp->tol->pfeas_tol);
00467   mpf_t *dtol = &(lp->tol->dfeas_tol);
00468   mpf_t *iptol = &(lp->tol->ip_tol);
00469   mpf_t *idtol = &(lp->tol->id_tol);
00470 
00471   fi->pstatus = -1;
00472   fi->dstatus = -1;
00473   if (fbid < bid - PARAM_PRIMAL_REFACTORGAP)
00474   {
00475     rval = mpf_ILLbasis_refactor (lp);
00476     CHECKRVALG (rval, CLEANUP);
00477   }
00478   if (fbid < bid - PARAM_PRIMAL_RESOLVEGAP)
00479     mpf_ILLfct_compute_xbz (lp);
00480 
00481   if (phase == PRIMAL_PHASEII)
00482   {
00483     if (fbid < bid - PARAM_PRIMAL_RESOLVEGAP)
00484     {
00485       mpf_ILLfct_compute_piz (lp);
00486       mpf_ILLfct_compute_dz (lp);
00487       if (p != NULL && p->p_strategy == COMPLETE_PRICING)
00488         mpf_ILLprice_compute_dual_inf (lp, p, NULL, 0, PRIMAL_PHASEII);
00489     }
00490     mpf_ILLfct_compute_pobj (lp);
00491     mpf_ILLfct_check_pfeasible (lp, fi, *ptol);
00492     mpf_ILLfct_check_dfeasible (lp, fi, *dtol);
00493   }
00494   else if (phase == PRIMAL_PHASEI)
00495   {
00496     mpf_ILLfct_check_pfeasible (lp, fi, *iptol);
00497     if (fi->pstatus != PRIMAL_FEASIBLE)
00498     {
00499       if (lp->pIpiz)
00500       {
00501         mpf_ILLfct_compute_phaseI_piz (lp);
00502         mpf_ILLfct_compute_phaseI_dz (lp);
00503         mpf_ILLfct_check_pIdfeasible (lp, fi, *idtol);
00504         if (p != NULL && p->p_strategy == COMPLETE_PRICING)
00505           mpf_ILLprice_compute_dual_inf (lp, p, NULL, 0, PRIMAL_PHASEI);
00506       }
00507     }
00508   }
00509 CLEANUP:
00510   if (rval == QS_LP_CHANGE_PREC)
00511   {
00512     MESSAGE (__QS_SB_VERB, "Changing mpf_precision");
00513     return rval;
00514   }
00515   EG_RETURN (rval);
00516 }
00517 
00518 int mpf_ILLsimplex_retest_dsolution (
00519   mpf_lpinfo * lp,
00520   mpf_price_info * p,
00521   int phase,
00522   mpf_feas_info * fi)
00523 {
00524   int rval = 0;
00525   int fbid = lp->fbasisid;
00526   int bid = lp->basisid;
00527   mpf_t *ptol = &(lp->tol->pfeas_tol);
00528   mpf_t *dtol = &(lp->tol->dfeas_tol);
00529   mpf_t *iptol = &(lp->tol->ip_tol);
00530   mpf_t *idtol = &(lp->tol->id_tol);
00531 
00532   //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00533 
00534   fi->pstatus = -1;
00535   fi->dstatus = -1;
00536   if (fbid < bid - PARAM_DUAL_REFACTORGAP)
00537   {
00538     //ILL_IFTRACE("Refactor: %s:%s:%d\n",__func__,__FILE__,__LINE__);
00539     rval = mpf_ILLbasis_refactor (lp);
00540     CHECKRVALG (rval, CLEANUP);
00541   }
00542   if (fbid < bid - PARAM_DUAL_RESOLVEGAP)
00543   {
00544     //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00545     mpf_ILLfct_compute_piz (lp);
00546     mpf_ILLfct_compute_dz (lp);
00547   }
00548 
00549   if (phase == DUAL_PHASEII)
00550   {
00551     //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00552     if (fbid < bid - PARAM_DUAL_RESOLVEGAP)
00553     {
00554       //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00555       mpf_ILLfct_compute_xbz (lp);
00556       CHECKRVALG (rval, CLEANUP);
00557       if (p != NULL)
00558       {
00559         //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00560         if (p->d_strategy == COMPLETE_PRICING)
00561         {
00562           //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00563           mpf_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEII);
00564         }
00565         else
00566         {
00567           //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00568           mpf_ILLprice_update_mpartial_price (lp, p, DUAL_PHASEII, ROW_PRICING);
00569         }
00570       }
00571     }
00572     //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
00573     mpf_ILLfct_compute_dobj (lp);
00574     mpf_ILLfct_check_dfeasible (lp, fi, *dtol);
00575     mpf_ILLfct_check_pfeasible (lp, fi, *ptol);
00576   }
00577   else if (phase == DUAL_PHASEI)
00578   {
00579     mpf_ILLfct_check_dfeasible (lp, fi, *idtol);
00580     if (fi->dstatus != DUAL_FEASIBLE)
00581     {
00582       mpf_ILLfct_compute_phaseI_xbz (lp);
00583       mpf_ILLfct_check_pIpfeasible (lp, fi, *iptol);
00584       if (p != NULL)
00585       {
00586         if (p->d_strategy == COMPLETE_PRICING)
00587           mpf_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEI);
00588         else
00589           mpf_ILLprice_update_mpartial_price (lp, p, DUAL_PHASEI, ROW_PRICING);
00590       }
00591     }
00592   }
00593 CLEANUP:
00594   EG_RETURN (rval);
00595 }
00596 
00597 int mpf_ILLsimplex_solution (
00598   mpf_lpinfo * lp,
00599   mpf_t * xz,
00600   mpf_t * piz,
00601   mpf_t * dz,
00602   mpf_t * objval)
00603 {
00604   int i, j;
00605   int col;
00606 
00607   if (xz != NULL)
00608   {
00609     if (lp->basisstat.optimal == 0)
00610     {
00611       EG_RETURN (1);
00612     }
00613     for (i = 0; i < lp->nrows; i++)
00614       mpf_EGlpNumCopy (xz[lp->baz[i]], lp->xbz[i]);
00615     for (j = 0; j < lp->nnbasic; j++)
00616     {
00617       col = lp->nbaz[j];
00618       if (lp->vstat[col] == STAT_UPPER)
00619         mpf_EGlpNumCopy (xz[col], lp->uz[col]);
00620       else if (lp->vstat[col] == STAT_LOWER)
00621         mpf_EGlpNumCopy (xz[col], lp->lz[col]);
00622       else
00623         mpf_EGlpNumZero (xz[col]);
00624     }
00625   }
00626   if (piz != NULL)
00627   {
00628     if (lp->basisstat.optimal == 0)
00629     {
00630       EG_RETURN (1);
00631     }
00632     for (i = 0; i < lp->nrows; i++)
00633       mpf_EGlpNumCopy (piz[i], lp->piz[i]);
00634   }
00635   if (dz != NULL)
00636   {
00637     if (lp->basisstat.optimal == 0)
00638     {
00639       EG_RETURN (1);
00640     }
00641     for (i = 0; i < lp->nrows; i++)
00642       mpf_EGlpNumZero (dz[lp->baz[i]]);
00643     for (j = 0; j < lp->nnbasic; j++)
00644       mpf_EGlpNumCopy (dz[lp->nbaz[j]], lp->dz[j]);
00645   }
00646   if (objval != NULL)
00647     mpf_EGlpNumCopy (*objval, lp->objval);
00648   return 0;
00649 }
00650 
00651 int mpf_ILLsimplex_infcertificate (
00652   mpf_lpinfo * lp,
00653   mpf_t * pi)
00654 {
00655   int i, col, nz;
00656   char *sense;
00657   mpf_t *x, *l, *u;
00658   mpf_lp_status_info *ls;
00659 
00660   if (pi == NULL)
00661     return 0;
00662 
00663   ls = &(lp->basisstat);
00664   if (ls->primal_infeasible == 0 && ls->dual_unbounded == 0)
00665   {
00666     EG_RETURN (1);
00667   }
00668 
00669   if (lp->final_phase == PRIMAL_PHASEI && lp->pIpiz != NULL)
00670   {
00671     for (i = 0; i < lp->nrows; i++)
00672       mpf_EGlpNumCopy (pi[i], lp->pIpiz[i]);
00673   }
00674   else if (lp->final_phase == DUAL_PHASEII && lp->infub_ix != -1)
00675   {
00676     col = lp->baz[lp->infub_ix];
00677     x = &(lp->xbz[lp->infub_ix]);
00678     l = &(lp->lz[col]);
00679     u = &(lp->uz[col]);
00680 
00681     for (i = 0; i < lp->nrows; i++)
00682       mpf_EGlpNumZero (pi[i]);
00683 
00684     if (mpf_EGlpNumIsNeqq (*l, mpf_NINFTY) && mpf_EGlpNumIsLess (*x, *l))
00685     {
00686       for (i = 0, nz = lp->zz.nzcnt; i < nz; i++)
00687         mpf_EGlpNumCopyNeg (pi[lp->zz.indx[i]], lp->zz.coef[i]);
00688     }
00689     else
00690     {
00691       for (i = 0, nz = lp->zz.nzcnt; i < nz; i++)
00692         mpf_EGlpNumCopy (pi[lp->zz.indx[i]], lp->zz.coef[i]);
00693     }
00694   }
00695   else
00696   {
00697     fprintf (stderr, "Invalid call to inf. certificate routine\n");
00698     EG_RETURN (1);
00699   }
00700 
00701   sense = lp->O->sense;
00702   for (i = 0; i < lp->nrows; i++)
00703   {
00704     if (sense[i] == 'G' && mpf_EGlpNumIsLessZero (pi[i]))
00705       mpf_EGlpNumZero (pi[i]);
00706     if (sense[i] == 'L' && mpf_EGlpNumIsGreatZero (pi[i]))
00707       mpf_EGlpNumZero (pi[i]);
00708   }
00709   return 0;
00710 }
00711 
00712 #if SIMPLEX_DEBUG > 1
00713 static void mpf_test_cert (
00714   mpf_lpinfo * lp,
00715   mpf_t * pi)
00716 {
00717   int i, j;
00718   int mcnt, mbeg;
00719   mpf_t fsum, sum;
00720 
00721   mpf_EGlpNumInitVar (fsum);
00722   mpf_EGlpNumInitVar (sum);
00723   mpf_EGlpNumZero (fsum);
00724 
00725   for (i = 0; i < lp->nrows; i++)
00726   {
00727     if (lp->O->sense[i] == 'G' && mpf_EGlpNumIsLessZero (pi[i]))
00728       printf ("compl \n");
00729     if (lp->O->sense[i] == 'L' && mpf_EGlpNumIsGreatZero (pi[i]))
00730       printf ("compll \n");
00731   }
00732 
00733   for (i = 0; i < lp->nrows; i++)
00734     mpf_EGlpNumAddInnProdTo (fsum, pi[i], lp->bz[i]);
00735 
00736   for (j = 0; j < lp->nnbasic; j++)
00737   {
00738     mpf_EGlpNumZero (sum);
00739     mcnt = lp->matcnt[j];
00740     mbeg = lp->matbeg[j];
00741     for (i = 0; i < mcnt; i++)
00742       mpf_EGlpNumAddInnProdTo (sum, pi[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00743 
00744     if (mpf_EGlpNumIsLess (mpf_PFEAS_TOLER, sum) &&
00745         (lp->vtype[j] == VLOWER || lp->vtype[j] == VFREE))
00746       printf ("compl2\n");
00747     else
00748     {
00749       mpf_EGlpNumSign (sum);
00750       if (mpf_EGlpNumIsLess (mpf_PFEAS_TOLER, sum) &&
00751           (lp->vtype[j] == VUPPER || lp->vtype[j] == VFREE))
00752         printf ("compl1\n");
00753       mpf_EGlpNumSign (sum);
00754     }
00755 
00756     if (mpf_EGlpNumIsLessZero (sum)
00757         && (lp->vtype[j] & (VFREE | VUPPER)) == 0)
00758       mpf_EGlpNumSubInnProdTo (fsum, sum, lp->lz[j]);
00759     else if (mpf_EGlpNumIsGreatZero (sum)
00760              && (lp->vtype[j] & (VFREE | VLOWER)) == 0)
00761       mpf_EGlpNumSubInnProdTo (fsum, sum, lp->uz[j]);
00762   }
00763   printf ("fsum = %.8f\n", mpf_EGlpNumToLf (fsum));
00764   mpf_EGlpNumClearVar (fsum);
00765   mpf_EGlpNumClearVar (sum);
00766 }
00767 #endif
00768 
00769 static void mpf_save_paraminfo (
00770   mpf_price_info * pinf,
00771   mpf_iter_info * it)
00772 {
00773   mpf_param_info *pr = &(it->oldinfo);
00774 
00775   pr->origalgo = it->algorithm;
00776   pr->pphaseI = pinf->pI_price;
00777   pr->pphaseII = pinf->pII_price;
00778   pr->dphaseI = pinf->dI_price;
00779   pr->dphaseII = pinf->dII_price;
00780   pr->p_strategy = pinf->p_strategy;
00781   pr->d_strategy = pinf->d_strategy;
00782 }
00783 
00784 static void mpf_restore_paraminfo (
00785   mpf_iter_info * it,
00786   mpf_price_info * pinf)
00787 {
00788   mpf_param_info *pr = &(it->oldinfo);
00789 
00790   it->algorithm = pr->origalgo;
00791   pinf->pI_price = pr->pphaseI;
00792   pinf->pII_price = pr->pphaseII;
00793   pinf->dI_price = pr->dphaseI;
00794   pinf->dII_price = pr->dphaseII;
00795   pinf->p_strategy = pr->p_strategy;
00796   pinf->d_strategy = pr->d_strategy;
00797 }
00798 
00799 int mpf_ILLsimplex (
00800   mpf_lpinfo * lp,
00801   int algorithm,
00802   mpf_ILLlp_basis * B,
00803   mpf_price_info * pinf,
00804   int *status,
00805   int sdisplay,
00806   itcnt_t*itcnt)
00807 {
00808   int phase = -1;
00809   int singular = -1;
00810   int rval = 0;
00811   int new_price = -1;
00812   mpf_svector wz;
00813   mpf_svector updz;
00814   mpf_feas_info fi;
00815   mpf_iter_info it;
00816 
00817   mpf_EGlpNumInitVar (fi.totinfeas);
00818   mpf_EGlpNumInitVar (it.prevobj);
00819   mpf_EGlpNumInitVar (it.objtol);
00820 
00821   it.newphase = -1;
00822   it.nextphase = -1;
00823   it.nextstep = -1;
00824   it.sdisplay = sdisplay;
00825   it.n_pivot_fail = 0;
00826   it.itercnt = 0;
00827   it.n_restart = 0;
00828   it.solstatus = ILL_LP_UNSOLVED;
00829   it.curtime = 0;
00830   it.rounds = 0;
00831   mpf_EGlpNumCopy (it.prevobj, mpf_INFTY);
00832   it.nosolve = 0;
00833   it.noprog = 0;
00834   mpf_EGlpNumCopy (it.objtol, mpf_OBJBND_TOLER);
00835   it.chkobj = PARAM_MAX_NOPROG;
00836   it.inner = 0;
00837   it.algorithm = algorithm;
00838   it.pricetype = -1;
00839   it.resumeid = -1;
00840   mpf_save_paraminfo (pinf, &it);
00841 
00842 #if SIMPLEX_DEBUG > 0
00843   if (lp->O->nrows > 1000)
00844     it.sdisplay = 1;
00845 #endif
00846   if (status)
00847     *status = QS_LP_UNSOLVED;
00848 
00849   mpf_free_internal_lpinfo (lp);
00850   mpf_init_internal_lpinfo (lp);
00851   rval = mpf_build_internal_lpinfo (lp);
00852   CHECKRVALG (rval, CLEANUP);
00853 
00854   mpf_ILLsvector_init (&wz);
00855   rval = mpf_ILLsvector_alloc (&wz, lp->nrows);
00856   CHECKRVALG (rval, CLEANUP);
00857   mpf_ILLsvector_init (&updz);
00858   rval = mpf_ILLsvector_alloc (&updz, lp->nrows);
00859   CHECKRVALG (rval, CLEANUP);
00860 
00861   if (it.sdisplay)
00862   {
00863     char buffer[256];
00864     int nonzero = 0;
00865     register int i = lp->ncols;
00866 
00867     while (i--)
00868       nonzero += lp->matcnt[i];
00869     sprintf (buffer, "starting mpf_ILLsimplex on %s...\n", lp->O->probname);
00870     /* depending on LP's reporter 
00871      * string is printed to stdout 
00872      * or handed to GUI */
00873     rval = rval || ILLstring_report (buffer, &lp->O->reporter);
00874     printf ("Problem has %d rows and %d cols and %d nonzeros\n", lp->nrows,
00875             lp->ncols, nonzero);
00876     fflush (stdout);
00877   }
00878   mpf_ILLfct_set_variable_type (lp);
00879 
00880   if (B != 0)
00881   {
00882     rval = mpf_ILLbasis_load (lp, B);
00883     CHECKRVALG (rval, CLEANUP);
00884     if (it.algorithm == DUAL_SIMPLEX)
00885     {
00886       if (B->rownorms)
00887       {
00888         rval = mpf_ILLprice_load_rownorms (lp, B->rownorms, pinf);
00889         CHECKRVALG (rval, CLEANUP);
00890       }
00891       else
00892         mpf_EGlpNumFreeArray (pinf->dsinfo.norms);
00893     }
00894     else if (it.algorithm == PRIMAL_SIMPLEX)
00895     {
00896       if (B->colnorms)
00897       {
00898         rval = mpf_ILLprice_load_colnorms (lp, B->colnorms, pinf);
00899         CHECKRVALG (rval, CLEANUP);
00900       }
00901       else
00902         mpf_EGlpNumFreeArray (pinf->psinfo.norms);
00903     }
00904     else if (it.algorithm != PRIMAL_OR_DUAL)
00905     {
00906       fprintf (stderr, "Unknown algorithm %d in mpf_ILLsimplex\n", it.algorithm);
00907       rval = 1;
00908       ILL_CLEANUP;
00909     }
00910   }
00911   else if (lp->basisid == -1)
00912   {
00913     if (lp->nrows < 200 && lp->ncols < 400)
00914       rval = mpf_ILLbasis_get_initial (lp, it.algorithm);
00915     else
00916       rval = mpf_ILLbasis_get_cinitial (lp, it.algorithm);
00917     CHECKRVALG (rval, CLEANUP);
00918     mpf_ILLprice_free_pricing_info (pinf);
00919   }
00920 
00921   if (lp->fbasisid != lp->basisid)
00922   {
00923     rval = mpf_ILLbasis_factor (lp, &singular);
00924     CHECKRVALG (rval, CLEANUP);
00925     if (singular)
00926     {
00927       MESSAGE (__QS_SB_VERB, "Singular basis found!");
00928       mpf_ILLprice_free_pricing_info (pinf);
00929     }
00930   }
00931 
00932 START:
00933 #if 0
00934   if (it.resumeid == SIMPLEX_RESUME_UNSHIFT)
00935     fprintf (stderr, "Resuming Unshift\n");
00936   else if (it.resumeid == SIMPLEX_RESUME_SING)
00937     fprintf (stderr, "Resuming Singular\n");
00938   else if (it.resumeid == SIMPLEX_RESUME_NUMER)
00939     fprintf (stderr, "Resuming Numer\n");
00940   else if (it.resumeid != -1)
00941     fprintf (stderr, "Resuming for other reason... %d\n", it.resumeid);
00942 #endif
00943   it.solstatus = ILL_LP_UNSOLVED;
00944   mpf_init_lp_status_info (&(lp->basisstat));
00945 
00946   mpf_ILLfct_compute_piz (lp);
00947   mpf_ILLfct_compute_dz (lp);
00948   if (it.algorithm == DUAL_SIMPLEX)
00949   {
00950     if (B != NULL || it.resumeid == SIMPLEX_RESUME_UNSHIFT)
00951       mpf_ILLfct_dual_adjust (lp, lp->tol->dfeas_tol);
00952     else
00953       mpf_ILLfct_dual_adjust (lp, mpf_zeroLpNum);
00954   }
00955   mpf_ILLfct_compute_xbz (lp);
00956 
00957   mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
00958   mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
00959   mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
00960 
00961   if (fi.dstatus == DUAL_FEASIBLE && mpf_EGlpNumIsNeqq (lp->objbound, mpf_INFTY))
00962   {
00963     mpf_ILLfct_compute_dobj (lp);
00964     if (mpf_EGlpNumIsLess (lp->objbound, lp->dobjval))
00965     {
00966       it.solstatus = ILL_BND_REACHED;
00967       printf ("solstatus = ILL_BND_REACHED = 5 %lf %lf\n",
00968               mpf_EGlpNumToLf (lp->objbound), mpf_EGlpNumToLf (lp->dobjval));
00969       goto TERMINATE;
00970     }
00971   }
00972   if (fi.pstatus == PRIMAL_FEASIBLE && fi.dstatus == DUAL_FEASIBLE)
00973   {
00974     it.solstatus = ILL_LP_SOLVED;
00975     mpf_ILLfct_compute_pobj (lp);
00976     goto TERMINATE;
00977   }
00978 
00979   if (it.algorithm == PRIMAL_OR_DUAL)
00980   {
00981     if (fi.pstatus == PRIMAL_FEASIBLE)
00982       it.algorithm = PRIMAL_SIMPLEX;
00983     else if (fi.dstatus == DUAL_FEASIBLE)
00984       it.algorithm = DUAL_SIMPLEX;
00985     else if (mpf_EGlpNumToLf (lp->pinfeas) < 10 * mpf_EGlpNumToLf (lp->dinfeas))
00986       it.algorithm = PRIMAL_SIMPLEX;
00987     else
00988       it.algorithm = DUAL_SIMPLEX;
00989   }
00990 
00991   if (it.algorithm == PRIMAL_SIMPLEX)
00992   {
00993     if (fi.pstatus == PRIMAL_FEASIBLE)
00994       phase = PRIMAL_PHASEII;
00995     else
00996       phase = PRIMAL_PHASEI;
00997   }
00998   else if (it.algorithm == DUAL_SIMPLEX)
00999   {
01000     if (fi.dstatus == DUAL_FEASIBLE)
01001       phase = DUAL_PHASEII;
01002     else
01003       phase = DUAL_PHASEI;
01004   }
01005 
01006   rval = mpf_ILLprice_build_pricing_info (lp, pinf, phase);
01007   CHECKRVALG (rval, CLEANUP);
01008 
01009   it.newphase = SIMPLEX_PHASE_NEW;
01010   it.nextstep = SIMPLEX_CONTINUE;
01011 
01012   while (it.nextstep == SIMPLEX_CONTINUE)
01013   {
01014     //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
01015 
01016     if (phase == PRIMAL_PHASEI)
01017     {
01018       rval = mpf_primal_phaseI_step (lp, pinf, &updz, &wz, &it);
01019       CHECKRVALG (rval, CLEANUP);
01020     }
01021 
01022     else if (phase == PRIMAL_PHASEII)
01023     {
01024       rval = mpf_primal_phaseII_step (lp, pinf, &updz, &wz, &it);
01025       CHECKRVALG (rval, CLEANUP);
01026     }
01027 
01028     else if (phase == DUAL_PHASEI)
01029     {
01030       rval = mpf_dual_phaseI_step (lp, pinf, &updz, &wz, &it);
01031       CHECKRVALG (rval, CLEANUP);
01032     }
01033 
01034     else if (phase == DUAL_PHASEII)
01035     {
01036       rval = mpf_dual_phaseII_step (lp, pinf, &updz, &wz, &it);
01037       CHECKRVALG (rval, CLEANUP);
01038     }
01039     if (it.nextstep == SIMPLEX_RESUME)
01040     {
01041       //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
01042       mpf_ILLprice_free_pricing_info (pinf);
01043       if (it.resumeid == SIMPLEX_RESUME_UNSHIFT)
01044       {
01045         if (it.pricetype == QS_PRICE_PDEVEX)
01046         {
01047           pinf->pI_price = QS_PRICE_PDEVEX;
01048           pinf->pII_price = QS_PRICE_PDEVEX;
01049         }
01050         else if (it.pricetype == QS_PRICE_DDEVEX)
01051         {
01052           pinf->dI_price = QS_PRICE_DDEVEX;
01053           pinf->dII_price = QS_PRICE_DDEVEX;
01054         }
01055       }
01056       else if (it.resumeid == SIMPLEX_RESUME_NUMER)
01057       {
01058         mpf_ILLfct_unroll_bound_change (lp);
01059         mpf_ILLfct_unroll_coef_change (lp);
01060         /* we are disabling re-do under this circunstances ! */
01061         rval = mpf_ILLbasis_get_initial (lp, it.algorithm);
01062         CHECKRVALG (rval, CLEANUP);
01063         rval = mpf_ILLbasis_factor (lp, &singular);
01064         if (singular)
01065           MESSAGE (__QS_SB_VERB, "Singular basis found!");
01066         CHECKRVALG (rval, CLEANUP);
01067       }
01068       it.pricetype = -1;
01069       if (it.n_restart > SIMPLEX_MAX_RESTART)
01070       {
01071         it.solstatus = ILL_MAX_ITER;
01072         goto LIMIT_TERMINATE;
01073       }
01074       goto START;
01075     }
01076     else if (it.nextstep == SIMPLEX_CONTINUE)
01077     {
01078       //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
01079       it.itercnt++;
01080 
01081       if (it.nextphase != phase)
01082       {
01083         //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
01084         it.newphase = SIMPLEX_PHASE_NEW;
01085         phase = it.nextphase;
01086         new_price = mpf_ILLprice_get_price (pinf, phase);
01087 
01088         if (pinf->cur_price != new_price)
01089         {
01090           //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
01091           mpf_ILLprice_free_pricing_info (pinf);
01092           rval = mpf_ILLprice_build_pricing_info (lp, pinf, phase);
01093           CHECKRVALG (rval, CLEANUP);
01094         }
01095       }
01096     }
01097   }
01098 
01099 #if SIMPLEX_DEBUG > 0
01100   mpf_ILLfct_print_counts (lp);
01101 #endif
01102 LIMIT_TERMINATE:
01103   rval = mpf_terminate_simplex (lp, phase, &it);
01104   CHECKRVALG (rval, CLEANUP);
01105 
01106 TERMINATE:
01107   mpf_restore_paraminfo (&it, pinf);
01108 
01109   if (it.sdisplay)
01110   {
01111     printf ("completed mpf_ILLsimplex\n");
01112     printf ("%s: ", lp->O->probname);
01113     fflush (stdout);
01114   }
01115 
01116   if (status)
01117   {
01118     if (it.solstatus == ILL_MAX_ITER)
01119     {
01120       *status = QS_LP_ITER_LIMIT;
01121     }
01122     else if (it.solstatus == ILL_MAX_TIME)
01123     {
01124       *status = QS_LP_TIME_LIMIT;
01125     }
01126     else if (it.solstatus == ILL_LP_ABORTED)
01127     {
01128       *status = QS_LP_ABORTED;
01129     }
01130     else if (it.solstatus == ILL_PPHASEI_ERROR ||
01131              it.solstatus == ILL_PPHASEII_ERROR ||
01132              it.solstatus == ILL_DPHASEI_ERROR ||
01133              it.solstatus == ILL_DPHASEII_ERROR)
01134     {
01135       *status = QS_LP_NUMERR;
01136     }
01137     else if(it.solstatus == ILL_LP_UNSOLVED)
01138     {
01139       *status = QS_LP_UNSOLVED;
01140     }
01141     else if (it.solstatus == ILL_BND_REACHED)
01142     {
01143       *status = QS_LP_OBJ_LIMIT;
01144     }
01145     else if (it.solstatus == ILL_LP_SOLVED)
01146     {
01147       if (lp->basisstat.optimal)
01148       {
01149         *status = QS_LP_OPTIMAL;
01150       }
01151       else if (lp->basisstat.primal_infeasible || lp->basisstat.dual_unbounded)
01152       {
01153         *status = QS_LP_INFEASIBLE;
01154         if (it.sdisplay)
01155         {
01156           if (lp->basisstat.primal_infeasible)
01157             fprintf (stdout, "Primal Infeasible\n");
01158           else
01159             fprintf (stdout, "Dual Unbounded\n");
01160         }
01161       }
01162       else if (lp->basisstat.primal_unbounded)
01163       {
01164         *status = QS_LP_UNBOUNDED;
01165       }
01166     }
01167     else
01168     {
01169       fprintf (stderr, "unknown solution status in mpf_ILLsimplex %d\n",
01170                it.solstatus);
01171       rval = 1;
01172       CHECKRVALG (rval, CLEANUP);
01173     }
01174   }
01175 
01176 #if SIMPLEX_DEBUG > 1
01177   {
01178     int rva = 0;
01179     mpf_t *pi = NULL;
01180 
01181     pi = mpf_EGlpNumAllocArray (lp->nrows);
01182     rva = mpf_ILLsimplex_infcertificate (lp, pi);
01183     printf ("rva = %d\n", rva);
01184     if (!rva)
01185     {
01186       mpf_test_cert (lp, pi);
01187     }
01188     mpf_EGlpNumFreeArray (pi);
01189   }
01190 #endif
01191    /* update counter */
01192    itcnt->pI_iter += lp->cnts->pI_iter;
01193    itcnt->pII_iter += lp->cnts->pII_iter;
01194    itcnt->dI_iter += lp->cnts->dI_iter;
01195    itcnt->dII_iter += lp->cnts->dII_iter;
01196    itcnt->tot_iter = itcnt->pI_iter + itcnt->pII_iter + itcnt->dI_iter +
01197                       itcnt->dII_iter;
01198    /* end update */
01199   if (it.sdisplay)
01200   {
01201     int bstat = 0;
01202 
01203     printf ("time = %.3f, pI = %d, pII = %d, dI = %d, dII = %d, ",
01204             ILLutil_zeit () - lp->starttime, lp->cnts->pI_iter,
01205             lp->cnts->pII_iter, lp->cnts->dI_iter, lp->cnts->dII_iter);
01206     fflush (stdout);
01207     mpf_get_current_stat (&(lp->basisstat), it.algorithm, &bstat);
01208     switch (bstat)
01209     {
01210     case OPTIMAL:
01211       printf ("opt = %f\n", mpf_EGlpNumToLf (lp->objval));
01212       break;
01213     case PRIMAL_INFEASIBLE:
01214       printf ("no primal soln\n");
01215       break;
01216     case PRIMAL_UNBOUNDED:
01217       printf ("primal unbounded\n");
01218       break;
01219     case PRIMAL_FEASIBLE:
01220       printf ("primal obj = %f\n", mpf_EGlpNumToLf (lp->pobjval));
01221       break;
01222     case DUAL_INFEASIBLE:
01223       printf ("no dual soln\n");
01224       break;
01225     case DUAL_UNBOUNDED:
01226       printf ("dual unbounded\n");
01227       break;
01228     case DUAL_FEASIBLE:
01229       printf ("dual obj = %f\n", mpf_EGlpNumToLf (lp->dobjval));
01230       break;
01231     }
01232     fflush (stdout);
01233 
01234     if (it.sdisplay > 1)
01235     {
01236       if (it.algorithm == PRIMAL_SIMPLEX && pinf->pI_price == QS_PRICE_PDEVEX)
01237         printf ("Devex norms initialised %d times\n", pinf->pdinfo.ninit);
01238       fflush (stdout);
01239     }
01240   }
01241 
01242 CLEANUP:
01243   mpf_ILLsvector_free (&wz);
01244   mpf_ILLsvector_free (&updz);
01245   mpf_EGlpNumClearVar (it.prevobj);
01246   mpf_EGlpNumClearVar (it.objtol);
01247   mpf_EGlpNumClearVar (fi.totinfeas);
01248   if (rval == QS_LP_CHANGE_PREC)
01249   {
01250     MESSAGE (__QS_SB_VERB, "Changing mpf_precision");
01251     return rval;
01252   }
01253   else
01254   {
01255     MESSAGE (rval ? 0 : 1000, "Error code %d", rval);
01256     EG_RETURN (rval);
01257   }
01258 }
01259 
01260 static int mpf_terminate_simplex (
01261   mpf_lpinfo * lp,
01262   int phase,
01263   mpf_iter_info * it)
01264 {
01265   int rval = 0;
01266   int sphase;
01267   mpf_feas_info fi;
01268 
01269   mpf_EGlpNumInitVar (fi.totinfeas);
01270 
01271   if (it->solstatus != ILL_MAX_TIME && it->solstatus != ILL_MAX_ITER)
01272     ILL_CLEANUP;
01273 
01274   if (it->algorithm == PRIMAL_SIMPLEX)
01275   {
01276     if (lp->nbchange != 0)
01277     {
01278       if (it->sdisplay > 1)
01279       {
01280         printf ("unrolling %d bound shifts\n", lp->nbchange);
01281         fflush (stdout);
01282       }
01283       mpf_ILLfct_unroll_bound_change (lp);
01284     }
01285     rval = mpf_ILLsimplex_retest_psolution (lp, NULL, phase, &fi);
01286     CHECKRVALG (rval, CLEANUP);
01287 
01288     sphase = (phase == PRIMAL_PHASEI) ? PHASEI : PHASEII;
01289     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, sphase);
01290   }
01291   else if (it->algorithm == DUAL_SIMPLEX)
01292   {
01293     if (lp->ncchange != 0)
01294     {
01295       if (it->sdisplay > 1)
01296       {
01297         printf ("unrolling %d coef shifts\n", lp->ncchange);
01298         fflush (stdout);
01299       }
01300       mpf_ILLfct_unroll_coef_change (lp);
01301     }
01302     rval = mpf_ILLsimplex_retest_dsolution (lp, NULL, phase, &fi);
01303     CHECKRVALG (rval, CLEANUP);
01304 
01305     sphase = (phase == DUAL_PHASEI) ? PHASEI : PHASEII;
01306     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, sphase, PHASEII);
01307   }
01308 
01309 CLEANUP:
01310   mpf_EGlpNumClearVar (fi.totinfeas);
01311   EG_RETURN (rval);
01312 }
01313 
01314 static int mpf_test_progress (
01315   mpf_t objval,
01316   mpf_t prevobj)
01317 {
01318   mpf_t denom;
01319 
01320   mpf_EGlpNumInitVar (denom);
01321   mpf_EGlpNumCopyDiff (denom, objval, prevobj);
01322   if (mpf_EGlpNumIsNeqZero (objval, mpf_PROGRESS_ZERO))
01323     mpf_EGlpNumDivTo (denom, objval);
01324   if (!mpf_EGlpNumIsNeqZero (denom, mpf_PROGRESS_THRESH))
01325   {
01326     mpf_EGlpNumClearVar (denom);
01327     return 0;
01328   }
01329   else
01330   {
01331     mpf_EGlpNumClearVar (denom);
01332     return 1;
01333   }
01334 }
01335 
01336 static void mpf_monitor_iter (
01337   mpf_lpinfo * lp,
01338   mpf_price_info * p,
01339   mpf_iter_info * it,
01340   int phase)
01341 {
01342   mpf_t print_val;
01343   double tottime = ILLutil_zeit () - lp->starttime;
01344   int curtime = mpf_ILLutil_our_floor (tottime);  /* MONIKA */
01345   char print_str[20];
01346   mpf_feas_info fi;
01347   int aborted = 0;
01348 
01349   mpf_EGlpNumInitVar (print_val);
01350   mpf_EGlpNumInitVar (fi.totinfeas);
01351   mpf_EGlpNumZero (fi.totinfeas);
01352   mpf_EGlpNumZero (print_val);
01353 
01354   /* one of the following two time display mechanisms */
01355   switch (phase)
01356   {
01357   case PRIMAL_PHASEI:
01358     mpf_EGlpNumCopy (print_val, lp->pinfeas);
01359     mpf_EGlpNumAddTo (fi.totinfeas, lp->pinfeas);
01360     strcpy (print_str, "primal infeas");
01361     if (mpf_EGlpNumIsLessZero (lp->pinfeas) &&
01362         (mpf_EGlpNumIsNeqZero (lp->pinfeas, mpf_oneLpNum)))
01363     {
01364       /*printf ("Negative Infeasibility! Imposible %lg %la, iter %d\n",
01365        * mpf_EGlpNumToLf (print_val), mpf_EGlpNumToLf (print_val), it->itercnt);
01366        */
01367       //exit(1);
01368     }
01369     break;
01370   case PRIMAL_PHASEII:
01371     mpf_EGlpNumCopy (print_val, lp->pobjval);
01372     strcpy (print_str, "primal objval");
01373     break;
01374   case DUAL_PHASEI:
01375     mpf_EGlpNumCopy (print_val, lp->dinfeas);
01376     mpf_EGlpNumAddTo (fi.totinfeas, lp->dinfeas);
01377     strcpy (print_str, "dual infeas");
01378     break;
01379   case DUAL_PHASEII:
01380     mpf_EGlpNumCopy (print_val, lp->dobjval);
01381     strcpy (print_str, "dual objval");
01382     break;
01383   }
01384 
01385   aborted = mpf_report_value (lp, it, print_str, print_val);
01386   /*if (it->sdisplay && it->itercnt % lp->iterskip == 0) {
01387    * // printf ("(%d): %s = %f\n", it->itercnt, print_str, print_val);
01388    * // fflush (stdout);
01389    * } */
01390   if (curtime != it->curtime)
01391   {
01392     it->curtime = curtime;
01393     /*
01394      * if (it->sdisplay){
01395      * printf ("time = %d.0, ", curtime);
01396      * printf ("(%d): %s = %f\n", it->itercnt, print_str, print_val);
01397      * fflush (stdout);
01398      * }
01399      */
01400   }
01401 
01402   mpf_EGlpNumAddUiTo (fi.totinfeas, 1000);
01403   if (mpf_EGlpNumIsLessZero (fi.totinfeas))
01404   {
01405     it->nextstep = SIMPLEX_TERMINATE;
01406     it->solstatus = ILL_MAX_ITER;
01407     MESSAGE (it->sdisplay ? 0 : __QS_SB_VERB,
01408              "early finish by excess infeasibility");
01409     ILL_CLEANUP;
01410   }
01411 
01412   if (phase == DUAL_PHASEII && mpf_EGlpNumIsNeqq (lp->objbound, mpf_INFTY))
01413   {
01414     /*if (lp->dobjval > lp->objbound + it->objtol) */
01415     mpf_EGlpNumCopyDiff (print_val, lp->dobjval, lp->objbound);
01416     if (mpf_EGlpNumIsLess (it->objtol, print_val))
01417     {
01418       mpf_ILLfct_unroll_coef_change (lp);
01419       mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01420       mpf_ILLfct_set_status_values (lp, -1, fi.dstatus, -1, PHASEII);
01421 
01422       if (fi.dstatus == DUAL_FEASIBLE)
01423       {
01424         mpf_ILLfct_compute_dobj (lp);
01425         if (mpf_EGlpNumIsLess (lp->objbound, lp->dobjval))
01426         {
01427           it->solstatus = ILL_BND_REACHED;
01428           it->nextstep = SIMPLEX_TERMINATE;
01429           /*if (it->sdisplay) */
01430           {
01431             printf ("bound reached %lf %lf\n", mpf_EGlpNumToLf (lp->objbound),
01432                     mpf_EGlpNumToLf (lp->dobjval));
01433             fflush (stdout);
01434           }
01435         }
01436         else
01437           mpf_EGlpNumMultUiTo (it->objtol, 10);
01438       }
01439       else
01440       {
01441         it->nextphase = DUAL_PHASEI;
01442         it->newphase = SIMPLEX_PHASE_NEW;
01443         mpf_EGlpNumMultUiTo (it->objtol, 5);
01444       }
01445     }
01446   }
01447   if (it->itercnt >= lp->maxiter)
01448   {
01449     it->solstatus = ILL_MAX_ITER;
01450     it->nextstep = SIMPLEX_TERMINATE;
01451     if (it->sdisplay)
01452     {
01453       printf ("iter limit reached\n");
01454       fflush (stdout);
01455     }
01456     ILL_CLEANUP;
01457   }
01458   else if (tottime >= lp->maxtime)
01459   {
01460     it->solstatus = ILL_MAX_TIME;
01461     it->nextstep = SIMPLEX_TERMINATE;
01462     if (it->sdisplay)
01463     {
01464       printf ("time limit reached\n");
01465       fflush (stdout);
01466     }
01467     ILL_CLEANUP;
01468   }
01469   else if (aborted)
01470   {
01471     it->solstatus = ILL_LP_ABORTED;
01472     it->nextstep = SIMPLEX_TERMINATE;
01473     if (it->sdisplay)
01474     {
01475       printf ("aborted\n");
01476       fflush (stdout);
01477     }
01478     ILL_CLEANUP;
01479   }
01480   /* why is this commented out? */
01481   if(0){
01482     if (it->rounds && it->inner){
01483       it->inner --;
01484       if (it->inner == 0){
01485         printf ("restoring ..\n");
01486         mpf_restore_paraminfo (it, p);
01487         it->newphase   = SIMPLEX_PHASE_NEW;
01488         it->nextstep   = SIMPLEX_RESUME;
01489         /*it->resumeid   = SIMPLEX_RESUME_OUTER;*/
01490         ILL_CLEANUP;
01491       }
01492     }
01493   }
01494   if (phase == DUAL_PHASEII)
01495   {
01496     if (it->noprog > it->chkobj)
01497     {
01498       mpf_ILLfct_perturb_coefs (lp);
01499       it->noprog = 0;
01500       mpf_EGlpNumCopy (it->prevobj, lp->dobjval);
01501     }
01502   }
01503   else if (phase == PRIMAL_PHASEII)
01504   {
01505     if (it->noprog > it->chkobj)
01506     {
01507       mpf_ILLfct_perturb_bounds (lp);
01508       it->noprog = 0;
01509       mpf_EGlpNumCopy (it->prevobj, lp->pobjval);
01510     }
01511   }
01512   else if (phase == PRIMAL_PHASEI)
01513   {
01514     if (it->noprog > it->chkobj)
01515     {
01516       it->algorithm = DUAL_SIMPLEX;
01517       it->nextstep = SIMPLEX_RESUME;
01518       it->resumeid = SIMPLEX_RESUME_NUMER;
01519       /* this is to force to exit in the case of bad basis */
01520       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01521       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01522       it->n_restart++;
01523       //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
01524     }
01525   }
01526   else if (phase == DUAL_PHASEI)
01527   {
01528     if (it->noprog > it->chkobj)
01529     {
01530       it->algorithm = PRIMAL_SIMPLEX;
01531       it->nextstep = SIMPLEX_RESUME;
01532       it->resumeid = SIMPLEX_RESUME_NUMER;
01533       /* this is to force to exit in the case of bad basis */
01534       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01535       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01536       it->n_restart++;
01537       //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
01538     }
01539   }
01540 CLEANUP:
01541   mpf_EGlpNumClearVar (fi.totinfeas);
01542   mpf_EGlpNumClearVar (print_val);
01543   return;
01544 }
01545 
01546 static int mpf_primal_phaseI_step (
01547   mpf_lpinfo * lp,
01548   mpf_price_info * pinf,
01549   mpf_svector * updz,
01550   mpf_svector * wz,
01551   mpf_iter_info * it)
01552 {
01553   int rval = 0;
01554   int singular = 0;
01555   int refactor = 0;
01556   int cphase = PRIMAL_PHASEI;
01557   mpf_t alpha;
01558   mpf_feas_info fi;
01559   mpf_ratio_res rs;
01560   mpf_price_res pr;
01561 
01562   mpf_EGlpNumInitVar (alpha);
01563   mpf_EGlpNumInitVar (fi.totinfeas);
01564   mpf_EGlpNumInitVar (pr.dinfeas);
01565   mpf_EGlpNumInitVar (pr.pinfeas);
01566   mpf_EGlpNumInitVar (rs.tz);
01567   mpf_EGlpNumInitVar (rs.lbound);
01568   mpf_EGlpNumInitVar (rs.ecoeff);
01569   mpf_EGlpNumInitVar (rs.pivotval);
01570   mpf_EGlpNumZero (alpha);
01571 
01572   mpf_ILLfct_update_counts (lp, CNT_PPHASE1ITER, 0, mpf_zeroLpNum);
01573   it->nextstep = SIMPLEX_CONTINUE;
01574   it->nextphase = PRIMAL_PHASEI;
01575   lp->final_phase = PRIMAL_PHASEI;
01576   it->nosolve++;
01577 
01578   if (it->newphase != 0)
01579   {
01580     mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->ip_tol);
01581     if (it->newphase == SIMPLEX_PHASE_NEW)
01582     {
01583       it->noprog = 0;
01584       if (it->sdisplay)
01585       {
01586         printf ("starting primal phase I, nosolve %d\n", it->nosolve);
01587         fflush (stdout);
01588       }
01589     }
01590     it->newphase = 0;
01591     it->nosolve = 0;
01592     mpf_EGlpNumCopy (it->prevobj, lp->pinfeas);
01593     lp->pIpiz = mpf_EGlpNumAllocArray (lp->nrows);
01594     lp->pIdz = mpf_EGlpNumAllocArray (lp->nnbasic);
01595 
01596     mpf_ILLfct_compute_phaseI_piz (lp);
01597     if (pinf->p_strategy == COMPLETE_PRICING)
01598     {
01599       mpf_ILLfct_compute_phaseI_dz (lp);
01600 #if USEHEAP > 0
01601       mpf_ILLprice_free_heap (pinf);
01602 #endif
01603       mpf_ILLprice_compute_dual_inf (lp, pinf, NULL, 0, PRIMAL_PHASEI);
01604 #if USEHEAP > 0
01605       rval = mpf_ILLprice_test_for_heap (lp, pinf, lp->nnbasic, pinf->d_scaleinf,
01606                                      PRIMAL_SIMPLEX, 0);
01607       CHECKRVALG (rval, CLEANUP);
01608 #endif
01609     }
01610     else if (pinf->p_strategy == MULTI_PART_PRICING)
01611       mpf_ILLprice_init_mpartial_price (lp, pinf, cphase, COL_PRICING);
01612   }
01613 
01614   mpf_monitor_iter (lp, pinf, it, cphase);
01615   if (it->nextstep == SIMPLEX_TERMINATE || it->nextstep == SIMPLEX_RESUME ||
01616       it->newphase != 0)
01617     ILL_CLEANUP;
01618 
01619   mpf_ILLprice_primal (lp, pinf, &pr, cphase);
01620   ILL_IFTRACE2 ("%s:after_price\n", __func__);
01621 
01622   if (pr.price_stat == PRICE_OPTIMAL)
01623   {
01624     if (it->sdisplay > 1)
01625     {
01626       printf ("primal phase I seemingly done\n");
01627       printf ("retesting soln\n");
01628       fflush (stdout);
01629     }
01630     rval = mpf_ILLsimplex_retest_psolution (lp, pinf, cphase, &fi);
01631 
01632     CHECKRVALG (rval, CLEANUP);
01633     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEI);
01634 
01635     if (fi.pstatus == PRIMAL_FEASIBLE)
01636     {
01637       it->nextphase = PRIMAL_PHASEII;
01638     }
01639     else if (fi.dstatus == DUAL_FEASIBLE)
01640     {
01641       it->solstatus = ILL_LP_SOLVED;
01642       it->nextstep = SIMPLEX_TERMINATE;
01643     }
01644     ILL_CLEANUP;
01645   }
01646 
01647   mpf_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[pr.eindex]);
01648   mpf_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, mpf_zeroLpNum);
01649   mpf_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, mpf_zeroLpNum);
01650 
01651   mpf_ILLratio_pI_test (lp, pr.eindex, pr.dir, &rs);
01652   //ILL_IFTRACE(":%d",rs.lindex);
01653 
01654   if (rs.ratio_stat == RATIO_FAILED)
01655   {
01656     /*
01657      * rval = E_SIMPLEX_ERROR;
01658      * it->solstatus = ILL_PPHASEI_ERROR;
01659      */
01660     //ILL_IFTRACE("ratio_failed\n");
01661     it->algorithm = DUAL_SIMPLEX;
01662     it->nextstep = SIMPLEX_RESUME;
01663     it->resumeid = SIMPLEX_RESUME_NUMER;
01664     /* this is to force to exit in the case of bad basis */
01665     it->n_restart++;
01666     mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01667     mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01668     //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
01669     ILL_CLEANUP;
01670   }
01671   else if (rs.ratio_stat == RATIO_NEGATIVE)
01672   {
01673     mpf_t itol;
01674 
01675     mpf_EGlpNumInitVar (itol);
01676     //ILL_IFTRACE("ratio_negative\n");
01677     mpf_EGlpNumCopy (itol, lp->tol->ip_tol);
01678     mpf_EGlpNumZero (lp->tol->ip_tol);
01679     mpf_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01680     if (!mpf_test_progress (lp->pinfeas, it->prevobj))
01681       it->noprog++;
01682     else
01683     {
01684       mpf_EGlpNumCopy (it->prevobj, lp->pinfeas);
01685       it->noprog = 0;
01686     }
01687     mpf_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01688     mpf_EGlpNumCopy (lp->tol->ip_tol, itol);
01689     mpf_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01690     mpf_ILLfct_update_ppI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.eindex,
01691                               rs.lindex, mpf_zeroLpNum);
01692     mpf_EGlpNumClearVar (itol);
01693   }
01694   else if (rs.ratio_stat == RATIO_NOBCHANGE)
01695   {
01696     //ILL_IFTRACE("ratio_nobchange\n");
01697     mpf_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01698     if (!mpf_test_progress (lp->pinfeas, it->prevobj))
01699       it->noprog++;
01700     else
01701     {
01702       mpf_EGlpNumCopy (it->prevobj, lp->pinfeas);
01703       it->noprog = 0;
01704     }
01705 
01706     //ILL_IFTRACE("%s:a\n",__func__);
01707     mpf_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
01708     mpf_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01709     mpf_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01710     mpf_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
01711 #if DENSE_PI > 0
01712     mpf_fct_test_workvector (lp);
01713     mpf_fct_test_pfeasible (lp);
01714 #endif
01715     mpf_ILLfct_update_ppI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.eindex,
01716                               rs.lindex, mpf_zeroLpNum);
01717   }
01718   else if (rs.ratio_stat == RATIO_BCHANGE)
01719   {
01720     //ILL_IFTRACE("ratio_bchange\n");
01721     mpf_EGlpNumCopyFrac (alpha, lp->pIdz[pr.eindex], rs.pivotval);
01722     mpf_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01723 
01724     if (!mpf_test_progress (lp->pinfeas, it->prevobj))
01725     {
01726       if (lp->vtype[lp->nbaz[pr.eindex]] == VFREE ||
01727           lp->vtype[lp->baz[rs.lindex]] == VARTIFICIAL)
01728       {
01729         if (it->noprog > 0)
01730           it->noprog--;
01731       }
01732       else
01733         it->noprog++;
01734     }
01735     else
01736     {
01737       mpf_EGlpNumCopy (it->prevobj, lp->pinfeas);
01738       it->noprog = 0;
01739     }
01740 
01741     mpf_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
01742     mpf_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, mpf_zeroLpNum);
01743     if (pinf->p_strategy == COMPLETE_PRICING)
01744     {
01745       //ILL_IFTRACE("%s:a\n",__func__);
01746       mpf_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
01747       mpf_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, mpf_zeroLpNum);
01748 
01749       if (pinf->pI_price == QS_PRICE_PSTEEP)
01750       {
01751         mpf_ILLfct_compute_psteep_upv (lp, wz);
01752       }
01753     }
01754 
01755     rval =
01756       mpf_ILLprice_update_pricing_info (lp, pinf, cphase, wz, pr.eindex, rs.lindex,
01757                                     rs.pivotval);
01758     CHECKRVALG (rval, CLEANUP);
01759 
01760     //ILL_IFTRACE("%s:b:%d\n",__func__,rs.lindex);
01761     mpf_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
01762     mpf_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01763     //ILL_IFTRACE("%s:%d:%d\n",__func__,rs.lindex,lp->srhs.nzcnt);
01764     mpf_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01765     mpf_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
01766 #if DENSE_PI > 0
01767     mpf_fct_test_workvector (lp);
01768     mpf_fct_test_pfeasible (lp);
01769 #endif
01770     rval = mpf_ILLbasis_update (lp, updz, rs.lindex, &refactor, &singular);
01771     CHECKRVALG (rval, CLEANUP);
01772 
01773     if (singular)
01774     {
01775       it->nextstep = SIMPLEX_RESUME;
01776       it->resumeid = SIMPLEX_RESUME_SING;
01777       /* this is to force to exit in the case of bad basis */
01778       it->n_restart++;
01779       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01780       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01781       //fprintf(stderr,"Resume Singular %s:%s:%d\n",__func__,__FILE__,__LINE__);
01782       ILL_CLEANUP;
01783     }
01784     if (!refactor)
01785     {
01786       mpf_ILLfct_update_ppI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.eindex,
01787                                 rs.lindex, alpha);
01788     }
01789     if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
01790     {
01791       mpf_ILLfct_compute_xbz (lp);
01792       mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->ip_tol);
01793       mpf_ILLfct_set_status_values (lp, fi.pstatus, -1, PHASEII, -1);
01794       if (fi.pstatus == PRIMAL_FEASIBLE)
01795         it->nextphase = PRIMAL_PHASEII;
01796 
01797       it->newphase = SIMPLEX_PHASE_RECOMP;
01798       ILL_CLEANUP;
01799     }
01800   }
01801 
01802 #if DENSE_PI > 1
01803   mpf_fct_test_workvector (lp);
01804   fct_test_pi_dz (lp, pinf);
01805 #endif
01806 
01807 CLEANUP:
01808   if (it->nextphase != PRIMAL_PHASEI || it->nextstep == SIMPLEX_RESUME ||
01809       it->newphase != 0 || rval != 0)
01810   {
01811     mpf_EGlpNumFreeArray (lp->pIpiz);
01812     mpf_EGlpNumFreeArray (lp->pIdz);
01813   }
01814   mpf_EGlpNumClearVar (alpha);
01815   mpf_EGlpNumClearVar (fi.totinfeas);
01816   mpf_EGlpNumClearVar (pr.dinfeas);
01817   mpf_EGlpNumClearVar (pr.pinfeas);
01818   mpf_EGlpNumClearVar (rs.tz);
01819   mpf_EGlpNumClearVar (rs.lbound);
01820   mpf_EGlpNumClearVar (rs.ecoeff);
01821   mpf_EGlpNumClearVar (rs.pivotval);
01822   //EG_RETURN(rval);
01823   return rval;
01824 }
01825 
01826 static int mpf_primal_phaseII_step (
01827   mpf_lpinfo * lp,
01828   mpf_price_info * pinf,
01829   mpf_svector * updz,
01830   mpf_svector * wz,
01831   mpf_iter_info * it)
01832 {
01833   int boundch;
01834   int rval = 0;
01835   int bndtype = 0;
01836   int singular = 0;
01837   int refactor = 0;
01838   int ratio_iter = 0;
01839   int cphase = PRIMAL_PHASEII;
01840   mpf_t lbound;
01841   mpf_t alpha;
01842   mpf_feas_info fi;
01843   mpf_ratio_res rs;
01844   mpf_price_res pr;
01845 
01846   mpf_EGlpNumInitVar (alpha);
01847   mpf_EGlpNumInitVar (lbound);
01848   mpf_EGlpNumInitVar (fi.totinfeas);
01849   mpf_EGlpNumInitVar (pr.dinfeas);
01850   mpf_EGlpNumInitVar (pr.pinfeas);
01851   mpf_EGlpNumInitVar (rs.tz);
01852   mpf_EGlpNumInitVar (rs.lbound);
01853   mpf_EGlpNumInitVar (rs.ecoeff);
01854   mpf_EGlpNumInitVar (rs.pivotval);
01855 
01856   mpf_ILLfct_update_counts (lp, CNT_PPHASE2ITER, 0, mpf_zeroLpNum);
01857   it->nextstep = SIMPLEX_CONTINUE;
01858   it->nextphase = PRIMAL_PHASEII;
01859   lp->final_phase = PRIMAL_PHASEII;
01860   it->nosolve++;
01861 
01862   if (it->newphase != 0)
01863   {
01864     mpf_ILLfct_compute_pobj (lp);
01865     if (it->newphase == SIMPLEX_PHASE_NEW)
01866     {
01867       it->noprog = 0;
01868       if (it->sdisplay)
01869       {
01870         printf ("starting primal phase II, nosolve %d\n", it->nosolve);
01871         fflush (stdout);
01872       }
01873     }
01874     it->newphase = 0;
01875     it->nosolve = 0;
01876     mpf_EGlpNumCopy (it->prevobj, lp->pobjval);
01877     mpf_ILLfct_compute_piz (lp);
01878     if (pinf->p_strategy == COMPLETE_PRICING)
01879     {
01880       mpf_ILLfct_compute_dz (lp);
01881 #if USEHEAP > 0
01882       mpf_ILLprice_free_heap (pinf);
01883 #endif
01884       mpf_ILLprice_compute_dual_inf (lp, pinf, NULL, 0, PRIMAL_PHASEII);
01885 #if USEHEAP > 0
01886       rval = mpf_ILLprice_test_for_heap (lp, pinf, lp->nnbasic, pinf->d_scaleinf,
01887                                      PRIMAL_SIMPLEX, 0);
01888       CHECKRVALG (rval, CLEANUP);
01889 #endif
01890     }
01891     else if (pinf->p_strategy == MULTI_PART_PRICING)
01892     {
01893       mpf_ILLprice_init_mpartial_price (lp, pinf, cphase, COL_PRICING);
01894     }
01895   }
01896 
01897   mpf_monitor_iter (lp, pinf, it, cphase);
01898   if (it->nextstep == SIMPLEX_TERMINATE || it->nextstep == SIMPLEX_RESUME ||
01899       it->newphase != 0)
01900     ILL_CLEANUP;
01901 
01902   mpf_ILLprice_primal (lp, pinf, &pr, cphase);
01903 
01904   if (pr.price_stat == PRICE_OPTIMAL)
01905   {
01906     //ILL_IFTRACE("%s:PRICE_OPTIMAL\n",__func__);
01907     if (lp->nbchange != 0)
01908     {
01909       if (it->sdisplay > 1)
01910       {
01911         printf ("unrolling %d bound shifts\n", lp->nbchange);
01912         fflush (stdout);
01913       }
01914       mpf_ILLfct_unroll_bound_change (lp);
01915       mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
01916       mpf_ILLfct_set_status_values (lp, fi.pstatus, -1, PHASEII, -1);
01917 
01918        /*HHH*/ mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01919       /*HHH* printf ("primal (opt) infeas %.6f\n", lp->pinfeas); fflush (stdout); 
01920        *HHH* printf ("dual (opt) infeas %.6f\n", lp->dinfeas); fflush (stdout);*/
01921 
01922       if (fi.pstatus != PRIMAL_FEASIBLE)
01923       {
01924         it->algorithm = DUAL_SIMPLEX;
01925         it->nextstep = SIMPLEX_RESUME;
01926         it->resumeid = SIMPLEX_RESUME_UNSHIFT;
01927         it->pricetype = QS_PRICE_DDEVEX;
01928         /* this is to force to exit in the case of bad basis */
01929         //fprintf(stderr,"Resume Unshift %s:%s:%d\n",__func__,__FILE__,__LINE__);
01930         mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01931         mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01932         it->n_restart++;
01933         ILL_CLEANUP;
01934         /*
01935          * it->nextphase = PRIMAL_PHASEI;
01936          * lp->tol->ip_tol /= 5.0;
01937          * lp->tol->id_tol /= 5.0;
01938          * ILL_CLEANUP;
01939          */
01940       }
01941     }
01942 
01943     if (it->sdisplay > 1)
01944     {
01945       printf ("problem seemingly solved\n");
01946       printf ("seemingly opt = %f\nretesting soln\n",
01947               mpf_EGlpNumToLf (lp->pobjval));
01948       fflush (stdout);
01949     }
01950     rval = mpf_ILLsimplex_retest_psolution (lp, pinf, cphase, &fi);
01951     CHECKRVALG (rval, CLEANUP);
01952     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
01953 
01954     if (fi.pstatus == PRIMAL_INFEASIBLE)
01955     {
01956       it->nextphase = PRIMAL_PHASEI;
01957       mpf_EGlpNumDivUiTo (lp->tol->ip_tol, 5);
01958       mpf_EGlpNumDivUiTo (lp->tol->id_tol, 5);
01959       ILL_IFTRACE ("%s:PINF:%lg\n", __func__, mpf_EGlpNumToLf (lp->tol->ip_tol));
01960     }
01961     else if (fi.dstatus == DUAL_FEASIBLE)
01962     {
01963       //ILL_IFTRACE("%s:PFEAS_DFEAS\n",__func__);
01964       it->solstatus = ILL_LP_SOLVED;
01965       mpf_EGlpNumCopy (lp->objval, lp->pobjval);
01966       it->nextstep = SIMPLEX_TERMINATE;
01967     }
01968     else
01969       ILL_IFTRACE ("%s:DINF:%la:%lf\n", __func__, mpf_EGlpNumToLf (lp->dinfeas),
01970                    mpf_EGlpNumToLf (lp->dinfeas));
01971     ILL_CLEANUP;
01972   }
01973 
01974   mpf_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[pr.eindex]);
01975   mpf_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, mpf_zeroLpNum);
01976   mpf_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, mpf_zeroLpNum);
01977   ratio_iter = 0;
01978   do
01979   {
01980     mpf_ILLratio_pII_test (lp, pr.eindex, pr.dir, &rs);
01981     //ILL_IFTRACE("all:%d",rs.lindex);
01982     mpf_EGlpNumCopy (lbound, rs.lbound);
01983     boundch = rs.boundch;
01984     ratio_iter++;
01985 
01986     if (boundch)
01987     {
01988       /*
01989        * if (ratio_iter > PARAM_PRATIOTESTS){
01990        * lbound = lp->xbz[rs.lindex];
01991        * boundch = 0;
01992        * }
01993        */
01994       boundch = 0;
01995       bndtype = (rs.lvstat == STAT_UPPER) ? BOUND_UPPER : BOUND_LOWER;
01996       rval = mpf_ILLfct_bound_shift (lp, lp->baz[rs.lindex], bndtype, lbound);
01997       CHECKRVALG (rval, CLEANUP);
01998     }
01999   } while (boundch);
02000 
02001   if (rs.ratio_stat == RATIO_FAILED)
02002   {
02003     //ILL_IFTRACE(":%d",rs.lindex);
02004     /*
02005      * rval = E_SIMPLEX_ERROR;
02006      * it->solstatus = ILL_PPHASEII_ERROR;
02007      */
02008     it->algorithm = DUAL_SIMPLEX;
02009     it->nextstep = SIMPLEX_RESUME;
02010     it->resumeid = SIMPLEX_RESUME_NUMER;
02011     /* this is to force to exit in the case of bad basis */
02012     it->n_restart++;
02013     mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02014     mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02015     //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
02016     ILL_CLEANUP;
02017   }
02018   else if (rs.ratio_stat == RATIO_UNBOUNDED)
02019   {
02020     //ILL_IFTRACE(":%d",rs.lindex);
02021     if (lp->nbchange != 0)
02022     {
02023       if (it->sdisplay > 1)
02024       {
02025         printf ("unrolling %d bound shifts\n", lp->nbchange);
02026         fflush (stdout);
02027       }
02028       mpf_ILLfct_unroll_bound_change (lp);
02029     }
02030     mpf_ILLfct_set_status_values (lp, PRIMAL_UNBOUNDED, -1, PHASEII, -1);
02031     it->solstatus = ILL_LP_SOLVED;
02032     it->nextstep = SIMPLEX_TERMINATE;
02033     ILL_CLEANUP;
02034   }
02035   else if (rs.ratio_stat == RATIO_NOBCHANGE)
02036   {
02037     //ILL_IFTRACE(":%d",rs.lindex);
02038     mpf_EGlpNumAddInnProdTo (lp->pobjval, rs.tz, lp->dz[pr.eindex]);
02039     mpf_EGlpNumCopy (lp->objval, lp->pobjval);
02040     if (!mpf_test_progress (lp->pobjval, it->prevobj))
02041       it->noprog++;
02042     else
02043     {
02044       mpf_EGlpNumCopy (it->prevobj, lp->pobjval);
02045       it->noprog = 0;
02046     }
02047 
02048     //ILL_IFTRACE("%s:c:%d\n",__func__,rs.lindex);
02049     mpf_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
02050     mpf_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
02051     if (pinf->p_strategy == COMPLETE_PRICING)
02052       mpf_ILLprice_compute_dual_inf (lp, pinf, &pr.eindex, 1, PRIMAL_PHASEII);
02053     else if (pinf->p_strategy == MULTI_PART_PRICING)
02054       mpf_ILLprice_update_mpartial_price (lp, pinf, cphase, COL_PRICING);
02055   }
02056   else if (rs.ratio_stat == RATIO_BCHANGE)
02057   {
02058     mpf_EGlpNumCopyFrac (alpha, lp->dz[pr.eindex], rs.pivotval);
02059     mpf_EGlpNumAddInnProdTo (lp->pobjval, rs.tz, lp->dz[pr.eindex]);
02060     mpf_EGlpNumCopy (lp->objval, lp->pobjval);
02061 
02062     if (!mpf_test_progress (lp->pobjval, it->prevobj))
02063     {
02064       //ILL_IFTRACE(":%d",rs.lindex);
02065       if (lp->vtype[lp->nbaz[pr.eindex]] == VFREE ||
02066           lp->vtype[lp->baz[rs.lindex]] == VARTIFICIAL)
02067       {
02068         if (it->noprog > 0)
02069           it->noprog--;
02070       }
02071       else
02072         it->noprog++;
02073     }
02074     else
02075     {
02076       mpf_EGlpNumCopy (it->prevobj, lp->pobjval);
02077       it->noprog = 0;
02078     }
02079 
02080     //ILL_IFTRACE(":%d",rs.lindex);
02081     mpf_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
02082     mpf_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, mpf_zeroLpNum);
02083     if (pinf->p_strategy == COMPLETE_PRICING)
02084     {
02085       //ILL_IFTRACE("%s:b\n",__func__);
02086       mpf_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02087       mpf_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, mpf_zeroLpNum);
02088       if (pinf->pII_price == QS_PRICE_PSTEEP)
02089         mpf_ILLfct_compute_psteep_upv (lp, wz);
02090     }
02091     rval =
02092       mpf_ILLprice_update_pricing_info (lp, pinf, cphase, wz, pr.eindex, rs.lindex,
02093                                     rs.pivotval);
02094     CHECKRVALG (rval, CLEANUP);
02095 
02096     //ILL_IFTRACE("%s:d:%d\n",__func__,rs.lindex);
02097     mpf_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
02098     mpf_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
02099     rval = mpf_ILLbasis_update (lp, updz, rs.lindex, &refactor, &singular);
02100     CHECKRVALG (rval, CLEANUP);
02101 
02102     if (singular)
02103     {
02104       it->nextstep = SIMPLEX_RESUME;
02105       it->resumeid = SIMPLEX_RESUME_SING;
02106       /* this is to force to exit in the case of bad basis */
02107       it->n_restart++;
02108       //fprintf(stderr,"Resume Singular %s:%s:%d\n",__func__,__FILE__,__LINE__);
02109       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02110       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02111       ILL_CLEANUP;
02112     }
02113     if (!refactor)
02114     {
02115       mpf_ILLfct_update_piz (lp, alpha);
02116 
02117       if (pinf->p_strategy == COMPLETE_PRICING)
02118       {
02119         mpf_ILLfct_update_dz (lp, pr.eindex, alpha);
02120         mpf_ILLprice_compute_dual_inf (lp, pinf, lp->zA.indx, lp->zA.nzcnt,
02121                                    PRIMAL_PHASEII);
02122         mpf_ILLfct_update_counts (lp, CNT_ZARAVG, lp->zA.nzcnt, mpf_zeroLpNum);
02123       }
02124       else if (pinf->p_strategy == MULTI_PART_PRICING)
02125       {
02126         mpf_ILLprice_update_mpartial_price (lp, pinf, cphase, COL_PRICING);
02127       }
02128     }
02129     if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02130     {
02131       mpf_ILLfct_compute_xbz (lp);
02132       it->newphase = SIMPLEX_PHASE_RECOMP;
02133     }
02134   }
02135 
02136 CLEANUP:
02137   mpf_EGlpNumClearVar (alpha);
02138   mpf_EGlpNumClearVar (lbound);
02139   mpf_EGlpNumClearVar (fi.totinfeas);
02140   mpf_EGlpNumClearVar (pr.dinfeas);
02141   mpf_EGlpNumClearVar (pr.pinfeas);
02142   mpf_EGlpNumClearVar (rs.tz);
02143   mpf_EGlpNumClearVar (rs.lbound);
02144   mpf_EGlpNumClearVar (rs.ecoeff);
02145   mpf_EGlpNumClearVar (rs.pivotval);
02146   //EG_RETURN(rval);
02147   return rval;
02148 }
02149 
02150 static int mpf_dual_phaseI_step (
02151   mpf_lpinfo * lp,
02152   mpf_price_info * pinf,
02153   mpf_svector * updz,
02154   mpf_svector * wz,
02155   mpf_iter_info * it)
02156 {
02157   int rval = 0;
02158   int singular = 0;
02159   int refactor = 0;
02160   int cphase = DUAL_PHASEI;
02161   mpf_t alpha;
02162   mpf_t alpha1;
02163   mpf_feas_info fi;
02164   mpf_ratio_res rs;
02165   mpf_price_res pr;
02166 
02167   mpf_EGlpNumInitVar (alpha);
02168   mpf_EGlpNumInitVar (alpha1);
02169   mpf_EGlpNumInitVar (fi.totinfeas);
02170   mpf_EGlpNumInitVar (pr.dinfeas);
02171   mpf_EGlpNumInitVar (pr.pinfeas);
02172   mpf_EGlpNumInitVar (rs.tz);
02173   mpf_EGlpNumInitVar (rs.lbound);
02174   mpf_EGlpNumInitVar (rs.ecoeff);
02175   mpf_EGlpNumInitVar (rs.pivotval);
02176   mpf_EGlpNumZero (alpha1);
02177 
02178   mpf_ILLfct_update_counts (lp, CNT_DPHASE1ITER, 0, mpf_zeroLpNum);
02179   it->nextstep = SIMPLEX_CONTINUE;
02180   it->nextphase = DUAL_PHASEI;
02181   lp->final_phase = DUAL_PHASEI;
02182   it->nosolve++;
02183 
02184   if (it->newphase != 0)
02185   {
02186     mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->id_tol);
02187     if (it->newphase == SIMPLEX_PHASE_NEW)
02188     {
02189       it->noprog = 0;
02190       if (it->sdisplay)
02191       {
02192         printf ("starting dual phase I, nosolve %d\n", it->nosolve);
02193         fflush (stdout);
02194       }
02195     }
02196     it->newphase = 0;
02197     it->nosolve = 0;
02198     mpf_EGlpNumCopy (it->prevobj, lp->dinfeas);
02199 
02200     mpf_ILLfct_compute_phaseI_xbz (lp);
02201     if (pinf->d_strategy == COMPLETE_PRICING)
02202     {
02203 #if USEHEAP > 0
02204       mpf_ILLprice_free_heap (pinf);
02205 #endif
02206       mpf_ILLprice_compute_primal_inf (lp, pinf, NULL, 0, DUAL_PHASEI);
02207 #if USEHEAP > 0
02208       rval = mpf_ILLprice_test_for_heap (lp, pinf, lp->nrows, pinf->p_scaleinf,
02209                                      DUAL_SIMPLEX, 0);
02210       CHECKRVALG (rval, CLEANUP);
02211 #endif
02212     }
02213     else if (pinf->d_strategy == MULTI_PART_PRICING)
02214     {
02215       mpf_ILLprice_init_mpartial_price (lp, pinf, cphase, ROW_PRICING);
02216     }
02217   }
02218 
02219   mpf_monitor_iter (lp, pinf, it, cphase);
02220   if (it->nextstep == SIMPLEX_TERMINATE || it->nextstep == SIMPLEX_RESUME ||
02221       it->newphase != 0)
02222     ILL_CLEANUP;
02223 
02224   mpf_ILLprice_dual (lp, pinf, cphase, &pr);
02225 
02226   if (pr.price_stat == PRICE_OPTIMAL)
02227   {
02228     if (it->sdisplay > 1)
02229     {
02230       printf ("dual phase I seemingly done\n");
02231       printf ("retesting soln\n");
02232       fflush (stdout);
02233     }
02234 
02235     rval = mpf_ILLsimplex_retest_dsolution (lp, pinf, cphase, &fi);
02236     CHECKRVALG (rval, CLEANUP);
02237     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEI, PHASEII);
02238 
02239     if (fi.dstatus == DUAL_FEASIBLE)
02240     {
02241       it->nextphase = DUAL_PHASEII;
02242     }
02243     else if (fi.pstatus == PRIMAL_FEASIBLE)
02244     {
02245       it->solstatus = ILL_LP_SOLVED;
02246       it->nextstep = SIMPLEX_TERMINATE;
02247     }
02248     it->newphase = SIMPLEX_PHASE_NEW;
02249     ILL_CLEANUP;
02250   }
02251 
02252   mpf_ILLfct_compute_zz (lp, &(lp->zz), pr.lindex);
02253   //ILL_IFTRACE("%s:c\n",__func__);
02254   mpf_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02255   mpf_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, mpf_zeroLpNum);
02256   mpf_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, mpf_zeroLpNum);
02257 
02258   mpf_ILLratio_dI_test (lp, pr.lindex, pr.lvstat, &rs);
02259 
02260   if (rs.ratio_stat == RATIO_FAILED)
02261   {
02262     /*
02263      * rval = E_SIMPLEX_ERROR;
02264      * it->solstatus = ILL_DPHASEI_ERROR;
02265      */
02266     it->algorithm = PRIMAL_SIMPLEX;
02267     it->nextstep = SIMPLEX_RESUME;
02268     it->resumeid = SIMPLEX_RESUME_NUMER;
02269     /* this is to force to exit in the case of bad basis */
02270     it->n_restart++;
02271     mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02272     mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02273     //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
02274     ILL_CLEANUP;
02275   }
02276   else if (rs.ratio_stat == RATIO_BCHANGE)
02277   {
02278     mpf_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[rs.eindex]);
02279     rval = mpf_ILLfct_test_pivot (lp, pr.lindex, ROW_PIVOT, rs.pivotval);
02280     if (rval)
02281     {
02282       it->n_pivot_fail++;
02283       if (it->n_pivot_fail > SIMPLEX_MAX_PIVOT_FAIL)
02284       {
02285         it->n_pivot_fail = 0;
02286         /* this is to force to exit in the case of bad basis */
02287         it->n_restart++;
02288         it->algorithm = PRIMAL_SIMPLEX;
02289         it->nextstep = SIMPLEX_RESUME;
02290         it->resumeid = SIMPLEX_RESUME_NUMER;
02291         mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02292         mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02293         //fprintf(stderr,"Resume Pivot %s:%s:%d\n",__func__,__FILE__,__LINE__);
02294         rval = 0;
02295         ILL_CLEANUP;
02296       }
02297       rval = mpf_ILLbasis_factor (lp, &singular);
02298       if (singular)
02299         MESSAGE (__QS_SB_VERB, "Singular basis found!");
02300       CHECKRVALG (rval, CLEANUP);
02301       if (singular == 0)
02302         refactor = 1;
02303       goto END;
02304     }
02305     mpf_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, mpf_zeroLpNum);
02306     mpf_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, mpf_zeroLpNum);
02307 
02308     if (pinf->dI_price == QS_PRICE_DSTEEP)
02309       mpf_ILLfct_compute_dsteep_upv (lp, wz);
02310     rval =
02311       mpf_ILLprice_update_pricing_info (lp, pinf, cphase, wz, rs.eindex, pr.lindex,
02312                                     rs.pivotval);
02313     CHECKRVALG (rval, CLEANUP);
02314 
02315     mpf_EGlpNumSubTo (lp->dinfeas, lp->upd.c_obj);
02316 
02317     if (!mpf_test_progress (lp->dinfeas, it->prevobj))
02318     {
02319       if (lp->vtype[lp->baz[pr.lindex]] == VARTIFICIAL ||
02320           lp->vtype[lp->nbaz[rs.eindex]] == VFREE)
02321       {
02322         if (it->noprog > 0)
02323           it->noprog--;
02324       }
02325       else
02326         it->noprog++;
02327     }
02328     else
02329     {
02330       mpf_EGlpNumCopy (it->prevobj, lp->dinfeas);
02331       it->noprog = 0;
02332     }
02333 
02334     mpf_EGlpNumCopyFrac (alpha, lp->dz[rs.eindex], rs.pivotval);
02335     mpf_EGlpNumCopyFrac (alpha1, lp->xbz[pr.lindex], rs.pivotval);
02336 
02337     mpf_ILLfct_update_piz (lp, alpha);
02338     mpf_ILLfct_update_dz (lp, rs.eindex, alpha);
02339     mpf_ILLfct_update_dfeas (lp, rs.eindex, &(lp->srhs));
02340     mpf_ILLfct_compute_dpIy (lp, &(lp->srhs), &(lp->ssoln));
02341     mpf_ILLfct_update_basis_info (lp, rs.eindex, pr.lindex, pr.lvstat);
02342 
02343 #if DENSE_PI > 0
02344     mpf_fct_test_workvector (lp);
02345     mpf_fct_test_dfeasible (lp);
02346 #endif
02347     rval = mpf_ILLbasis_update (lp, updz, pr.lindex, &refactor, &singular);
02348     CHECKRVALG (rval, CLEANUP);
02349 
02350 #if DENSE_NORM > 0
02351     mpf_test_dsteep_norms (lp, pinf);
02352 #endif
02353 
02354     mpf_ILLfct_update_dpI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.lindex,
02355                               alpha1);
02356 
02357   END:
02358     if (singular)
02359     {
02360       it->nextstep = SIMPLEX_RESUME;
02361       it->resumeid = SIMPLEX_RESUME_SING;
02362       /* this is to force to exit in the case of bad basis */
02363       it->n_restart++;
02364       //fprintf(stderr,"Resume Singular %s:%s:%d\n",__func__,__FILE__,__LINE__);
02365       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02366       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02367       ILL_CLEANUP;
02368     }
02369     if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02370     {
02371       mpf_ILLfct_compute_piz (lp);
02372       mpf_ILLfct_compute_dz (lp);
02373       mpf_ILLfct_dual_adjust (lp, mpf_zeroLpNum);
02374       mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->id_tol);
02375       mpf_ILLfct_set_status_values (lp, -1, fi.dstatus, -1, PHASEII);
02376       if (fi.dstatus == DUAL_FEASIBLE)
02377         it->nextphase = DUAL_PHASEII;
02378 
02379       it->newphase = SIMPLEX_PHASE_RECOMP;
02380       ILL_CLEANUP;
02381     }
02382   }
02383 
02384 #if DENSE_PI > 1
02385   mpf_fct_test_workvector (lp);
02386   mpf_fct_test_pI_x (lp, pinf);
02387 #endif
02388 
02389 CLEANUP:
02390   mpf_EGlpNumClearVar (alpha);
02391   mpf_EGlpNumClearVar (alpha1);
02392   mpf_EGlpNumClearVar (fi.totinfeas);
02393   mpf_EGlpNumClearVar (pr.dinfeas);
02394   mpf_EGlpNumClearVar (pr.pinfeas);
02395   mpf_EGlpNumClearVar (rs.tz);
02396   mpf_EGlpNumClearVar (rs.lbound);
02397   mpf_EGlpNumClearVar (rs.ecoeff);
02398   mpf_EGlpNumClearVar (rs.pivotval);
02399   //EG_RETURN(rval);
02400   return rval;
02401 }
02402 
02403 static int mpf_dual_phaseII_step (
02404   mpf_lpinfo * lp,
02405   mpf_price_info * pinf,
02406   mpf_svector * updz,
02407   mpf_svector * wz,
02408   mpf_iter_info * it)
02409 {
02410   int coeffch;
02411   int rval = 0;
02412   int singular = 0;
02413   int refactor = 0;
02414   int ratio_iter = 0;
02415   int cphase = DUAL_PHASEII;
02416   int lcol, ecol;
02417   int estat, newphase;
02418   mpf_t x_bi, v_l, eval;
02419   mpf_t ecoeff;
02420   mpf_t alpha;
02421   mpf_t alpha1;
02422   mpf_feas_info fi;
02423   mpf_ratio_res rs;
02424   mpf_price_res pr;
02425 
02426   mpf_EGlpNumInitVar (x_bi);
02427   mpf_EGlpNumInitVar (v_l);
02428   mpf_EGlpNumInitVar (eval);
02429   mpf_EGlpNumInitVar (ecoeff);
02430   mpf_EGlpNumInitVar (alpha);
02431   mpf_EGlpNumInitVar (alpha1);
02432   mpf_EGlpNumInitVar (fi.totinfeas);
02433   mpf_EGlpNumInitVar (pr.dinfeas);
02434   mpf_EGlpNumInitVar (pr.pinfeas);
02435   mpf_EGlpNumInitVar (rs.tz);
02436   mpf_EGlpNumInitVar (rs.lbound);
02437   mpf_EGlpNumInitVar (rs.ecoeff);
02438   mpf_EGlpNumInitVar (rs.pivotval);
02439   mpf_EGlpNumZero (rs.ecoeff);
02440   mpf_EGlpNumZero (alpha1);
02441 
02442   mpf_ILLfct_update_counts (lp, CNT_DPHASE2ITER, 0, mpf_zeroLpNum);
02443   it->nextstep = SIMPLEX_CONTINUE;
02444   it->nextphase = DUAL_PHASEII;
02445   lp->final_phase = DUAL_PHASEII;
02446   newphase = it->newphase;
02447   it->nosolve++;
02448 
02449   if (it->newphase != 0)
02450   {
02451     mpf_ILLfct_compute_dobj (lp);
02452     if (it->newphase == SIMPLEX_PHASE_NEW)
02453     {
02454       it->noprog = 0;
02455       if (it->sdisplay)
02456       {
02457         printf ("starting dual phase II, nosolve %d\n", it->nosolve);
02458         fflush (stdout);
02459       }
02460     }
02461     it->newphase = 0;
02462     it->nosolve = 0;
02463     mpf_EGlpNumCopy (it->prevobj, lp->dobjval);
02464     mpf_ILLfct_compute_xbz (lp);
02465 
02466     if (pinf->d_strategy == COMPLETE_PRICING)
02467     {
02468 #if USEHEAP > 0
02469       mpf_ILLprice_free_heap (pinf);
02470 #endif
02471       mpf_ILLprice_compute_primal_inf (lp, pinf, NULL, 0, DUAL_PHASEII);
02472 #if USEHEAP > 0
02473       rval = mpf_ILLprice_test_for_heap (lp, pinf, lp->nrows, pinf->p_scaleinf,
02474                                      DUAL_SIMPLEX, 0);
02475       CHECKRVALG (rval, CLEANUP);
02476 #endif
02477     }
02478     else if (pinf->d_strategy == MULTI_PART_PRICING)
02479     {
02480       mpf_ILLprice_init_mpartial_price (lp, pinf, cphase, ROW_PRICING);
02481     }
02482   }
02483 
02484   mpf_monitor_iter (lp, pinf, it, cphase);
02485   if (it->nextstep == SIMPLEX_TERMINATE || it->nextstep == SIMPLEX_RESUME ||
02486       it->newphase != 0)
02487     ILL_CLEANUP;
02488 
02489   mpf_ILLprice_dual (lp, pinf, cphase, &pr);
02490 
02491   if (pr.price_stat == PRICE_OPTIMAL)
02492   {
02493     if (lp->ncchange != 0)
02494     {
02495       if (it->sdisplay > 1)
02496       {
02497         printf ("unrolling %d coef shifts\n", lp->ncchange);
02498         fflush (stdout);
02499       }
02500       mpf_ILLfct_unroll_coef_change (lp);
02501       mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
02502       mpf_ILLfct_set_status_values (lp, -1, fi.dstatus, -1, PHASEII);
02503 
02504        /*HHH*/ mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
02505       /*HHH* printf ("dual (opt) infeas %.6f\n", lp->dinfeas); fflush (stdout);
02506        *HHH* printf ("primal (opt) infeas %.6f\n", lp->pinfeas); fflush (stdout);*/
02507 
02508       if (fi.dstatus != DUAL_FEASIBLE)
02509       {
02510         it->algorithm = PRIMAL_SIMPLEX;
02511         it->nextstep = SIMPLEX_RESUME;
02512         it->resumeid = SIMPLEX_RESUME_UNSHIFT;
02513         it->pricetype = QS_PRICE_PDEVEX;
02514         /* this is to force to exit in the case of bad basis */
02515         it->n_restart++;
02516         mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02517         mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02518         //fprintf(stderr,"Resume Unshift %s:%s:%d\n",__func__,__FILE__,__LINE__);
02519         ILL_CLEANUP;
02520         /*
02521          * it->nextphase = DUAL_PHASEI;
02522          * lp->tol->ip_tol /= 5.0;
02523          * lp->tol->id_tol /= 5.0;
02524          * ILL_CLEANUP;
02525          */
02526       }
02527     }
02528     if (it->sdisplay > 1)
02529     {
02530       //ILL_IFTRACE("\t%s:%s:%d\n",__func__,__FILE__,__LINE__);
02531       printf ("problem seemingly solved\n");
02532       printf ("seemingly dual opt = %f\n", mpf_EGlpNumToLf (lp->dobjval));
02533       printf ("retesting soln\n");
02534       fflush (stdout);
02535     }
02536 
02537     rval = mpf_ILLsimplex_retest_dsolution (lp, pinf, cphase, &fi);
02538     CHECKRVALG (rval, CLEANUP);
02539     mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
02540 
02541     if (fi.dstatus == DUAL_INFEASIBLE)
02542     {
02543       ILL_IFTRACE ("DUAL_INFEAS: %s\n", __func__);
02544       it->nextphase = DUAL_PHASEI;
02545       mpf_EGlpNumDivUiTo (lp->tol->ip_tol, 5);
02546       mpf_EGlpNumDivUiTo (lp->tol->id_tol, 5);
02547     }
02548     else if (fi.pstatus == PRIMAL_FEASIBLE)
02549     {
02550       ILL_IFTRACE ("PRIM_FEAS: %s\n", __func__);
02551       mpf_EGlpNumCopy (lp->objval, lp->dobjval);
02552       it->solstatus = ILL_LP_SOLVED;
02553       it->nextstep = SIMPLEX_TERMINATE;
02554     }
02555     else
02556       ILL_IFTRACE ("PRIM_INFEAS: %s\n", __func__);
02557     ILL_CLEANUP;
02558   }
02559 
02560   mpf_ILLfct_compute_zz (lp, &(lp->zz), pr.lindex);
02561   //ILL_IFTRACE("%s:d\n",__func__);
02562   mpf_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02563   mpf_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, mpf_zeroLpNum);
02564   mpf_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, mpf_zeroLpNum);
02565 
02566   ratio_iter = 0;
02567   do
02568   {
02569     mpf_ILLratio_longdII_test (lp, pr.lindex, pr.lvstat, &rs);
02570     if (rs.ratio_stat == RATIO_NEGATIVE)
02571     {
02572       if (it->sdisplay > 1)
02573       {
02574         printf ("adjust coefs to remove negative ratio tests\n");
02575         fflush (stdout);
02576       }
02577       mpf_ILLfct_adjust_viol_coefs (lp);
02578       mpf_ILLratio_longdII_test (lp, pr.lindex, pr.lvstat, &rs);
02579       if (rs.ratio_stat == RATIO_NEGATIVE)
02580       {
02581         MESSAGE (__QS_SB_VERB, "internal error: bad ratio test");
02582         fflush (stdout);
02583         rs.ratio_stat = RATIO_FAILED;
02584         break;
02585       }
02586     }
02587 
02588     coeffch = rs.coeffch;
02589     mpf_EGlpNumCopy (ecoeff, rs.ecoeff);
02590     ratio_iter++;
02591 
02592     if (coeffch)
02593     {
02594       /*
02595        * if (ratio_iter > PARAM_DRATIOTESTS){
02596        * ecoeff = lp->cz[lp->nbaz[rs.eindex]] - lp->dz[rs.eindex];
02597        * coeffch = 0;
02598        * }
02599        */
02600       coeffch = 0;
02601       rval = mpf_ILLfct_coef_shift (lp, lp->nbaz[rs.eindex], ecoeff);
02602       CHECKRVALG (rval, CLEANUP);
02603     }
02604     if (rs.ratio_stat == RATIO_BCHANGE)
02605       if (lp->vstat[lp->nbaz[rs.eindex]] == STAT_ZERO)
02606         break;
02607 
02608   } while (coeffch);
02609 
02610   if (rs.ratio_stat == RATIO_FAILED)
02611   {
02612     /*
02613      * rval = E_SIMPLEX_ERROR;
02614      * it->solstatus = ILL_DPHASEII_ERROR;
02615      */
02616     it->algorithm = PRIMAL_SIMPLEX;
02617     it->nextstep = SIMPLEX_RESUME;
02618     it->resumeid = SIMPLEX_RESUME_NUMER;
02619     /* this is to force to exit in the case of bad basis */
02620     it->n_restart++;
02621     mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02622     mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02623     //fprintf(stderr,"Resume Numerical %s:%s:%d\n",__func__,__FILE__,__LINE__);
02624     ILL_CLEANUP;
02625   }
02626   else if (rs.ratio_stat == RATIO_UNBOUNDED)
02627   {
02628     lp->infub_ix = pr.lindex;
02629     if (lp->ncchange != 0)
02630     {
02631       if (it->sdisplay > 1)
02632       {
02633         printf ("unrolling %d coef shifts\n", lp->ncchange);
02634         fflush (stdout);
02635       }
02636       mpf_ILLfct_unroll_coef_change (lp);
02637     }
02638     mpf_ILLfct_set_status_values (lp, -1, DUAL_UNBOUNDED, -1, PHASEII);
02639     it->solstatus = ILL_LP_SOLVED;
02640     it->nextstep = SIMPLEX_TERMINATE;
02641   }
02642   else if (rs.ratio_stat == RATIO_BCHANGE)
02643   {
02644     lcol = lp->baz[pr.lindex];
02645     ecol = lp->nbaz[rs.eindex];
02646 
02647     mpf_ILLfct_compute_yz (lp, &(lp->yjz), updz, ecol);
02648     mpf_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, mpf_zeroLpNum);
02649     mpf_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, mpf_zeroLpNum);
02650     rval = mpf_ILLfct_test_pivot (lp, pr.lindex, ROW_PIVOT, rs.pivotval);
02651     if (rval != 0)
02652     {
02653       it->n_pivot_fail++;
02654       if (it->n_pivot_fail > SIMPLEX_MAX_PIVOT_FAIL)
02655       {
02656         it->n_pivot_fail = 0;
02657         /* this is to force to exit in the case of bad basis */
02658         it->algorithm = PRIMAL_SIMPLEX;
02659         it->nextstep = SIMPLEX_RESUME;
02660         it->resumeid = SIMPLEX_RESUME_NUMER;
02661         it->n_restart++;
02662         mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02663         mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02664         //fprintf(stderr,"Resume Pivot %s:%s:%d\n",__func__,__FILE__,__LINE__);
02665         rval = 0;
02666         ILL_CLEANUP;
02667       }
02668       if (newphase == 0)
02669       {
02670         rval = mpf_ILLbasis_factor (lp, &singular);
02671         CHECKRVALG (rval, CLEANUP);
02672         if (singular)
02673           MESSAGE (__QS_SB_VERB, "Singular basis found!");
02674 #ifdef dbl_QSOPT_CURRENT_PRECICION
02675         if (singular)
02676         {
02677           MESSAGE (__QS_SB_VERB, "Forcing fail!");
02678           rval = QS_LP_CHANGE_PREC;
02679         }
02680 #endif
02681         if (singular == 0)
02682           refactor = 1;
02683         goto END;
02684       }
02685       else
02686       {
02687         if (it->sdisplay > 1)
02688         {
02689           printf ("warning: bad step\n");
02690           fflush (stdout);
02691         }
02692       }
02693     }
02694 
02695     mpf_EGlpNumAddTo (lp->dobjval, lp->upd.c_obj);
02696     mpf_EGlpNumCopy (lp->objval, lp->dobjval);
02697 
02698     if (!mpf_test_progress (lp->dobjval, it->prevobj))
02699     {
02700       if (lp->vtype[lcol] == VARTIFICIAL || lp->vtype[ecol] == VFREE)
02701       {
02702         if (it->noprog > 0)
02703           it->noprog--;
02704       }
02705       else
02706         it->noprog++;
02707     }
02708     else
02709     {
02710       mpf_EGlpNumCopy (it->prevobj, lp->dobjval);
02711       it->noprog = 0;
02712     }
02713 
02714     if (pinf->dII_price == QS_PRICE_DSTEEP)
02715       mpf_ILLfct_compute_dsteep_upv (lp, wz);
02716     rval =
02717       mpf_ILLprice_update_pricing_info (lp, pinf, cphase, wz, rs.eindex, pr.lindex,
02718                                     rs.pivotval);
02719     CHECKRVALG (rval, CLEANUP);
02720 
02721     mpf_EGlpNumCopy (x_bi, lp->xbz[pr.lindex]);
02722     if (pr.lvstat == STAT_LOWER)
02723       mpf_EGlpNumCopy (v_l, lp->lz[lcol]);
02724     else
02725       mpf_EGlpNumCopy (v_l, lp->uz[lcol]);
02726     mpf_EGlpNumCopy (alpha, rs.tz);
02727     if (pr.lvstat == STAT_LOWER)
02728       mpf_EGlpNumSign (alpha);
02729     estat = lp->vstat[ecol];
02730     if (estat == STAT_LOWER)
02731       mpf_EGlpNumCopy (eval, lp->lz[ecol]);
02732     else if (estat == STAT_ZERO)
02733       mpf_EGlpNumZero (eval);
02734     else
02735       mpf_EGlpNumCopy (eval, lp->uz[ecol]);
02736 
02737     mpf_ILLfct_update_piz (lp, alpha);
02738     mpf_ILLfct_update_dz (lp, rs.eindex, alpha);
02739     mpf_ILLfct_update_dIIfeas (lp, rs.eindex, &(lp->srhs));
02740     mpf_ILLfct_compute_dpIIy (lp, &(lp->srhs), &(lp->ssoln));
02741     mpf_EGlpNumCopyDiff (alpha1, x_bi, v_l);
02742     mpf_EGlpNumSubTo (alpha1, lp->upd.dty);
02743     mpf_EGlpNumDivTo (alpha1, rs.pivotval);
02744     mpf_ILLfct_update_basis_info (lp, rs.eindex, pr.lindex, pr.lvstat);
02745     rval = mpf_ILLbasis_update (lp, updz, pr.lindex, &refactor, &singular);
02746     CHECKRVALG (rval, CLEANUP);
02747 
02748     mpf_ILLfct_update_dpII_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), /*rs.eindex,*/
02749                                pr.lindex, eval, alpha1);
02750 
02751 #if DENSE_NORM > 0
02752     mpf_test_dsteep_norms (lp, pinf);
02753 #endif
02754 
02755   END:
02756     if (singular)
02757     {
02758       it->nextstep = SIMPLEX_RESUME;
02759       it->resumeid = SIMPLEX_RESUME_SING;
02760       /* this is to force to exit in the case of bad basis */
02761       it->n_restart++;
02762       mpf_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02763       mpf_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02764       //fprintf(stderr,"Resume Singular %s:%s:%d\n",__func__,__FILE__,__LINE__);
02765       ILL_CLEANUP;
02766     }
02767     if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02768     {
02769       mpf_ILLfct_compute_piz (lp);
02770       mpf_ILLfct_compute_dz (lp);
02771       mpf_ILLfct_dual_adjust (lp, mpf_zeroLpNum);
02772       it->newphase = SIMPLEX_PHASE_RECOMP;
02773     }
02774   }
02775 
02776 #if DENSE_PIIPI > 0
02777   mpf_fct_test_workvector (lp);
02778   if (!refactor)
02779   {
02780     mpf_fct_test_pII_x (lp, pinf);
02781     mpf_fct_test_pII_pi_dz (lp, pinf);
02782   }
02783 #endif
02784 
02785 CLEANUP:
02786   mpf_EGlpNumClearVar (x_bi);
02787   mpf_EGlpNumClearVar (v_l);
02788   mpf_EGlpNumClearVar (eval);
02789   mpf_EGlpNumClearVar (ecoeff);
02790   mpf_EGlpNumClearVar (alpha);
02791   mpf_EGlpNumClearVar (alpha1);
02792   mpf_EGlpNumClearVar (fi.totinfeas);
02793   mpf_EGlpNumClearVar (pr.dinfeas);
02794   mpf_EGlpNumClearVar (pr.pinfeas);
02795   mpf_EGlpNumClearVar (rs.tz);
02796   mpf_EGlpNumClearVar (rs.lbound);
02797   mpf_EGlpNumClearVar (rs.ecoeff);
02798   mpf_EGlpNumClearVar (rs.pivotval);
02799   //EG_RETURN(rval);
02800   return rval;
02801 }
02802 
02803 static void mpf_get_current_stat (
02804   mpf_lp_status_info * p,
02805   int algorithm,
02806   int *bstat)
02807 {
02808   if (p->optimal)
02809     *bstat = OPTIMAL;
02810   else if (algorithm == PRIMAL_SIMPLEX)
02811   {
02812     if (p->primal_feasible)
02813       *bstat = PRIMAL_FEASIBLE;
02814     else if (p->primal_infeasible)
02815       *bstat = PRIMAL_INFEASIBLE;
02816     else if (p->primal_unbounded)
02817       *bstat = PRIMAL_UNBOUNDED;
02818     else
02819       *bstat = NONOPTIMAL;
02820   }
02821   else if (algorithm == DUAL_SIMPLEX)
02822   {
02823     if (p->dual_feasible)
02824       *bstat = DUAL_FEASIBLE;
02825     else if (p->dual_infeasible)
02826       *bstat = DUAL_INFEASIBLE;
02827     else if (p->dual_unbounded)
02828       *bstat = DUAL_UNBOUNDED;
02829     else
02830       *bstat = NONOPTIMAL;
02831   }
02832 }
02833 
02834 int mpf_ILLsimplex_pivotin (
02835   mpf_lpinfo * lp,
02836   mpf_price_info * pinf,
02837   int rcnt,
02838   int *rlist,
02839   int pivot_opt,
02840   int *basis_mod)
02841 {
02842   int i, npiv = 0;
02843   int eindex;
02844   int rval = 0;
02845   int singular = 0;
02846   int refactor = 0;
02847   int *rowmap = lp->O->rowmap;
02848   int *clist = NULL;
02849   mpf_svector wz;
02850   mpf_svector updz;
02851   mpf_t alpha;
02852   mpf_ratio_res rs;
02853   mpf_feas_info fi;
02854 
02855   mpf_EGlpNumInitVar (alpha);
02856   mpf_EGlpNumInitVar (fi.totinfeas);
02857   mpf_EGlpNumInitVar (rs.tz);
02858   mpf_EGlpNumInitVar (rs.lbound);
02859   mpf_EGlpNumInitVar (rs.ecoeff);
02860   mpf_EGlpNumInitVar (rs.pivotval);
02861   mpf_EGlpNumZero (alpha);
02862 
02863   *basis_mod = 0;
02864   if (rcnt <= 0)
02865   {
02866     EG_RETURN (rval);
02867   }
02868 
02869   if (pivot_opt == SIMPLEX_PIVOTINROW)
02870   {
02871     ILL_SAFE_MALLOC (clist, rcnt, int);
02872 
02873     for (i = 0; i < rcnt; i++)
02874       clist[i] = rowmap[rlist[i]];
02875   }
02876   else
02877     clist = rlist;
02878 
02879   for (i = 0; i < rcnt; i++)
02880   {
02881     if (lp->vstat[clist[i]] != STAT_BASIC)
02882     {
02883       *basis_mod = 1;
02884       break;
02885     }
02886   }
02887   if (*basis_mod == 0)
02888   {
02889     if (pivot_opt == SIMPLEX_PIVOTINROW)
02890     {
02891       ILL_IFFREE (clist, int);
02892     }
02893     EG_RETURN (rval);
02894   }
02895 
02896   /* printf ("Forcing vars into basis in mpf_ILLsimplex_pivotin \n"); */
02897   mpf_ILLsvector_init (&wz);
02898   rval = mpf_ILLsvector_alloc (&wz, lp->nrows);
02899   CHECKRVALG (rval, CLEANUP);
02900   mpf_ILLsvector_init (&updz);
02901   rval = mpf_ILLsvector_alloc (&updz, lp->nrows);
02902   CHECKRVALG (rval, CLEANUP);
02903 
02904   mpf_EGlpNumCopy (lp->pobjval, lp->dobjval);
02905   for (i = 0; i < rcnt; i++)
02906   {
02907     if (lp->vstat[clist[i]] == STAT_BASIC)
02908       continue;
02909     npiv++;
02910 
02911     eindex = lp->vindex[clist[i]];
02912     mpf_ILLfct_compute_yz (lp, &(lp->yjz), &updz, lp->nbaz[eindex]);
02913     mpf_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, mpf_zeroLpNum);
02914     mpf_ILLfct_update_counts (lp, CNT_UPNZ, updz.nzcnt, mpf_zeroLpNum);
02915 
02916     mpf_ILLratio_pivotin_test (lp, clist, rcnt, &rs);
02917 
02918     if (rs.ratio_stat == RATIO_UNBOUNDED || rs.ratio_stat == RATIO_FAILED)
02919     {
02920       fprintf (stderr, "Pivot_in failed\n");
02921       rval = E_SIMPLEX_ERROR;
02922       ILL_CLEANUP;
02923     }
02924     else if (rs.ratio_stat == RATIO_BCHANGE)
02925     {
02926       if (rs.lvstat == STAT_LOWER)
02927       {
02928         mpf_EGlpNumCopyDiff (alpha, lp->lz[lp->baz[rs.lindex]], lp->xbz[rs.lindex]);
02929         mpf_EGlpNumAddInnProdTo (lp->dobjval, rs.tz, alpha);
02930       }
02931       else
02932       {
02933         mpf_EGlpNumCopyDiff (alpha, lp->xbz[rs.lindex], lp->uz[lp->baz[rs.lindex]]);
02934         mpf_EGlpNumAddInnProdTo (lp->dobjval, rs.tz, alpha);
02935       }
02936       mpf_EGlpNumCopyFrac (alpha, lp->dz[eindex], rs.pivotval);
02937       mpf_EGlpNumCopy (lp->objval, lp->dobjval);
02938 
02939       mpf_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
02940       //ILL_IFTRACE("%s:e\n",__func__);
02941       mpf_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02942       mpf_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, mpf_zeroLpNum);
02943       mpf_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, mpf_zeroLpNum);
02944 
02945       if (pinf->dsinfo.norms && pinf->dII_price == QS_PRICE_DSTEEP)
02946       {
02947         mpf_ILLfct_compute_dsteep_upv (lp, &wz);
02948         rval = mpf_ILLprice_update_pricing_info (lp, pinf, DUAL_PHASEII, &wz,
02949                                              eindex, rs.lindex, rs.pivotval);
02950         CHECKRVALG (rval, CLEANUP);
02951       }
02952       else if (pinf->psinfo.norms && pinf->pII_price == QS_PRICE_PSTEEP)
02953       {
02954         mpf_ILLfct_compute_psteep_upv (lp, &wz);
02955         rval = mpf_ILLprice_update_pricing_info (lp, pinf, PRIMAL_PHASEII, &wz,
02956                                              eindex, rs.lindex, rs.pivotval);
02957         CHECKRVALG (rval, CLEANUP);
02958       }
02959 
02960       //ILL_IFTRACE("%s:e\n",__func__);
02961       mpf_ILLfct_update_xz (lp, rs.tz, eindex, rs.lindex);
02962       mpf_ILLfct_update_basis_info (lp, eindex, rs.lindex, rs.lvstat);
02963       rval = mpf_ILLbasis_update (lp, &updz, rs.lindex, &refactor, &singular);
02964       CHECKRVALG (rval, CLEANUP);
02965 
02966       if (singular)
02967       {
02968         fprintf (stderr, "singular matrix in pivot_in\n");
02969         rval = E_SIMPLEX_ERROR;
02970         ILL_CLEANUP;
02971       }
02972       if (!refactor)
02973       {
02974         mpf_ILLfct_update_piz (lp, alpha);
02975         mpf_ILLfct_update_dz (lp, eindex, alpha);
02976       }
02977       else
02978       {
02979         mpf_ILLfct_compute_xbz (lp);
02980         mpf_ILLfct_compute_piz (lp);
02981         mpf_ILLfct_compute_dz (lp);
02982         mpf_ILLfct_compute_dobj (lp);
02983       }
02984     }
02985   }
02986   /*
02987    * mpf_ILLfct_dphaseI_simple_update (lp, lp->tol->dfeas_tol);
02988    * mpf_ILLfct_compute_xbz (lp);
02989    * mpf_ILLfct_compute_dobj (lp);
02990    */
02991 
02992   mpf_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
02993   mpf_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
02994   mpf_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
02995 
02996 CLEANUP:
02997   if (pivot_opt == SIMPLEX_PIVOTINROW)
02998     ILL_IFFREE (clist, int);
02999 
03000   mpf_ILLsvector_free (&wz);
03001   mpf_ILLsvector_free (&updz);
03002   mpf_EGlpNumClearVar (alpha);
03003   mpf_EGlpNumClearVar (fi.totinfeas);
03004   mpf_EGlpNumClearVar (rs.tz);
03005   mpf_EGlpNumClearVar (rs.lbound);
03006   mpf_EGlpNumClearVar (rs.ecoeff);
03007   mpf_EGlpNumClearVar (rs.pivotval);
03008   EG_RETURN (rval);
03009 }
03010 
03011 static int mpf_report_value (
03012   mpf_lpinfo * lp,
03013   mpf_iter_info * it,
03014   const char *value_name,
03015   mpf_t value)
03016 {
03017   int rval = 0;
03018 
03019   if (it->sdisplay && it->itercnt % lp->iterskip == 0)
03020   {
03021     char buffer[1024];
03022 
03023     snprintf (buffer, (size_t) 1023, "(%d): %s = %10.7lf\n", it->itercnt,
03024               value_name, mpf_EGlpNumToLf (value));
03025     buffer[1022] = '\n';
03026     buffer[1023] = '\0';
03027     rval = ILLstring_report (buffer, &lp->O->reporter);
03028     fflush (stdout);
03029   }
03030   else
03031   {
03032     /* make sure ILLstring_report is called at least every 10 iterations */
03033     if (it->itercnt % (lp->iterskip / 10))
03034     {
03035       rval = ILLstring_report (NULL, &lp->O->reporter);
03036     }
03037   }
03038   if (rval != 0)
03039   {                             /* ILLstring_report was called and failed, which means we should abort */
03040     it->solstatus = QS_LP_ABORTED;
03041   }
03042   return rval;
03043 }
03044 
03045 #undef mpf_QSOPT_CURRENT_PRECICION

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