00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static int TRACE = 0;
00025
00026 #define ldbl_QSOPT_CURRENT_PRECICION
00027 #include "basicdefs.h"
00028 #include "config.h"
00029 #include "ldbl_iqsutil.h"
00030 #include "ldbl_lpdata.h"
00031 #include "ldbl_lpdefs.h"
00032
00033 #include "stddefs.h"
00034 #include "ldbl_fct.h"
00035 #include "ldbl_ratio.h"
00036 #include "ldbl_price.h"
00037 #include "ldbl_basis.h"
00038 #include "ldbl_simplex.h"
00039 #include "ldbl_dstruct.h"
00040 #include "ldbl_qstruct.h"
00041 #include "ldbl_qsopt.h"
00042 #include "ldbl_lib.h"
00043 #include "ldbl_lp.h"
00044
00045 #ifdef USEDMALLOC
00046 #include "dmalloc.h"
00047 #endif
00048
00049 static void ldbl_init_lp_status_info (
00050 ldbl_lp_status_info * ls),
00051 ldbl_init_simplex_tols (
00052 ldbl_lpinfo * lp),
00053 ldbl_monitor_iter (
00054 ldbl_lpinfo * lp,
00055 ldbl_price_info * p,
00056 ldbl_iter_info * it,
00057 int cphase),
00058 ldbl_get_current_stat (
00059 ldbl_lp_status_info * p,
00060 int algorithm,
00061 int *bstat);
00062
00063 static int ldbl_terminate_simplex (
00064 ldbl_lpinfo * lp,
00065 int phase,
00066 ldbl_iter_info * it),
00067 ldbl_primal_phaseI_step (
00068 ldbl_lpinfo * lp,
00069 ldbl_price_info * pinf,
00070 ldbl_svector * updz,
00071 ldbl_svector * wz,
00072 ldbl_iter_info * it),
00073 ldbl_primal_phaseII_step (
00074 ldbl_lpinfo * lp,
00075 ldbl_price_info * pinf,
00076 ldbl_svector * updz,
00077 ldbl_svector * wz,
00078 ldbl_iter_info * it),
00079 ldbl_dual_phaseI_step (
00080 ldbl_lpinfo * lp,
00081 ldbl_price_info * pinf,
00082 ldbl_svector * updz,
00083 ldbl_svector * wz,
00084 ldbl_iter_info * it),
00085 ldbl_dual_phaseII_step (
00086 ldbl_lpinfo * lp,
00087 ldbl_price_info * pinf,
00088 ldbl_svector * updz,
00089 ldbl_svector * wz,
00090 ldbl_iter_info * it),
00091 ldbl_report_value (
00092 ldbl_lpinfo * lp,
00093 ldbl_iter_info * it,
00094 const char *value_name,
00095 long double value);
00096
00097
00098 void ldbl_ILLsimplex_init_lpinfo (
00099 ldbl_lpinfo * lp)
00100 {
00101 ldbl_ILLbasis_init_basisinfo (lp);
00102 ldbl_init_internal_lpinfo (lp);
00103 }
00104
00105 void ldbl_ILLsimplex_free_lpinfo (
00106 ldbl_lpinfo * lp)
00107 {
00108 if (lp)
00109 {
00110 ldbl_EGlpNumFreeArray (lp->lz);
00111 ldbl_EGlpNumFreeArray (lp->uz);
00112 ldbl_EGlpNumFreeArray (lp->cz);
00113 ldbl_ILLbasis_free_basisinfo (lp);
00114 ldbl_free_internal_lpinfo (lp);
00115 }
00116 }
00117
00118 void ldbl_ILLsimplex_load_lpinfo (
00119 ldbl_ILLlpdata * qslp,
00120 ldbl_lpinfo * lp)
00121 {
00122 lp->basisid = -1;
00123 lp->maxiter = 500000;
00124 lp->maxtime = 300000;
00125
00126 lp->iterskip = 100;
00127 ldbl_EGlpNumCopy (lp->objbound, ldbl_INFTY);
00128 lp->O = qslp;
00129 }
00130
00131 void ldbl_ILLsimplex_set_bound (
00132 ldbl_lpinfo * lp,
00133 const long double * objbound,
00134 int sense)
00135 {
00136 ldbl_EGlpNumCopy (lp->objbound, *objbound);
00137 if (sense == ldbl_ILL_MAX)
00138 ldbl_EGlpNumSign (lp->objbound);
00139 }
00140
00141 static void ldbl_init_lp_status_info (
00142 ldbl_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 ldbl_init_simplex_tols (
00154 ldbl_lpinfo * lp)
00155 {
00156 ldbl_EGlpNumCopy (lp->tol->pfeas_tol, ldbl_PFEAS_TOLER);
00157 ldbl_EGlpNumCopy (lp->tol->dfeas_tol, ldbl_DFEAS_TOLER);
00158 ldbl_EGlpNumCopy (lp->tol->pivot_tol, ldbl_PIVOT_TOLER);
00159 ldbl_EGlpNumCopy (lp->tol->szero_tol, ldbl_SZERO_TOLER);
00160 ldbl_EGlpNumCopy (lp->tol->ip_tol, lp->tol->pfeas_tol);
00161 ldbl_EGlpNumCopy (lp->tol->id_tol, lp->tol->dfeas_tol);
00162 if (ldbl_EGlpNumIsNeqqZero (lp->tol->ip_tol))
00163 {
00164 #if VERBOSE_LEVEL <= DEBUG
00165 MESSAGE (VERBOSE_LEVEL, "ip_tol %lg", ldbl_EGlpNumToLf (lp->tol->ip_tol));
00166 MESSAGE (VERBOSE_LEVEL, "eps %lg", ldbl_EGlpNumToLf (ldbl_epsLpNum));
00167 MESSAGE (VERBOSE_LEVEL, "ldbl_PFEAS_TOLER %lg", ldbl_EGlpNumToLf (ldbl_PFEAS_TOLER));
00168 #endif
00169 ldbl_EGlpNumDivUiTo (lp->tol->ip_tol, 2UL);
00170 }
00171 if (ldbl_EGlpNumIsNeqqZero (lp->tol->id_tol))
00172 {
00173 #if VERBOSE_LEVEL <= DEBUG
00174 MESSAGE (VERBOSE_LEVEL, "id_tol %lg", ldbl_EGlpNumToLf (lp->tol->id_tol));
00175 #endif
00176 ldbl_EGlpNumDivUiTo (lp->tol->id_tol, 2UL);
00177 }
00178 }
00179
00180 void ldbl_init_internal_lpinfo (
00181 ldbl_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 ldbl_ILLsvector_init (&(lp->zz));
00214 ldbl_ILLsvector_init (&(lp->yjz));
00215 ldbl_ILLsvector_init (&(lp->zA));
00216 ldbl_ILLsvector_init (&(lp->work));
00217 ldbl_ILLsvector_init (&(lp->srhs));
00218 ldbl_ILLsvector_init (&(lp->ssoln));
00219 ILL_SAFE_MALLOC (lp->tol, 1, ldbl_tol_struct);
00220 ldbl_EGlpNumInitVar (lp->tol->pfeas_tol);
00221 ldbl_EGlpNumInitVar (lp->tol->dfeas_tol);
00222 ldbl_EGlpNumInitVar (lp->tol->pivot_tol);
00223 ldbl_EGlpNumInitVar (lp->tol->szero_tol);
00224 ldbl_EGlpNumInitVar (lp->tol->ip_tol);
00225 ldbl_EGlpNumInitVar (lp->tol->id_tol);
00226 ILL_SAFE_MALLOC (lp->cnts, 1, ldbl_count_struct);
00227 ldbl_EGlpNumInitVar (lp->cnts->y_ravg);
00228 ldbl_EGlpNumInitVar (lp->cnts->z_ravg);
00229 ldbl_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 ldbl_free_internal_lpinfo (
00239 ldbl_lpinfo * lp)
00240 {
00241 ldbl_bndinfo *binfo = 0;
00242 ldbl_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 ldbl_EGlpNumFreeArray (lp->rowval);
00251 lp->localrows = 0;
00252 }
00253 ldbl_EGlpNumFreeArray (lp->lz);
00254 ldbl_EGlpNumFreeArray (lp->uz);
00255 ldbl_EGlpNumFreeArray (lp->cz);
00256 ldbl_EGlpNumFreeArray (lp->xbz);
00257 ldbl_EGlpNumFreeArray (lp->piz);
00258 ldbl_EGlpNumFreeArray (lp->pIpiz);
00259 ldbl_EGlpNumFreeArray (lp->dz);
00260 ldbl_EGlpNumFreeArray (lp->pIdz);
00261 ldbl_EGlpNumFreeArray (lp->pIxbz);
00262
00263 ILL_IFFREE (lp->vtype, int);
00264 ILL_IFFREE (lp->vclass, char);
00265
00266 ldbl_ILLsvector_free (&(lp->zz));
00267 ldbl_ILLsvector_free (&(lp->yjz));
00268 ldbl_ILLsvector_free (&(lp->zA));
00269 ldbl_ILLsvector_free (&(lp->work));
00270 ldbl_ILLsvector_free (&(lp->srhs));
00271 ldbl_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 ldbl_EGlpNumFreeArray (lp->upd.t);
00277
00278 ILL_IFFREE (lp->bfeas, int);
00279 ILL_IFFREE (lp->dfeas, int);
00280
00281 if (lp->tol)
00282 {
00283 ldbl_EGlpNumClearVar (lp->tol->pfeas_tol);
00284 ldbl_EGlpNumClearVar (lp->tol->dfeas_tol);
00285 ldbl_EGlpNumClearVar (lp->tol->pivot_tol);
00286 ldbl_EGlpNumClearVar (lp->tol->szero_tol);
00287 ldbl_EGlpNumClearVar (lp->tol->ip_tol);
00288 ldbl_EGlpNumClearVar (lp->tol->id_tol);
00289 ILL_IFFREE (lp->tol, ldbl_tol_struct);
00290 }
00291 if (lp->cnts)
00292 {
00293 ldbl_EGlpNumClearVar (lp->cnts->y_ravg);
00294 ldbl_EGlpNumClearVar (lp->cnts->z_ravg);
00295 ldbl_EGlpNumClearVar (lp->cnts->za_ravg);
00296 ILL_IFFREE (lp->cnts, ldbl_count_struct);
00297 }
00298
00299 while (lp->bchanges)
00300 {
00301 binfo = lp->bchanges;
00302 ldbl_EGlpNumClearVar (binfo->pbound);
00303 ldbl_EGlpNumClearVar (binfo->cbound);
00304 lp->bchanges = binfo->next;
00305 ILL_IFFREE (binfo, ldbl_bndinfo);
00306 }
00307
00308 while (lp->cchanges)
00309 {
00310 cinfo = lp->cchanges;
00311 ldbl_EGlpNumClearVar (cinfo->pcoef);
00312 ldbl_EGlpNumClearVar (cinfo->ccoef);
00313 lp->cchanges = cinfo->next;
00314 ILL_IFFREE (cinfo, ldbl_coefinfo);
00315 }
00316 }
00317
00318 int ldbl_build_internal_lpinfo (
00319 ldbl_lpinfo * lp)
00320 {
00321 int rval = 0;
00322 int i, n;
00323 ldbl_ILLlpdata *qslp = lp->O;
00324 ldbl_ILLlp_sinfo *S = lp->O->sinfo;
00325 long double *lower, *upper, *obj;
00326 ldbl_ILLlp_rows lprows;
00327 ldbl_ILLmatrix *A;
00328
00329 ldbl_init_lp_status_info (&(lp->probstat));
00330 ldbl_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 = ldbl_EGlpNumAllocArray (lp->ncols);
00361 lp->cz = ldbl_EGlpNumAllocArray (lp->ncols);
00362 lp->uz = ldbl_EGlpNumAllocArray (lp->ncols);
00363 if (!lp->lz || !lp->uz || !lp->cz)
00364 {
00365 fprintf (stderr, "ldbl_build_internal_lpinfo\n");
00366 rval = 1;
00367 goto CLEANUP;
00368 }
00369 for (i = 0; i < lp->ncols; i++)
00370 {
00371 ldbl_EGlpNumCopy (lp->lz[i], lower[i]);
00372 ldbl_EGlpNumCopy (lp->uz[i], upper[i]);
00373 ldbl_EGlpNumCopy (lp->cz[i], obj[i]);
00374 if (qslp->objsense == ldbl_ILL_MAX)
00375 {
00376 ldbl_EGlpNumSign (lp->cz[i]);
00377 }
00378 }
00379
00380 if (!lp->O->rA)
00381 {
00382 rval = ldbl_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
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 = ldbl_EGlpNumAllocArray (lp->nrows);
00401 lp->piz = ldbl_EGlpNumAllocArray (lp->nrows);
00402 lp->dz = ldbl_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 = ldbl_ILLsvector_alloc (&(lp->zz), lp->nrows);
00410 CHECKRVALG (rval, CLEANUP);
00411 rval = ldbl_ILLsvector_alloc (&(lp->yjz), lp->nrows);
00412 CHECKRVALG (rval, CLEANUP);
00413 rval = ldbl_ILLsvector_alloc (&(lp->zA), lp->nnbasic);
00414 CHECKRVALG (rval, CLEANUP);
00415 rval = ldbl_ILLsvector_alloc (&(lp->work), lp->ncols);
00416 CHECKRVALG (rval, CLEANUP);
00417 rval = ldbl_ILLsvector_alloc (&(lp->srhs), lp->nrows);
00418 CHECKRVALG (rval, CLEANUP);
00419 rval = ldbl_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 ldbl_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 = ldbl_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 ldbl_init_simplex_tols (lp);
00439 ldbl_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 ldbl_free_internal_lpinfo (lp);
00454 EG_RETURN (rval);
00455 }
00456
00457 int ldbl_ILLsimplex_retest_psolution (
00458 ldbl_lpinfo * lp,
00459 ldbl_price_info * p,
00460 int phase,
00461 ldbl_feas_info * fi)
00462 {
00463 int rval = 0;
00464 int fbid = lp->fbasisid;
00465 int bid = lp->basisid;
00466 long double *ptol = &(lp->tol->pfeas_tol);
00467 long double *dtol = &(lp->tol->dfeas_tol);
00468 long double *iptol = &(lp->tol->ip_tol);
00469 long double *idtol = &(lp->tol->id_tol);
00470
00471 fi->pstatus = -1;
00472 fi->dstatus = -1;
00473 if (fbid < bid - PARAM_PRIMAL_REFACTORGAP)
00474 {
00475 rval = ldbl_ILLbasis_refactor (lp);
00476 CHECKRVALG (rval, CLEANUP);
00477 }
00478 if (fbid < bid - PARAM_PRIMAL_RESOLVEGAP)
00479 ldbl_ILLfct_compute_xbz (lp);
00480
00481 if (phase == PRIMAL_PHASEII)
00482 {
00483 if (fbid < bid - PARAM_PRIMAL_RESOLVEGAP)
00484 {
00485 ldbl_ILLfct_compute_piz (lp);
00486 ldbl_ILLfct_compute_dz (lp);
00487 if (p != NULL && p->p_strategy == COMPLETE_PRICING)
00488 ldbl_ILLprice_compute_dual_inf (lp, p, NULL, 0, PRIMAL_PHASEII);
00489 }
00490 ldbl_ILLfct_compute_pobj (lp);
00491 ldbl_ILLfct_check_pfeasible (lp, fi, *ptol);
00492 ldbl_ILLfct_check_dfeasible (lp, fi, *dtol);
00493 }
00494 else if (phase == PRIMAL_PHASEI)
00495 {
00496 ldbl_ILLfct_check_pfeasible (lp, fi, *iptol);
00497 if (fi->pstatus != PRIMAL_FEASIBLE)
00498 {
00499 if (lp->pIpiz)
00500 {
00501 ldbl_ILLfct_compute_phaseI_piz (lp);
00502 ldbl_ILLfct_compute_phaseI_dz (lp);
00503 ldbl_ILLfct_check_pIdfeasible (lp, fi, *idtol);
00504 if (p != NULL && p->p_strategy == COMPLETE_PRICING)
00505 ldbl_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 ldbl_precision");
00513 return rval;
00514 }
00515 EG_RETURN (rval);
00516 }
00517
00518 int ldbl_ILLsimplex_retest_dsolution (
00519 ldbl_lpinfo * lp,
00520 ldbl_price_info * p,
00521 int phase,
00522 ldbl_feas_info * fi)
00523 {
00524 int rval = 0;
00525 int fbid = lp->fbasisid;
00526 int bid = lp->basisid;
00527 long double *ptol = &(lp->tol->pfeas_tol);
00528 long double *dtol = &(lp->tol->dfeas_tol);
00529 long double *iptol = &(lp->tol->ip_tol);
00530 long double *idtol = &(lp->tol->id_tol);
00531
00532
00533
00534 fi->pstatus = -1;
00535 fi->dstatus = -1;
00536 if (fbid < bid - PARAM_DUAL_REFACTORGAP)
00537 {
00538
00539 rval = ldbl_ILLbasis_refactor (lp);
00540 CHECKRVALG (rval, CLEANUP);
00541 }
00542 if (fbid < bid - PARAM_DUAL_RESOLVEGAP)
00543 {
00544
00545 ldbl_ILLfct_compute_piz (lp);
00546 ldbl_ILLfct_compute_dz (lp);
00547 }
00548
00549 if (phase == DUAL_PHASEII)
00550 {
00551
00552 if (fbid < bid - PARAM_DUAL_RESOLVEGAP)
00553 {
00554
00555 ldbl_ILLfct_compute_xbz (lp);
00556 CHECKRVALG (rval, CLEANUP);
00557 if (p != NULL)
00558 {
00559
00560 if (p->d_strategy == COMPLETE_PRICING)
00561 {
00562
00563 ldbl_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEII);
00564 }
00565 else
00566 {
00567
00568 ldbl_ILLprice_update_mpartial_price (lp, p, DUAL_PHASEII, ROW_PRICING);
00569 }
00570 }
00571 }
00572
00573 ldbl_ILLfct_compute_dobj (lp);
00574 ldbl_ILLfct_check_dfeasible (lp, fi, *dtol);
00575 ldbl_ILLfct_check_pfeasible (lp, fi, *ptol);
00576 }
00577 else if (phase == DUAL_PHASEI)
00578 {
00579 ldbl_ILLfct_check_dfeasible (lp, fi, *idtol);
00580 if (fi->dstatus != DUAL_FEASIBLE)
00581 {
00582 ldbl_ILLfct_compute_phaseI_xbz (lp);
00583 ldbl_ILLfct_check_pIpfeasible (lp, fi, *iptol);
00584 if (p != NULL)
00585 {
00586 if (p->d_strategy == COMPLETE_PRICING)
00587 ldbl_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEI);
00588 else
00589 ldbl_ILLprice_update_mpartial_price (lp, p, DUAL_PHASEI, ROW_PRICING);
00590 }
00591 }
00592 }
00593 CLEANUP:
00594 EG_RETURN (rval);
00595 }
00596
00597 int ldbl_ILLsimplex_solution (
00598 ldbl_lpinfo * lp,
00599 long double * xz,
00600 long double * piz,
00601 long double * dz,
00602 long double * 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 ldbl_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 ldbl_EGlpNumCopy (xz[col], lp->uz[col]);
00620 else if (lp->vstat[col] == STAT_LOWER)
00621 ldbl_EGlpNumCopy (xz[col], lp->lz[col]);
00622 else
00623 ldbl_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 ldbl_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 ldbl_EGlpNumZero (dz[lp->baz[i]]);
00643 for (j = 0; j < lp->nnbasic; j++)
00644 ldbl_EGlpNumCopy (dz[lp->nbaz[j]], lp->dz[j]);
00645 }
00646 if (objval != NULL)
00647 ldbl_EGlpNumCopy (*objval, lp->objval);
00648 return 0;
00649 }
00650
00651 int ldbl_ILLsimplex_infcertificate (
00652 ldbl_lpinfo * lp,
00653 long double * pi)
00654 {
00655 int i, col, nz;
00656 char *sense;
00657 long double *x, *l, *u;
00658 ldbl_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 ldbl_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 ldbl_EGlpNumZero (pi[i]);
00683
00684 if (ldbl_EGlpNumIsNeqq (*l, ldbl_NINFTY) && ldbl_EGlpNumIsLess (*x, *l))
00685 {
00686 for (i = 0, nz = lp->zz.nzcnt; i < nz; i++)
00687 ldbl_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 ldbl_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' && ldbl_EGlpNumIsLessZero (pi[i]))
00705 ldbl_EGlpNumZero (pi[i]);
00706 if (sense[i] == 'L' && ldbl_EGlpNumIsGreatZero (pi[i]))
00707 ldbl_EGlpNumZero (pi[i]);
00708 }
00709 return 0;
00710 }
00711
00712 #if SIMPLEX_DEBUG > 1
00713 static void ldbl_test_cert (
00714 ldbl_lpinfo * lp,
00715 long double * pi)
00716 {
00717 int i, j;
00718 int mcnt, mbeg;
00719 long double fsum, sum;
00720
00721 ldbl_EGlpNumInitVar (fsum);
00722 ldbl_EGlpNumInitVar (sum);
00723 ldbl_EGlpNumZero (fsum);
00724
00725 for (i = 0; i < lp->nrows; i++)
00726 {
00727 if (lp->O->sense[i] == 'G' && ldbl_EGlpNumIsLessZero (pi[i]))
00728 printf ("compl \n");
00729 if (lp->O->sense[i] == 'L' && ldbl_EGlpNumIsGreatZero (pi[i]))
00730 printf ("compll \n");
00731 }
00732
00733 for (i = 0; i < lp->nrows; i++)
00734 ldbl_EGlpNumAddInnProdTo (fsum, pi[i], lp->bz[i]);
00735
00736 for (j = 0; j < lp->nnbasic; j++)
00737 {
00738 ldbl_EGlpNumZero (sum);
00739 mcnt = lp->matcnt[j];
00740 mbeg = lp->matbeg[j];
00741 for (i = 0; i < mcnt; i++)
00742 ldbl_EGlpNumAddInnProdTo (sum, pi[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00743
00744 if (ldbl_EGlpNumIsLess (ldbl_PFEAS_TOLER, sum) &&
00745 (lp->vtype[j] == VLOWER || lp->vtype[j] == VFREE))
00746 printf ("compl2\n");
00747 else
00748 {
00749 ldbl_EGlpNumSign (sum);
00750 if (ldbl_EGlpNumIsLess (ldbl_PFEAS_TOLER, sum) &&
00751 (lp->vtype[j] == VUPPER || lp->vtype[j] == VFREE))
00752 printf ("compl1\n");
00753 ldbl_EGlpNumSign (sum);
00754 }
00755
00756 if (ldbl_EGlpNumIsLessZero (sum)
00757 && (lp->vtype[j] & (VFREE | VUPPER)) == 0)
00758 ldbl_EGlpNumSubInnProdTo (fsum, sum, lp->lz[j]);
00759 else if (ldbl_EGlpNumIsGreatZero (sum)
00760 && (lp->vtype[j] & (VFREE | VLOWER)) == 0)
00761 ldbl_EGlpNumSubInnProdTo (fsum, sum, lp->uz[j]);
00762 }
00763 printf ("fsum = %.8f\n", ldbl_EGlpNumToLf (fsum));
00764 ldbl_EGlpNumClearVar (fsum);
00765 ldbl_EGlpNumClearVar (sum);
00766 }
00767 #endif
00768
00769 static void ldbl_save_paraminfo (
00770 ldbl_price_info * pinf,
00771 ldbl_iter_info * it)
00772 {
00773 ldbl_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 ldbl_restore_paraminfo (
00785 ldbl_iter_info * it,
00786 ldbl_price_info * pinf)
00787 {
00788 ldbl_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 ldbl_ILLsimplex (
00800 ldbl_lpinfo * lp,
00801 int algorithm,
00802 ldbl_ILLlp_basis * B,
00803 ldbl_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 ldbl_svector wz;
00813 ldbl_svector updz;
00814 ldbl_feas_info fi;
00815 ldbl_iter_info it;
00816
00817 ldbl_EGlpNumInitVar (fi.totinfeas);
00818 ldbl_EGlpNumInitVar (it.prevobj);
00819 ldbl_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 ldbl_EGlpNumCopy (it.prevobj, ldbl_INFTY);
00832 it.nosolve = 0;
00833 it.noprog = 0;
00834 ldbl_EGlpNumCopy (it.objtol, ldbl_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 ldbl_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 ldbl_free_internal_lpinfo (lp);
00850 ldbl_init_internal_lpinfo (lp);
00851 rval = ldbl_build_internal_lpinfo (lp);
00852 CHECKRVALG (rval, CLEANUP);
00853
00854 ldbl_ILLsvector_init (&wz);
00855 rval = ldbl_ILLsvector_alloc (&wz, lp->nrows);
00856 CHECKRVALG (rval, CLEANUP);
00857 ldbl_ILLsvector_init (&updz);
00858 rval = ldbl_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 ldbl_ILLsimplex on %s...\n", lp->O->probname);
00870
00871
00872
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 ldbl_ILLfct_set_variable_type (lp);
00879
00880 if (B != 0)
00881 {
00882 rval = ldbl_ILLbasis_load (lp, B);
00883 CHECKRVALG (rval, CLEANUP);
00884 if (it.algorithm == DUAL_SIMPLEX)
00885 {
00886 if (B->rownorms)
00887 {
00888 rval = ldbl_ILLprice_load_rownorms (lp, B->rownorms, pinf);
00889 CHECKRVALG (rval, CLEANUP);
00890 }
00891 else
00892 ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);
00893 }
00894 else if (it.algorithm == PRIMAL_SIMPLEX)
00895 {
00896 if (B->colnorms)
00897 {
00898 rval = ldbl_ILLprice_load_colnorms (lp, B->colnorms, pinf);
00899 CHECKRVALG (rval, CLEANUP);
00900 }
00901 else
00902 ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
00903 }
00904 else if (it.algorithm != PRIMAL_OR_DUAL)
00905 {
00906 fprintf (stderr, "Unknown algorithm %d in ldbl_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 = ldbl_ILLbasis_get_initial (lp, it.algorithm);
00915 else
00916 rval = ldbl_ILLbasis_get_cinitial (lp, it.algorithm);
00917 CHECKRVALG (rval, CLEANUP);
00918 ldbl_ILLprice_free_pricing_info (pinf);
00919 }
00920
00921 if (lp->fbasisid != lp->basisid)
00922 {
00923 rval = ldbl_ILLbasis_factor (lp, &singular);
00924 CHECKRVALG (rval, CLEANUP);
00925 if (singular)
00926 {
00927 MESSAGE (__QS_SB_VERB, "Singular basis found!");
00928 ldbl_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 ldbl_init_lp_status_info (&(lp->basisstat));
00945
00946 ldbl_ILLfct_compute_piz (lp);
00947 ldbl_ILLfct_compute_dz (lp);
00948 if (it.algorithm == DUAL_SIMPLEX)
00949 {
00950 if (B != NULL || it.resumeid == SIMPLEX_RESUME_UNSHIFT)
00951 ldbl_ILLfct_dual_adjust (lp, lp->tol->dfeas_tol);
00952 else
00953 ldbl_ILLfct_dual_adjust (lp, ldbl_zeroLpNum);
00954 }
00955 ldbl_ILLfct_compute_xbz (lp);
00956
00957 ldbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
00958 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
00959 ldbl_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, PHASEII, PHASEII);
00960
00961 if (fi.dstatus == DUAL_FEASIBLE && ldbl_EGlpNumIsNeqq (lp->objbound, ldbl_INFTY))
00962 {
00963 ldbl_ILLfct_compute_dobj (lp);
00964 if (ldbl_EGlpNumIsLess (lp->objbound, lp->dobjval))
00965 {
00966 it.solstatus = ILL_BND_REACHED;
00967 printf ("solstatus = ILL_BND_REACHED = 5 %lf %lf\n",
00968 ldbl_EGlpNumToLf (lp->objbound), ldbl_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 ldbl_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 (ldbl_EGlpNumToLf (lp->pinfeas) < 10 * ldbl_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 = ldbl_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
01015
01016 if (phase == PRIMAL_PHASEI)
01017 {
01018 rval = ldbl_primal_phaseI_step (lp, pinf, &updz, &wz, &it);
01019 CHECKRVALG (rval, CLEANUP);
01020 }
01021
01022 else if (phase == PRIMAL_PHASEII)
01023 {
01024 rval = ldbl_primal_phaseII_step (lp, pinf, &updz, &wz, &it);
01025 CHECKRVALG (rval, CLEANUP);
01026 }
01027
01028 else if (phase == DUAL_PHASEI)
01029 {
01030 rval = ldbl_dual_phaseI_step (lp, pinf, &updz, &wz, &it);
01031 CHECKRVALG (rval, CLEANUP);
01032 }
01033
01034 else if (phase == DUAL_PHASEII)
01035 {
01036 rval = ldbl_dual_phaseII_step (lp, pinf, &updz, &wz, &it);
01037 CHECKRVALG (rval, CLEANUP);
01038 }
01039 if (it.nextstep == SIMPLEX_RESUME)
01040 {
01041
01042 ldbl_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 ldbl_ILLfct_unroll_bound_change (lp);
01059 ldbl_ILLfct_unroll_coef_change (lp);
01060
01061 rval = ldbl_ILLbasis_get_initial (lp, it.algorithm);
01062 CHECKRVALG (rval, CLEANUP);
01063 rval = ldbl_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
01079 it.itercnt++;
01080
01081 if (it.nextphase != phase)
01082 {
01083
01084 it.newphase = SIMPLEX_PHASE_NEW;
01085 phase = it.nextphase;
01086 new_price = ldbl_ILLprice_get_price (pinf, phase);
01087
01088 if (pinf->cur_price != new_price)
01089 {
01090
01091 ldbl_ILLprice_free_pricing_info (pinf);
01092 rval = ldbl_ILLprice_build_pricing_info (lp, pinf, phase);
01093 CHECKRVALG (rval, CLEANUP);
01094 }
01095 }
01096 }
01097 }
01098
01099 #if SIMPLEX_DEBUG > 0
01100 ldbl_ILLfct_print_counts (lp);
01101 #endif
01102 LIMIT_TERMINATE:
01103 rval = ldbl_terminate_simplex (lp, phase, &it);
01104 CHECKRVALG (rval, CLEANUP);
01105
01106 TERMINATE:
01107 ldbl_restore_paraminfo (&it, pinf);
01108
01109 if (it.sdisplay)
01110 {
01111 printf ("completed ldbl_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 ldbl_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 long double *pi = NULL;
01180
01181 pi = ldbl_EGlpNumAllocArray (lp->nrows);
01182 rva = ldbl_ILLsimplex_infcertificate (lp, pi);
01183 printf ("rva = %d\n", rva);
01184 if (!rva)
01185 {
01186 ldbl_test_cert (lp, pi);
01187 }
01188 ldbl_EGlpNumFreeArray (pi);
01189 }
01190 #endif
01191
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
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 ldbl_get_current_stat (&(lp->basisstat), it.algorithm, &bstat);
01208 switch (bstat)
01209 {
01210 case OPTIMAL:
01211 printf ("opt = %f\n", ldbl_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", ldbl_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", ldbl_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 ldbl_ILLsvector_free (&wz);
01244 ldbl_ILLsvector_free (&updz);
01245 ldbl_EGlpNumClearVar (it.prevobj);
01246 ldbl_EGlpNumClearVar (it.objtol);
01247 ldbl_EGlpNumClearVar (fi.totinfeas);
01248 if (rval == QS_LP_CHANGE_PREC)
01249 {
01250 MESSAGE (__QS_SB_VERB, "Changing ldbl_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 ldbl_terminate_simplex (
01261 ldbl_lpinfo * lp,
01262 int phase,
01263 ldbl_iter_info * it)
01264 {
01265 int rval = 0;
01266 int sphase;
01267 ldbl_feas_info fi;
01268
01269 ldbl_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 ldbl_ILLfct_unroll_bound_change (lp);
01284 }
01285 rval = ldbl_ILLsimplex_retest_psolution (lp, NULL, phase, &fi);
01286 CHECKRVALG (rval, CLEANUP);
01287
01288 sphase = (phase == PRIMAL_PHASEI) ? PHASEI : PHASEII;
01289 ldbl_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 ldbl_ILLfct_unroll_coef_change (lp);
01301 }
01302 rval = ldbl_ILLsimplex_retest_dsolution (lp, NULL, phase, &fi);
01303 CHECKRVALG (rval, CLEANUP);
01304
01305 sphase = (phase == DUAL_PHASEI) ? PHASEI : PHASEII;
01306 ldbl_ILLfct_set_status_values (lp, fi.pstatus, fi.dstatus, sphase, PHASEII);
01307 }
01308
01309 CLEANUP:
01310 ldbl_EGlpNumClearVar (fi.totinfeas);
01311 EG_RETURN (rval);
01312 }
01313
01314 static int ldbl_test_progress (
01315 long double objval,
01316 long double prevobj)
01317 {
01318 long double denom;
01319
01320 ldbl_EGlpNumInitVar (denom);
01321 ldbl_EGlpNumCopyDiff (denom, objval, prevobj);
01322 if (ldbl_EGlpNumIsNeqZero (objval, ldbl_PROGRESS_ZERO))
01323 ldbl_EGlpNumDivTo (denom, objval);
01324 if (!ldbl_EGlpNumIsNeqZero (denom, ldbl_PROGRESS_THRESH))
01325 {
01326 ldbl_EGlpNumClearVar (denom);
01327 return 0;
01328 }
01329 else
01330 {
01331 ldbl_EGlpNumClearVar (denom);
01332 return 1;
01333 }
01334 }
01335
01336 static void ldbl_monitor_iter (
01337 ldbl_lpinfo * lp,
01338 ldbl_price_info * p,
01339 ldbl_iter_info * it,
01340 int phase)
01341 {
01342 long double print_val;
01343 double tottime = ILLutil_zeit () - lp->starttime;
01344 int curtime = ldbl_ILLutil_our_floor (tottime);
01345 char print_str[20];
01346 ldbl_feas_info fi;
01347 int aborted = 0;
01348
01349 ldbl_EGlpNumInitVar (print_val);
01350 ldbl_EGlpNumInitVar (fi.totinfeas);
01351 ldbl_EGlpNumZero (fi.totinfeas);
01352 ldbl_EGlpNumZero (print_val);
01353
01354
01355 switch (phase)
01356 {
01357 case PRIMAL_PHASEI:
01358 ldbl_EGlpNumCopy (print_val, lp->pinfeas);
01359 ldbl_EGlpNumAddTo (fi.totinfeas, lp->pinfeas);
01360 strcpy (print_str, "primal infeas");
01361 if (ldbl_EGlpNumIsLessZero (lp->pinfeas) &&
01362 (ldbl_EGlpNumIsNeqZero (lp->pinfeas, ldbl_oneLpNum)))
01363 {
01364
01365
01366
01367
01368 }
01369 break;
01370 case PRIMAL_PHASEII:
01371 ldbl_EGlpNumCopy (print_val, lp->pobjval);
01372 strcpy (print_str, "primal objval");
01373 break;
01374 case DUAL_PHASEI:
01375 ldbl_EGlpNumCopy (print_val, lp->dinfeas);
01376 ldbl_EGlpNumAddTo (fi.totinfeas, lp->dinfeas);
01377 strcpy (print_str, "dual infeas");
01378 break;
01379 case DUAL_PHASEII:
01380 ldbl_EGlpNumCopy (print_val, lp->dobjval);
01381 strcpy (print_str, "dual objval");
01382 break;
01383 }
01384
01385 aborted = ldbl_report_value (lp, it, print_str, print_val);
01386
01387
01388
01389
01390 if (curtime != it->curtime)
01391 {
01392 it->curtime = curtime;
01393
01394
01395
01396
01397
01398
01399
01400 }
01401
01402 ldbl_EGlpNumAddUiTo (fi.totinfeas, 1000);
01403 if (ldbl_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 && ldbl_EGlpNumIsNeqq (lp->objbound, ldbl_INFTY))
01413 {
01414
01415 ldbl_EGlpNumCopyDiff (print_val, lp->dobjval, lp->objbound);
01416 if (ldbl_EGlpNumIsLess (it->objtol, print_val))
01417 {
01418 ldbl_ILLfct_unroll_coef_change (lp);
01419 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01420 ldbl_ILLfct_set_status_values (lp, -1, fi.dstatus, -1, PHASEII);
01421
01422 if (fi.dstatus == DUAL_FEASIBLE)
01423 {
01424 ldbl_ILLfct_compute_dobj (lp);
01425 if (ldbl_EGlpNumIsLess (lp->objbound, lp->dobjval))
01426 {
01427 it->solstatus = ILL_BND_REACHED;
01428 it->nextstep = SIMPLEX_TERMINATE;
01429
01430 {
01431 printf ("bound reached %lf %lf\n", ldbl_EGlpNumToLf (lp->objbound),
01432 ldbl_EGlpNumToLf (lp->dobjval));
01433 fflush (stdout);
01434 }
01435 }
01436 else
01437 ldbl_EGlpNumMultUiTo (it->objtol, 10);
01438 }
01439 else
01440 {
01441 it->nextphase = DUAL_PHASEI;
01442 it->newphase = SIMPLEX_PHASE_NEW;
01443 ldbl_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
01481 if(0){
01482 if (it->rounds && it->inner){
01483 it->inner --;
01484 if (it->inner == 0){
01485 printf ("restoring ..\n");
01486 ldbl_restore_paraminfo (it, p);
01487 it->newphase = SIMPLEX_PHASE_NEW;
01488 it->nextstep = SIMPLEX_RESUME;
01489
01490 ILL_CLEANUP;
01491 }
01492 }
01493 }
01494 if (phase == DUAL_PHASEII)
01495 {
01496 if (it->noprog > it->chkobj)
01497 {
01498 ldbl_ILLfct_perturb_coefs (lp);
01499 it->noprog = 0;
01500 ldbl_EGlpNumCopy (it->prevobj, lp->dobjval);
01501 }
01502 }
01503 else if (phase == PRIMAL_PHASEII)
01504 {
01505 if (it->noprog > it->chkobj)
01506 {
01507 ldbl_ILLfct_perturb_bounds (lp);
01508 it->noprog = 0;
01509 ldbl_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
01520 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01521 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01522 it->n_restart++;
01523
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
01534 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01535 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01536 it->n_restart++;
01537
01538 }
01539 }
01540 CLEANUP:
01541 ldbl_EGlpNumClearVar (fi.totinfeas);
01542 ldbl_EGlpNumClearVar (print_val);
01543 return;
01544 }
01545
01546 static int ldbl_primal_phaseI_step (
01547 ldbl_lpinfo * lp,
01548 ldbl_price_info * pinf,
01549 ldbl_svector * updz,
01550 ldbl_svector * wz,
01551 ldbl_iter_info * it)
01552 {
01553 int rval = 0;
01554 int singular = 0;
01555 int refactor = 0;
01556 int cphase = PRIMAL_PHASEI;
01557 long double alpha;
01558 ldbl_feas_info fi;
01559 ldbl_ratio_res rs;
01560 ldbl_price_res pr;
01561
01562 ldbl_EGlpNumInitVar (alpha);
01563 ldbl_EGlpNumInitVar (fi.totinfeas);
01564 ldbl_EGlpNumInitVar (pr.dinfeas);
01565 ldbl_EGlpNumInitVar (pr.pinfeas);
01566 ldbl_EGlpNumInitVar (rs.tz);
01567 ldbl_EGlpNumInitVar (rs.lbound);
01568 ldbl_EGlpNumInitVar (rs.ecoeff);
01569 ldbl_EGlpNumInitVar (rs.pivotval);
01570 ldbl_EGlpNumZero (alpha);
01571
01572 ldbl_ILLfct_update_counts (lp, CNT_PPHASE1ITER, 0, ldbl_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 ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->pinfeas);
01593 lp->pIpiz = ldbl_EGlpNumAllocArray (lp->nrows);
01594 lp->pIdz = ldbl_EGlpNumAllocArray (lp->nnbasic);
01595
01596 ldbl_ILLfct_compute_phaseI_piz (lp);
01597 if (pinf->p_strategy == COMPLETE_PRICING)
01598 {
01599 ldbl_ILLfct_compute_phaseI_dz (lp);
01600 #if USEHEAP > 0
01601 ldbl_ILLprice_free_heap (pinf);
01602 #endif
01603 ldbl_ILLprice_compute_dual_inf (lp, pinf, NULL, 0, PRIMAL_PHASEI);
01604 #if USEHEAP > 0
01605 rval = ldbl_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 ldbl_ILLprice_init_mpartial_price (lp, pinf, cphase, COL_PRICING);
01612 }
01613
01614 ldbl_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 ldbl_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 = ldbl_ILLsimplex_retest_psolution (lp, pinf, cphase, &fi);
01631
01632 CHECKRVALG (rval, CLEANUP);
01633 ldbl_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 ldbl_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[pr.eindex]);
01648 ldbl_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, ldbl_zeroLpNum);
01649 ldbl_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, ldbl_zeroLpNum);
01650
01651 ldbl_ILLratio_pI_test (lp, pr.eindex, pr.dir, &rs);
01652
01653
01654 if (rs.ratio_stat == RATIO_FAILED)
01655 {
01656
01657
01658
01659
01660
01661 it->algorithm = DUAL_SIMPLEX;
01662 it->nextstep = SIMPLEX_RESUME;
01663 it->resumeid = SIMPLEX_RESUME_NUMER;
01664
01665 it->n_restart++;
01666 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01667 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01668
01669 ILL_CLEANUP;
01670 }
01671 else if (rs.ratio_stat == RATIO_NEGATIVE)
01672 {
01673 long double itol;
01674
01675 ldbl_EGlpNumInitVar (itol);
01676
01677 ldbl_EGlpNumCopy (itol, lp->tol->ip_tol);
01678 ldbl_EGlpNumZero (lp->tol->ip_tol);
01679 ldbl_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01680 if (!ldbl_test_progress (lp->pinfeas, it->prevobj))
01681 it->noprog++;
01682 else
01683 {
01684 ldbl_EGlpNumCopy (it->prevobj, lp->pinfeas);
01685 it->noprog = 0;
01686 }
01687 ldbl_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01688 ldbl_EGlpNumCopy (lp->tol->ip_tol, itol);
01689 ldbl_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01690 ldbl_ILLfct_update_ppI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.eindex,
01691 rs.lindex, ldbl_zeroLpNum);
01692 ldbl_EGlpNumClearVar (itol);
01693 }
01694 else if (rs.ratio_stat == RATIO_NOBCHANGE)
01695 {
01696
01697 ldbl_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01698 if (!ldbl_test_progress (lp->pinfeas, it->prevobj))
01699 it->noprog++;
01700 else
01701 {
01702 ldbl_EGlpNumCopy (it->prevobj, lp->pinfeas);
01703 it->noprog = 0;
01704 }
01705
01706
01707 ldbl_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
01708 ldbl_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01709 ldbl_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01710 ldbl_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
01711 #if DENSE_PI > 0
01712 ldbl_fct_test_workvector (lp);
01713 ldbl_fct_test_pfeasible (lp);
01714 #endif
01715 ldbl_ILLfct_update_ppI_prices (lp, pinf, &(lp->srhs), &(lp->ssoln), pr.eindex,
01716 rs.lindex, ldbl_zeroLpNum);
01717 }
01718 else if (rs.ratio_stat == RATIO_BCHANGE)
01719 {
01720
01721 ldbl_EGlpNumCopyFrac (alpha, lp->pIdz[pr.eindex], rs.pivotval);
01722 ldbl_EGlpNumAddTo (lp->pinfeas, lp->upd.c_obj);
01723
01724 if (!ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->pinfeas);
01738 it->noprog = 0;
01739 }
01740
01741 ldbl_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
01742 ldbl_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, ldbl_zeroLpNum);
01743 if (pinf->p_strategy == COMPLETE_PRICING)
01744 {
01745
01746 ldbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
01747 ldbl_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, ldbl_zeroLpNum);
01748
01749 if (pinf->pI_price == QS_PRICE_PSTEEP)
01750 {
01751 ldbl_ILLfct_compute_psteep_upv (lp, wz);
01752 }
01753 }
01754
01755 rval =
01756 ldbl_ILLprice_update_pricing_info (lp, pinf, cphase, wz, pr.eindex, rs.lindex,
01757 rs.pivotval);
01758 CHECKRVALG (rval, CLEANUP);
01759
01760
01761 ldbl_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
01762 ldbl_ILLfct_update_pfeas (lp, rs.lindex, &(lp->srhs));
01763
01764 ldbl_ILLfct_compute_ppIzz (lp, &(lp->srhs), &(lp->ssoln));
01765 ldbl_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
01766 #if DENSE_PI > 0
01767 ldbl_fct_test_workvector (lp);
01768 ldbl_fct_test_pfeasible (lp);
01769 #endif
01770 rval = ldbl_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
01778 it->n_restart++;
01779 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01780 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01781
01782 ILL_CLEANUP;
01783 }
01784 if (!refactor)
01785 {
01786 ldbl_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 ldbl_ILLfct_compute_xbz (lp);
01792 ldbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->ip_tol);
01793 ldbl_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 ldbl_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 ldbl_EGlpNumFreeArray (lp->pIpiz);
01812 ldbl_EGlpNumFreeArray (lp->pIdz);
01813 }
01814 ldbl_EGlpNumClearVar (alpha);
01815 ldbl_EGlpNumClearVar (fi.totinfeas);
01816 ldbl_EGlpNumClearVar (pr.dinfeas);
01817 ldbl_EGlpNumClearVar (pr.pinfeas);
01818 ldbl_EGlpNumClearVar (rs.tz);
01819 ldbl_EGlpNumClearVar (rs.lbound);
01820 ldbl_EGlpNumClearVar (rs.ecoeff);
01821 ldbl_EGlpNumClearVar (rs.pivotval);
01822
01823 return rval;
01824 }
01825
01826 static int ldbl_primal_phaseII_step (
01827 ldbl_lpinfo * lp,
01828 ldbl_price_info * pinf,
01829 ldbl_svector * updz,
01830 ldbl_svector * wz,
01831 ldbl_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 long double lbound;
01841 long double alpha;
01842 ldbl_feas_info fi;
01843 ldbl_ratio_res rs;
01844 ldbl_price_res pr;
01845
01846 ldbl_EGlpNumInitVar (alpha);
01847 ldbl_EGlpNumInitVar (lbound);
01848 ldbl_EGlpNumInitVar (fi.totinfeas);
01849 ldbl_EGlpNumInitVar (pr.dinfeas);
01850 ldbl_EGlpNumInitVar (pr.pinfeas);
01851 ldbl_EGlpNumInitVar (rs.tz);
01852 ldbl_EGlpNumInitVar (rs.lbound);
01853 ldbl_EGlpNumInitVar (rs.ecoeff);
01854 ldbl_EGlpNumInitVar (rs.pivotval);
01855
01856 ldbl_ILLfct_update_counts (lp, CNT_PPHASE2ITER, 0, ldbl_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 ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->pobjval);
01877 ldbl_ILLfct_compute_piz (lp);
01878 if (pinf->p_strategy == COMPLETE_PRICING)
01879 {
01880 ldbl_ILLfct_compute_dz (lp);
01881 #if USEHEAP > 0
01882 ldbl_ILLprice_free_heap (pinf);
01883 #endif
01884 ldbl_ILLprice_compute_dual_inf (lp, pinf, NULL, 0, PRIMAL_PHASEII);
01885 #if USEHEAP > 0
01886 rval = ldbl_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 ldbl_ILLprice_init_mpartial_price (lp, pinf, cphase, COL_PRICING);
01894 }
01895 }
01896
01897 ldbl_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 ldbl_ILLprice_primal (lp, pinf, &pr, cphase);
01903
01904 if (pr.price_stat == PRICE_OPTIMAL)
01905 {
01906
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 ldbl_ILLfct_unroll_bound_change (lp);
01915 ldbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
01916 ldbl_ILLfct_set_status_values (lp, fi.pstatus, -1, PHASEII, -1);
01917
01918 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01919
01920
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
01929
01930 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
01931 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
01932 it->n_restart++;
01933 ILL_CLEANUP;
01934
01935
01936
01937
01938
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 ldbl_EGlpNumToLf (lp->pobjval));
01948 fflush (stdout);
01949 }
01950 rval = ldbl_ILLsimplex_retest_psolution (lp, pinf, cphase, &fi);
01951 CHECKRVALG (rval, CLEANUP);
01952 ldbl_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 ldbl_EGlpNumDivUiTo (lp->tol->ip_tol, 5);
01958 ldbl_EGlpNumDivUiTo (lp->tol->id_tol, 5);
01959 ILL_IFTRACE ("%s:PINF:%lg\n", __func__, ldbl_EGlpNumToLf (lp->tol->ip_tol));
01960 }
01961 else if (fi.dstatus == DUAL_FEASIBLE)
01962 {
01963
01964 it->solstatus = ILL_LP_SOLVED;
01965 ldbl_EGlpNumCopy (lp->objval, lp->pobjval);
01966 it->nextstep = SIMPLEX_TERMINATE;
01967 }
01968 else
01969 ILL_IFTRACE ("%s:DINF:%la:%lf\n", __func__, ldbl_EGlpNumToLf (lp->dinfeas),
01970 ldbl_EGlpNumToLf (lp->dinfeas));
01971 ILL_CLEANUP;
01972 }
01973
01974 ldbl_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[pr.eindex]);
01975 ldbl_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, ldbl_zeroLpNum);
01976 ldbl_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, ldbl_zeroLpNum);
01977 ratio_iter = 0;
01978 do
01979 {
01980 ldbl_ILLratio_pII_test (lp, pr.eindex, pr.dir, &rs);
01981
01982 ldbl_EGlpNumCopy (lbound, rs.lbound);
01983 boundch = rs.boundch;
01984 ratio_iter++;
01985
01986 if (boundch)
01987 {
01988
01989
01990
01991
01992
01993
01994 boundch = 0;
01995 bndtype = (rs.lvstat == STAT_UPPER) ? BOUND_UPPER : BOUND_LOWER;
01996 rval = ldbl_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
02004
02005
02006
02007
02008 it->algorithm = DUAL_SIMPLEX;
02009 it->nextstep = SIMPLEX_RESUME;
02010 it->resumeid = SIMPLEX_RESUME_NUMER;
02011
02012 it->n_restart++;
02013 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02014 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02015
02016 ILL_CLEANUP;
02017 }
02018 else if (rs.ratio_stat == RATIO_UNBOUNDED)
02019 {
02020
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 ldbl_ILLfct_unroll_bound_change (lp);
02029 }
02030 ldbl_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
02038 ldbl_EGlpNumAddInnProdTo (lp->pobjval, rs.tz, lp->dz[pr.eindex]);
02039 ldbl_EGlpNumCopy (lp->objval, lp->pobjval);
02040 if (!ldbl_test_progress (lp->pobjval, it->prevobj))
02041 it->noprog++;
02042 else
02043 {
02044 ldbl_EGlpNumCopy (it->prevobj, lp->pobjval);
02045 it->noprog = 0;
02046 }
02047
02048
02049 ldbl_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
02050 ldbl_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
02051 if (pinf->p_strategy == COMPLETE_PRICING)
02052 ldbl_ILLprice_compute_dual_inf (lp, pinf, &pr.eindex, 1, PRIMAL_PHASEII);
02053 else if (pinf->p_strategy == MULTI_PART_PRICING)
02054 ldbl_ILLprice_update_mpartial_price (lp, pinf, cphase, COL_PRICING);
02055 }
02056 else if (rs.ratio_stat == RATIO_BCHANGE)
02057 {
02058 ldbl_EGlpNumCopyFrac (alpha, lp->dz[pr.eindex], rs.pivotval);
02059 ldbl_EGlpNumAddInnProdTo (lp->pobjval, rs.tz, lp->dz[pr.eindex]);
02060 ldbl_EGlpNumCopy (lp->objval, lp->pobjval);
02061
02062 if (!ldbl_test_progress (lp->pobjval, it->prevobj))
02063 {
02064
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 ldbl_EGlpNumCopy (it->prevobj, lp->pobjval);
02077 it->noprog = 0;
02078 }
02079
02080
02081 ldbl_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
02082 ldbl_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, ldbl_zeroLpNum);
02083 if (pinf->p_strategy == COMPLETE_PRICING)
02084 {
02085
02086 ldbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02087 ldbl_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, ldbl_zeroLpNum);
02088 if (pinf->pII_price == QS_PRICE_PSTEEP)
02089 ldbl_ILLfct_compute_psteep_upv (lp, wz);
02090 }
02091 rval =
02092 ldbl_ILLprice_update_pricing_info (lp, pinf, cphase, wz, pr.eindex, rs.lindex,
02093 rs.pivotval);
02094 CHECKRVALG (rval, CLEANUP);
02095
02096
02097 ldbl_ILLfct_update_xz (lp, rs.tz, pr.eindex, rs.lindex);
02098 ldbl_ILLfct_update_basis_info (lp, pr.eindex, rs.lindex, rs.lvstat);
02099 rval = ldbl_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
02107 it->n_restart++;
02108
02109 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02110 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02111 ILL_CLEANUP;
02112 }
02113 if (!refactor)
02114 {
02115 ldbl_ILLfct_update_piz (lp, alpha);
02116
02117 if (pinf->p_strategy == COMPLETE_PRICING)
02118 {
02119 ldbl_ILLfct_update_dz (lp, pr.eindex, alpha);
02120 ldbl_ILLprice_compute_dual_inf (lp, pinf, lp->zA.indx, lp->zA.nzcnt,
02121 PRIMAL_PHASEII);
02122 ldbl_ILLfct_update_counts (lp, CNT_ZARAVG, lp->zA.nzcnt, ldbl_zeroLpNum);
02123 }
02124 else if (pinf->p_strategy == MULTI_PART_PRICING)
02125 {
02126 ldbl_ILLprice_update_mpartial_price (lp, pinf, cphase, COL_PRICING);
02127 }
02128 }
02129 if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02130 {
02131 ldbl_ILLfct_compute_xbz (lp);
02132 it->newphase = SIMPLEX_PHASE_RECOMP;
02133 }
02134 }
02135
02136 CLEANUP:
02137 ldbl_EGlpNumClearVar (alpha);
02138 ldbl_EGlpNumClearVar (lbound);
02139 ldbl_EGlpNumClearVar (fi.totinfeas);
02140 ldbl_EGlpNumClearVar (pr.dinfeas);
02141 ldbl_EGlpNumClearVar (pr.pinfeas);
02142 ldbl_EGlpNumClearVar (rs.tz);
02143 ldbl_EGlpNumClearVar (rs.lbound);
02144 ldbl_EGlpNumClearVar (rs.ecoeff);
02145 ldbl_EGlpNumClearVar (rs.pivotval);
02146
02147 return rval;
02148 }
02149
02150 static int ldbl_dual_phaseI_step (
02151 ldbl_lpinfo * lp,
02152 ldbl_price_info * pinf,
02153 ldbl_svector * updz,
02154 ldbl_svector * wz,
02155 ldbl_iter_info * it)
02156 {
02157 int rval = 0;
02158 int singular = 0;
02159 int refactor = 0;
02160 int cphase = DUAL_PHASEI;
02161 long double alpha;
02162 long double alpha1;
02163 ldbl_feas_info fi;
02164 ldbl_ratio_res rs;
02165 ldbl_price_res pr;
02166
02167 ldbl_EGlpNumInitVar (alpha);
02168 ldbl_EGlpNumInitVar (alpha1);
02169 ldbl_EGlpNumInitVar (fi.totinfeas);
02170 ldbl_EGlpNumInitVar (pr.dinfeas);
02171 ldbl_EGlpNumInitVar (pr.pinfeas);
02172 ldbl_EGlpNumInitVar (rs.tz);
02173 ldbl_EGlpNumInitVar (rs.lbound);
02174 ldbl_EGlpNumInitVar (rs.ecoeff);
02175 ldbl_EGlpNumInitVar (rs.pivotval);
02176 ldbl_EGlpNumZero (alpha1);
02177
02178 ldbl_ILLfct_update_counts (lp, CNT_DPHASE1ITER, 0, ldbl_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 ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->dinfeas);
02199
02200 ldbl_ILLfct_compute_phaseI_xbz (lp);
02201 if (pinf->d_strategy == COMPLETE_PRICING)
02202 {
02203 #if USEHEAP > 0
02204 ldbl_ILLprice_free_heap (pinf);
02205 #endif
02206 ldbl_ILLprice_compute_primal_inf (lp, pinf, NULL, 0, DUAL_PHASEI);
02207 #if USEHEAP > 0
02208 rval = ldbl_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 ldbl_ILLprice_init_mpartial_price (lp, pinf, cphase, ROW_PRICING);
02216 }
02217 }
02218
02219 ldbl_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 ldbl_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 = ldbl_ILLsimplex_retest_dsolution (lp, pinf, cphase, &fi);
02236 CHECKRVALG (rval, CLEANUP);
02237 ldbl_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 ldbl_ILLfct_compute_zz (lp, &(lp->zz), pr.lindex);
02253
02254 ldbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02255 ldbl_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, ldbl_zeroLpNum);
02256 ldbl_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, ldbl_zeroLpNum);
02257
02258 ldbl_ILLratio_dI_test (lp, pr.lindex, pr.lvstat, &rs);
02259
02260 if (rs.ratio_stat == RATIO_FAILED)
02261 {
02262
02263
02264
02265
02266 it->algorithm = PRIMAL_SIMPLEX;
02267 it->nextstep = SIMPLEX_RESUME;
02268 it->resumeid = SIMPLEX_RESUME_NUMER;
02269
02270 it->n_restart++;
02271 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02272 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02273
02274 ILL_CLEANUP;
02275 }
02276 else if (rs.ratio_stat == RATIO_BCHANGE)
02277 {
02278 ldbl_ILLfct_compute_yz (lp, &(lp->yjz), updz, lp->nbaz[rs.eindex]);
02279 rval = ldbl_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
02287 it->n_restart++;
02288 it->algorithm = PRIMAL_SIMPLEX;
02289 it->nextstep = SIMPLEX_RESUME;
02290 it->resumeid = SIMPLEX_RESUME_NUMER;
02291 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02292 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02293
02294 rval = 0;
02295 ILL_CLEANUP;
02296 }
02297 rval = ldbl_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 ldbl_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, ldbl_zeroLpNum);
02306 ldbl_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, ldbl_zeroLpNum);
02307
02308 if (pinf->dI_price == QS_PRICE_DSTEEP)
02309 ldbl_ILLfct_compute_dsteep_upv (lp, wz);
02310 rval =
02311 ldbl_ILLprice_update_pricing_info (lp, pinf, cphase, wz, rs.eindex, pr.lindex,
02312 rs.pivotval);
02313 CHECKRVALG (rval, CLEANUP);
02314
02315 ldbl_EGlpNumSubTo (lp->dinfeas, lp->upd.c_obj);
02316
02317 if (!ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->dinfeas);
02331 it->noprog = 0;
02332 }
02333
02334 ldbl_EGlpNumCopyFrac (alpha, lp->dz[rs.eindex], rs.pivotval);
02335 ldbl_EGlpNumCopyFrac (alpha1, lp->xbz[pr.lindex], rs.pivotval);
02336
02337 ldbl_ILLfct_update_piz (lp, alpha);
02338 ldbl_ILLfct_update_dz (lp, rs.eindex, alpha);
02339 ldbl_ILLfct_update_dfeas (lp, rs.eindex, &(lp->srhs));
02340 ldbl_ILLfct_compute_dpIy (lp, &(lp->srhs), &(lp->ssoln));
02341 ldbl_ILLfct_update_basis_info (lp, rs.eindex, pr.lindex, pr.lvstat);
02342
02343 #if DENSE_PI > 0
02344 ldbl_fct_test_workvector (lp);
02345 ldbl_fct_test_dfeasible (lp);
02346 #endif
02347 rval = ldbl_ILLbasis_update (lp, updz, pr.lindex, &refactor, &singular);
02348 CHECKRVALG (rval, CLEANUP);
02349
02350 #if DENSE_NORM > 0
02351 ldbl_test_dsteep_norms (lp, pinf);
02352 #endif
02353
02354 ldbl_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
02363 it->n_restart++;
02364
02365 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02366 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02367 ILL_CLEANUP;
02368 }
02369 if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02370 {
02371 ldbl_ILLfct_compute_piz (lp);
02372 ldbl_ILLfct_compute_dz (lp);
02373 ldbl_ILLfct_dual_adjust (lp, ldbl_zeroLpNum);
02374 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->id_tol);
02375 ldbl_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 ldbl_fct_test_workvector (lp);
02386 ldbl_fct_test_pI_x (lp, pinf);
02387 #endif
02388
02389 CLEANUP:
02390 ldbl_EGlpNumClearVar (alpha);
02391 ldbl_EGlpNumClearVar (alpha1);
02392 ldbl_EGlpNumClearVar (fi.totinfeas);
02393 ldbl_EGlpNumClearVar (pr.dinfeas);
02394 ldbl_EGlpNumClearVar (pr.pinfeas);
02395 ldbl_EGlpNumClearVar (rs.tz);
02396 ldbl_EGlpNumClearVar (rs.lbound);
02397 ldbl_EGlpNumClearVar (rs.ecoeff);
02398 ldbl_EGlpNumClearVar (rs.pivotval);
02399
02400 return rval;
02401 }
02402
02403 static int ldbl_dual_phaseII_step (
02404 ldbl_lpinfo * lp,
02405 ldbl_price_info * pinf,
02406 ldbl_svector * updz,
02407 ldbl_svector * wz,
02408 ldbl_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 long double x_bi, v_l, eval;
02419 long double ecoeff;
02420 long double alpha;
02421 long double alpha1;
02422 ldbl_feas_info fi;
02423 ldbl_ratio_res rs;
02424 ldbl_price_res pr;
02425
02426 ldbl_EGlpNumInitVar (x_bi);
02427 ldbl_EGlpNumInitVar (v_l);
02428 ldbl_EGlpNumInitVar (eval);
02429 ldbl_EGlpNumInitVar (ecoeff);
02430 ldbl_EGlpNumInitVar (alpha);
02431 ldbl_EGlpNumInitVar (alpha1);
02432 ldbl_EGlpNumInitVar (fi.totinfeas);
02433 ldbl_EGlpNumInitVar (pr.dinfeas);
02434 ldbl_EGlpNumInitVar (pr.pinfeas);
02435 ldbl_EGlpNumInitVar (rs.tz);
02436 ldbl_EGlpNumInitVar (rs.lbound);
02437 ldbl_EGlpNumInitVar (rs.ecoeff);
02438 ldbl_EGlpNumInitVar (rs.pivotval);
02439 ldbl_EGlpNumZero (rs.ecoeff);
02440 ldbl_EGlpNumZero (alpha1);
02441
02442 ldbl_ILLfct_update_counts (lp, CNT_DPHASE2ITER, 0, ldbl_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 ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->dobjval);
02464 ldbl_ILLfct_compute_xbz (lp);
02465
02466 if (pinf->d_strategy == COMPLETE_PRICING)
02467 {
02468 #if USEHEAP > 0
02469 ldbl_ILLprice_free_heap (pinf);
02470 #endif
02471 ldbl_ILLprice_compute_primal_inf (lp, pinf, NULL, 0, DUAL_PHASEII);
02472 #if USEHEAP > 0
02473 rval = ldbl_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 ldbl_ILLprice_init_mpartial_price (lp, pinf, cphase, ROW_PRICING);
02481 }
02482 }
02483
02484 ldbl_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 ldbl_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 ldbl_ILLfct_unroll_coef_change (lp);
02501 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
02502 ldbl_ILLfct_set_status_values (lp, -1, fi.dstatus, -1, PHASEII);
02503
02504 ldbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
02505
02506
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
02515 it->n_restart++;
02516 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02517 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02518
02519 ILL_CLEANUP;
02520
02521
02522
02523
02524
02525
02526 }
02527 }
02528 if (it->sdisplay > 1)
02529 {
02530
02531 printf ("problem seemingly solved\n");
02532 printf ("seemingly dual opt = %f\n", ldbl_EGlpNumToLf (lp->dobjval));
02533 printf ("retesting soln\n");
02534 fflush (stdout);
02535 }
02536
02537 rval = ldbl_ILLsimplex_retest_dsolution (lp, pinf, cphase, &fi);
02538 CHECKRVALG (rval, CLEANUP);
02539 ldbl_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 ldbl_EGlpNumDivUiTo (lp->tol->ip_tol, 5);
02546 ldbl_EGlpNumDivUiTo (lp->tol->id_tol, 5);
02547 }
02548 else if (fi.pstatus == PRIMAL_FEASIBLE)
02549 {
02550 ILL_IFTRACE ("PRIM_FEAS: %s\n", __func__);
02551 ldbl_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 ldbl_ILLfct_compute_zz (lp, &(lp->zz), pr.lindex);
02561
02562 ldbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02563 ldbl_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, ldbl_zeroLpNum);
02564 ldbl_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, ldbl_zeroLpNum);
02565
02566 ratio_iter = 0;
02567 do
02568 {
02569 ldbl_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 ldbl_ILLfct_adjust_viol_coefs (lp);
02578 ldbl_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 ldbl_EGlpNumCopy (ecoeff, rs.ecoeff);
02590 ratio_iter++;
02591
02592 if (coeffch)
02593 {
02594
02595
02596
02597
02598
02599
02600 coeffch = 0;
02601 rval = ldbl_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
02614
02615
02616 it->algorithm = PRIMAL_SIMPLEX;
02617 it->nextstep = SIMPLEX_RESUME;
02618 it->resumeid = SIMPLEX_RESUME_NUMER;
02619
02620 it->n_restart++;
02621 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02622 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02623
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 ldbl_ILLfct_unroll_coef_change (lp);
02637 }
02638 ldbl_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 ldbl_ILLfct_compute_yz (lp, &(lp->yjz), updz, ecol);
02648 ldbl_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, ldbl_zeroLpNum);
02649 ldbl_ILLfct_update_counts (lp, CNT_UPNZ, updz->nzcnt, ldbl_zeroLpNum);
02650 rval = ldbl_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
02658 it->algorithm = PRIMAL_SIMPLEX;
02659 it->nextstep = SIMPLEX_RESUME;
02660 it->resumeid = SIMPLEX_RESUME_NUMER;
02661 it->n_restart++;
02662 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02663 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02664
02665 rval = 0;
02666 ILL_CLEANUP;
02667 }
02668 if (newphase == 0)
02669 {
02670 rval = ldbl_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 ldbl_EGlpNumAddTo (lp->dobjval, lp->upd.c_obj);
02696 ldbl_EGlpNumCopy (lp->objval, lp->dobjval);
02697
02698 if (!ldbl_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 ldbl_EGlpNumCopy (it->prevobj, lp->dobjval);
02711 it->noprog = 0;
02712 }
02713
02714 if (pinf->dII_price == QS_PRICE_DSTEEP)
02715 ldbl_ILLfct_compute_dsteep_upv (lp, wz);
02716 rval =
02717 ldbl_ILLprice_update_pricing_info (lp, pinf, cphase, wz, rs.eindex, pr.lindex,
02718 rs.pivotval);
02719 CHECKRVALG (rval, CLEANUP);
02720
02721 ldbl_EGlpNumCopy (x_bi, lp->xbz[pr.lindex]);
02722 if (pr.lvstat == STAT_LOWER)
02723 ldbl_EGlpNumCopy (v_l, lp->lz[lcol]);
02724 else
02725 ldbl_EGlpNumCopy (v_l, lp->uz[lcol]);
02726 ldbl_EGlpNumCopy (alpha, rs.tz);
02727 if (pr.lvstat == STAT_LOWER)
02728 ldbl_EGlpNumSign (alpha);
02729 estat = lp->vstat[ecol];
02730 if (estat == STAT_LOWER)
02731 ldbl_EGlpNumCopy (eval, lp->lz[ecol]);
02732 else if (estat == STAT_ZERO)
02733 ldbl_EGlpNumZero (eval);
02734 else
02735 ldbl_EGlpNumCopy (eval, lp->uz[ecol]);
02736
02737 ldbl_ILLfct_update_piz (lp, alpha);
02738 ldbl_ILLfct_update_dz (lp, rs.eindex, alpha);
02739 ldbl_ILLfct_update_dIIfeas (lp, rs.eindex, &(lp->srhs));
02740 ldbl_ILLfct_compute_dpIIy (lp, &(lp->srhs), &(lp->ssoln));
02741 ldbl_EGlpNumCopyDiff (alpha1, x_bi, v_l);
02742 ldbl_EGlpNumSubTo (alpha1, lp->upd.dty);
02743 ldbl_EGlpNumDivTo (alpha1, rs.pivotval);
02744 ldbl_ILLfct_update_basis_info (lp, rs.eindex, pr.lindex, pr.lvstat);
02745 rval = ldbl_ILLbasis_update (lp, updz, pr.lindex, &refactor, &singular);
02746 CHECKRVALG (rval, CLEANUP);
02747
02748 ldbl_ILLfct_update_dpII_prices (lp, pinf, &(lp->srhs), &(lp->ssoln),
02749 pr.lindex, eval, alpha1);
02750
02751 #if DENSE_NORM > 0
02752 ldbl_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
02761 it->n_restart++;
02762 ldbl_EGlpNumMultUiTo (lp->tol->pfeas_tol, SIMPLEX_FACTOR);
02763 ldbl_EGlpNumMultUiTo (lp->tol->dfeas_tol, SIMPLEX_FACTOR);
02764
02765 ILL_CLEANUP;
02766 }
02767 if (refactor != 0 || it->nosolve > PARAM_MAX_NOSOLVE)
02768 {
02769 ldbl_ILLfct_compute_piz (lp);
02770 ldbl_ILLfct_compute_dz (lp);
02771 ldbl_ILLfct_dual_adjust (lp, ldbl_zeroLpNum);
02772 it->newphase = SIMPLEX_PHASE_RECOMP;
02773 }
02774 }
02775
02776 #if DENSE_PIIPI > 0
02777 ldbl_fct_test_workvector (lp);
02778 if (!refactor)
02779 {
02780 ldbl_fct_test_pII_x (lp, pinf);
02781 ldbl_fct_test_pII_pi_dz (lp, pinf);
02782 }
02783 #endif
02784
02785 CLEANUP:
02786 ldbl_EGlpNumClearVar (x_bi);
02787 ldbl_EGlpNumClearVar (v_l);
02788 ldbl_EGlpNumClearVar (eval);
02789 ldbl_EGlpNumClearVar (ecoeff);
02790 ldbl_EGlpNumClearVar (alpha);
02791 ldbl_EGlpNumClearVar (alpha1);
02792 ldbl_EGlpNumClearVar (fi.totinfeas);
02793 ldbl_EGlpNumClearVar (pr.dinfeas);
02794 ldbl_EGlpNumClearVar (pr.pinfeas);
02795 ldbl_EGlpNumClearVar (rs.tz);
02796 ldbl_EGlpNumClearVar (rs.lbound);
02797 ldbl_EGlpNumClearVar (rs.ecoeff);
02798 ldbl_EGlpNumClearVar (rs.pivotval);
02799
02800 return rval;
02801 }
02802
02803 static void ldbl_get_current_stat (
02804 ldbl_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 ldbl_ILLsimplex_pivotin (
02835 ldbl_lpinfo * lp,
02836 ldbl_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 ldbl_svector wz;
02850 ldbl_svector updz;
02851 long double alpha;
02852 ldbl_ratio_res rs;
02853 ldbl_feas_info fi;
02854
02855 ldbl_EGlpNumInitVar (alpha);
02856 ldbl_EGlpNumInitVar (fi.totinfeas);
02857 ldbl_EGlpNumInitVar (rs.tz);
02858 ldbl_EGlpNumInitVar (rs.lbound);
02859 ldbl_EGlpNumInitVar (rs.ecoeff);
02860 ldbl_EGlpNumInitVar (rs.pivotval);
02861 ldbl_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
02897 ldbl_ILLsvector_init (&wz);
02898 rval = ldbl_ILLsvector_alloc (&wz, lp->nrows);
02899 CHECKRVALG (rval, CLEANUP);
02900 ldbl_ILLsvector_init (&updz);
02901 rval = ldbl_ILLsvector_alloc (&updz, lp->nrows);
02902 CHECKRVALG (rval, CLEANUP);
02903
02904 ldbl_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 ldbl_ILLfct_compute_yz (lp, &(lp->yjz), &updz, lp->nbaz[eindex]);
02913 ldbl_ILLfct_update_counts (lp, CNT_YNZ, lp->yjz.nzcnt, ldbl_zeroLpNum);
02914 ldbl_ILLfct_update_counts (lp, CNT_UPNZ, updz.nzcnt, ldbl_zeroLpNum);
02915
02916 ldbl_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 ldbl_EGlpNumCopyDiff (alpha, lp->lz[lp->baz[rs.lindex]], lp->xbz[rs.lindex]);
02929 ldbl_EGlpNumAddInnProdTo (lp->dobjval, rs.tz, alpha);
02930 }
02931 else
02932 {
02933 ldbl_EGlpNumCopyDiff (alpha, lp->xbz[rs.lindex], lp->uz[lp->baz[rs.lindex]]);
02934 ldbl_EGlpNumAddInnProdTo (lp->dobjval, rs.tz, alpha);
02935 }
02936 ldbl_EGlpNumCopyFrac (alpha, lp->dz[eindex], rs.pivotval);
02937 ldbl_EGlpNumCopy (lp->objval, lp->dobjval);
02938
02939 ldbl_ILLfct_compute_zz (lp, &(lp->zz), rs.lindex);
02940
02941 ldbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
02942 ldbl_ILLfct_update_counts (lp, CNT_ZNZ, lp->zz.nzcnt, ldbl_zeroLpNum);
02943 ldbl_ILLfct_update_counts (lp, CNT_ZANZ, lp->zA.nzcnt, ldbl_zeroLpNum);
02944
02945 if (pinf->dsinfo.norms && pinf->dII_price == QS_PRICE_DSTEEP)
02946 {
02947 ldbl_ILLfct_compute_dsteep_upv (lp, &wz);
02948 rval = ldbl_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 ldbl_ILLfct_compute_psteep_upv (lp, &wz);
02955 rval = ldbl_ILLprice_update_pricing_info (lp, pinf, PRIMAL_PHASEII, &wz,
02956 eindex, rs.lindex, rs.pivotval);
02957 CHECKRVALG (rval, CLEANUP);
02958 }
02959
02960
02961 ldbl_ILLfct_update_xz (lp, rs.tz, eindex, rs.lindex);
02962 ldbl_ILLfct_update_basis_info (lp, eindex, rs.lindex, rs.lvstat);
02963 rval = ldbl_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 ldbl_ILLfct_update_piz (lp, alpha);
02975 ldbl_ILLfct_update_dz (lp, eindex, alpha);
02976 }
02977 else
02978 {
02979 ldbl_ILLfct_compute_xbz (lp);
02980 ldbl_ILLfct_compute_piz (lp);
02981 ldbl_ILLfct_compute_dz (lp);
02982 ldbl_ILLfct_compute_dobj (lp);
02983 }
02984 }
02985 }
02986
02987
02988
02989
02990
02991
02992 ldbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
02993 ldbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
02994 ldbl_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 ldbl_ILLsvector_free (&wz);
03001 ldbl_ILLsvector_free (&updz);
03002 ldbl_EGlpNumClearVar (alpha);
03003 ldbl_EGlpNumClearVar (fi.totinfeas);
03004 ldbl_EGlpNumClearVar (rs.tz);
03005 ldbl_EGlpNumClearVar (rs.lbound);
03006 ldbl_EGlpNumClearVar (rs.ecoeff);
03007 ldbl_EGlpNumClearVar (rs.pivotval);
03008 EG_RETURN (rval);
03009 }
03010
03011 static int ldbl_report_value (
03012 ldbl_lpinfo * lp,
03013 ldbl_iter_info * it,
03014 const char *value_name,
03015 long double 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, ldbl_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
03033 if (it->itercnt % (lp->iterskip / 10))
03034 {
03035 rval = ILLstring_report (NULL, &lp->O->reporter);
03036 }
03037 }
03038 if (rval != 0)
03039 {
03040 it->solstatus = QS_LP_ABORTED;
03041 }
03042 return rval;
03043 }
03044
03045 #undef ldbl_QSOPT_CURRENT_PRECICION
