#include <string.h>
#include "commonlib.h"
#include "lp_lib.h"
#include "lp_BFP.h"
#include "lp_simplex.h"
#include "lp_crash.h"
#include "lp_presolve.h"
#include "lp_price.h"
#include "lp_pricePSE.h"
#include "lp_report.h"
#ifdef FORTIFY
# include "lp_fortify.h"
#endif
STATIC void stallMonitor_update(lprec *lp, REAL newOF)
{
int newpos;
OBJmonrec *monitor = lp->monitor;
if(monitor->countstep < OBJ_STEPS)
monitor->countstep++;
else
monitor->startstep = mod(monitor->startstep + 1, OBJ_STEPS);
newpos = mod(monitor->startstep + monitor->countstep - 1, OBJ_STEPS);
monitor->objstep[newpos] = newOF;
monitor->idxstep[newpos] = monitor->Icount;
monitor->currentstep = newpos;
}
STATIC MYBOOL stallMonitor_creepingObj(lprec *lp)
{
OBJmonrec *monitor = lp->monitor;
if(monitor->countstep > 1) {
REAL deltaOF = (monitor->objstep[monitor->currentstep] -
monitor->objstep[monitor->startstep]) / monitor->countstep;
deltaOF /= MAX(1, (monitor->idxstep[monitor->currentstep] -
monitor->idxstep[monitor->startstep]));
deltaOF = my_chsign(monitor->isdual, deltaOF);
return( (MYBOOL) (deltaOF < monitor->epsvalue) );
}
else
return( FALSE );
}
STATIC MYBOOL stallMonitor_shortSteps(lprec *lp)
{
OBJmonrec *monitor = lp->monitor;
if(monitor->countstep == OBJ_STEPS) {
REAL deltaOF = MAX(1, (monitor->idxstep[monitor->currentstep] -
monitor->idxstep[monitor->startstep])) / monitor->countstep;
deltaOF = pow(deltaOF*OBJ_STEPS, 0.66);
return( (MYBOOL) (deltaOF > monitor->limitstall[TRUE]) );
}
else
return( FALSE );
}
STATIC void stallMonitor_reset(lprec *lp)
{
OBJmonrec *monitor = lp->monitor;
monitor->ruleswitches = 0;
monitor->Ncycle = 0;
monitor->Mcycle = 0;
monitor->Icount = 0;
monitor->startstep = 0;
monitor->objstep[monitor->startstep] = lp->infinite;
monitor->idxstep[monitor->startstep] = monitor->Icount;
monitor->prevobj = 0;
monitor->countstep = 1;
}
STATIC MYBOOL stallMonitor_create(lprec *lp, MYBOOL isdual, char *funcname)
{
OBJmonrec *monitor = NULL;
if(lp->monitor != NULL)
return( FALSE );
monitor = (OBJmonrec *) calloc(sizeof(*monitor), 1);
if(monitor == NULL)
return( FALSE );
monitor->lp = lp;
strcpy(monitor->spxfunc, funcname);
monitor->isdual = isdual;
monitor->pivdynamic = is_piv_mode(lp, PRICE_ADAPTIVE);
monitor->oldpivstrategy = lp->piv_strategy;
monitor->oldpivrule = get_piv_rule(lp);
if(MAX_STALLCOUNT <= 1)
monitor->limitstall[FALSE] = 0;
else
monitor->limitstall[FALSE] = MAX(MAX_STALLCOUNT,
(int) pow((REAL) (lp->rows+lp->columns)/2, 0.667));
#if 1
monitor->limitstall[FALSE] *= 2+2;
#endif
monitor->limitstall[TRUE] = monitor->limitstall[FALSE];
if(monitor->oldpivrule == PRICER_DEVEX)
monitor->limitstall[TRUE] *= 2;
if(MAX_RULESWITCH <= 0)
monitor->limitruleswitches = MAXINT32;
else
monitor->limitruleswitches = MAX(MAX_RULESWITCH,
lp->rows/MAX_RULESWITCH);
monitor->epsvalue = lp->epsprimal;
lp->monitor = monitor;
stallMonitor_reset(lp);
lp->suminfeas = lp->infinite;
return( TRUE );
}
STATIC MYBOOL stallMonitor_check(lprec *lp, int rownr, int colnr, int lastnr,
MYBOOL minit, MYBOOL approved, MYBOOL *forceoutEQ)
{
OBJmonrec *monitor = lp->monitor;
MYBOOL isStalled, isCreeping, acceptance = TRUE;
int altrule,
#ifdef Paranoia
msglevel = NORMAL;
#else
msglevel = DETAILED;
#endif
REAL deltaobj = lp->suminfeas;
monitor->active = FALSE;
if(monitor->Icount <= 1) {
if(monitor->Icount == 1) {
monitor->prevobj = lp->rhs[0];
monitor->previnfeas = deltaobj;
}
monitor->Icount++;
return( acceptance );
}
monitor->thisobj = lp->rhs[0];
monitor->thisinfeas = deltaobj;
if(lp->spx_trace &&
(lastnr > 0))
report(lp, NORMAL, "%s: Objective at iter %10.0f is " RESULTVALUEMASK " (%4d: %4d %s- %4d)\n",
monitor->spxfunc,
(double) get_total_iter(lp), monitor->thisobj, rownr, lastnr,
my_if(minit == ITERATE_MAJORMAJOR, "<","|"), colnr);
monitor->pivrule = get_piv_rule(lp);
deltaobj = my_reldiff(monitor->thisobj, monitor->prevobj);
deltaobj = fabs(deltaobj);
isStalled = (MYBOOL) (deltaobj < monitor->epsvalue);
if(isStalled) {
REAL testvalue, refvalue = monitor->epsvalue;
#if 1
if(monitor->isdual)
refvalue *= 1000*log10(9.0+lp->rows);
else
refvalue *= 1000*log10(9.0+lp->columns);
#else#endif
testvalue = my_reldiff(monitor->thisinfeas, monitor->previnfeas);
isStalled &= (fabs(testvalue) < refvalue);
#if !defined _PRICE_NOBOUNDFLIP
if(!isStalled && (testvalue > 0) && is_action(lp->anti_degen, ANTIDEGEN_BOUNDFLIP))
acceptance = AUTOMATIC;
}
#else
if(!isStalled && (testvalue > 0) && !ISMASKSET(lp->piv_strategy, PRICE_NOBOUNDFLIP)) {
SETMASK(lp->piv_strategy, PRICE_NOBOUNDFLIP);
acceptance = AUTOMATIC;
}
}
else
CLEARMASK(lp->piv_strategy, PRICE_NOBOUNDFLIP);
#endif
#if 1
isCreeping = FALSE;
#else#endif
if(isStalled || isCreeping) {
#if 1
if(minit != ITERATE_MAJORMAJOR) {
if(++monitor->Mcycle > 2) {
monitor->Mcycle = 0;
monitor->Ncycle++;
}
}
else
#endif
monitor->Ncycle++;
if(monitor->Ncycle <= 1) {
monitor->Ccycle = colnr;
monitor->Rcycle = rownr;
}
else if(isCreeping ||
(monitor->Ncycle > monitor->limitstall[monitor->isdual]) ||
((monitor->Ccycle == rownr) && (monitor->Rcycle == colnr))) {
monitor->active = TRUE;
if((lp->fixedvars > 0) && (*forceoutEQ != TRUE)) {
*forceoutEQ = TRUE;
goto Proceed;
}
approved &= monitor->pivdynamic && (monitor->ruleswitches < monitor->limitruleswitches);
if(!approved && !is_anti_degen(lp, ANTIDEGEN_STALLING)) {
lp->spx_status = DEGENERATE;
report(lp, msglevel, "%s: Stalling at iter %10.0f; no alternative strategy left.\n",
monitor->spxfunc, (double) get_total_iter(lp));
acceptance = FALSE;
return( acceptance );
}
switch (monitor->oldpivrule) {
case PRICER_FIRSTINDEX: altrule = PRICER_DEVEX;
break;
case PRICER_DANTZIG: altrule = PRICER_DEVEX;
break;
case PRICER_DEVEX: altrule = PRICER_STEEPESTEDGE;
break;
case PRICER_STEEPESTEDGE: altrule = PRICER_DEVEX;
break;
default: altrule = PRICER_FIRSTINDEX;
}
if(approved &&
(monitor->pivrule != altrule) && (monitor->pivrule == monitor->oldpivrule)) {
monitor->ruleswitches++;
lp->piv_strategy = altrule;
monitor->Ccycle = 0;
monitor->Rcycle = 0;
monitor->Ncycle = 0;
monitor->Mcycle = 0;
report(lp, msglevel, "%s: Stalling at iter %10.0f; changed to '%s' rule.\n",
monitor->spxfunc, (double) get_total_iter(lp),
get_str_piv_rule(get_piv_rule(lp)));
if((altrule == PRICER_DEVEX) || (altrule == PRICER_STEEPESTEDGE))
restartPricer(lp, AUTOMATIC);
}
else {
report(lp, msglevel, "%s: Stalling at iter %10.0f; proceed to bound relaxation.\n",
monitor->spxfunc, (double) get_total_iter(lp));
acceptance = FALSE;
lp->spx_status = DEGENERATE;
return( acceptance );
}
}
}
else {
if(monitor->pivrule != monitor->oldpivrule) {
lp->piv_strategy = monitor->oldpivstrategy;
altrule = monitor->oldpivrule;
if((altrule == PRICER_DEVEX) || (altrule == PRICER_STEEPESTEDGE))
restartPricer(lp, AUTOMATIC);
report(lp, msglevel, "...returned to original pivot selection rule at iter %.0f.\n",
(double) get_total_iter(lp));
}
stallMonitor_update(lp, monitor->thisobj);
monitor->Ccycle = 0;
monitor->Rcycle = 0;
monitor->Ncycle = 0;
monitor->Mcycle = 0;
}
Proceed:
monitor->Icount++;
if(deltaobj >= monitor->epsvalue)
monitor->prevobj = monitor->thisobj;
monitor->previnfeas = monitor->thisinfeas;
return( acceptance );
}
STATIC void stallMonitor_finish(lprec *lp)
{
OBJmonrec *monitor = lp->monitor;
if(monitor == NULL)
return;
if(lp->piv_strategy != monitor->oldpivstrategy)
lp->piv_strategy = monitor->oldpivstrategy;
FREE(monitor);
lp->monitor = NULL;
}
STATIC MYBOOL add_artificial(lprec *lp, int forrownr, REAL *nzarray, int *idxarray)
{
MYBOOL add;
add = !isBasisVarFeasible(lp, lp->epspivot, forrownr);
if(add) {
int *rownr = NULL, i, bvar, ii;
REAL *avalue = NULL, rhscoef, acoef;
MATrec *mat = lp->matA;
for(i = 1; i <= lp->rows; i++) {
if(lp->var_basic[i] == forrownr)
break;
}
acoef = 1;
if(i > lp->rows) {
for(i = 1; i <= lp->rows; i++) {
ii = lp->var_basic[i] - lp->rows;
if((ii <= 0) || (ii > (lp->columns-lp->P1extraDim)))
continue;
ii = mat_findelm(mat, forrownr, ii);
if(ii >= 0) {
acoef = COL_MAT_VALUE(ii);
break;
}
}
}
#if 0#endif
bvar = i;
add = (MYBOOL) (bvar <= lp->rows);
if(add) {
rhscoef = lp->rhs[forrownr];
if(nzarray == NULL)
allocREAL(lp, &avalue, 2, FALSE);
else
avalue = nzarray;
if(idxarray == NULL)
allocINT(lp, &rownr, 2, FALSE);
else
rownr = idxarray;
rownr[0] = 0;
avalue[0] = my_chsign(is_chsign(lp, 0), 1);
rownr[1] = forrownr;
avalue[1] = my_chsign(is_chsign(lp, forrownr), my_sign(rhscoef/acoef));
add_columnex(lp, 2, avalue, rownr);
if(idxarray == NULL)
FREE(rownr);
if(nzarray == NULL)
FREE(avalue);
set_basisvar(lp, bvar, lp->sum);
lp->P1extraDim++;
}
else {
report(lp, CRITICAL, "add_artificial: Could not find replacement basis variable for row %d\n",
forrownr);
lp->basis_valid = FALSE;
}
}
return(add);
}
STATIC int get_artificialRow(lprec *lp, int colnr)
{
MATrec *mat = lp->matA;
#ifdef Paranoia
if((colnr <= lp->columns-abs(lp->P1extraDim)) || (colnr > lp->columns))
report(lp, SEVERE, "get_artificialRow: Invalid column index %d\n", colnr);
if(mat->col_end[colnr] - mat->col_end[colnr-1] != 1)
report(lp, SEVERE, "get_artificialRow: Invalid column non-zero count\n");
#endif
colnr = mat->col_end[colnr-1];
colnr = COL_MAT_ROWNR(colnr);
return( colnr );
}
STATIC int findAnti_artificial(lprec *lp, int colnr)
{
int i, k, rownr = 0, P1extraDim = abs(lp->P1extraDim);
if((P1extraDim == 0) || (colnr > lp->rows) || !lp->is_basic[colnr])
return( rownr );
for(i = 1; i <= lp->rows; i++) {
k = lp->var_basic[i];
if((k > lp->sum-P1extraDim) && (lp->rhs[i] == 0)) {
rownr = get_artificialRow(lp, k-lp->rows);
if(rownr == colnr)
break;
rownr = 0;
}
}
return( rownr );
}
STATIC int findBasicArtificial(lprec *lp, int before)
{
int i = 0, P1extraDim = abs(lp->P1extraDim);
if(P1extraDim > 0) {
if(before > lp->rows || before <= 1)
i = lp->rows;
else
i = before;
while((i > 0) && (lp->var_basic[i] <= lp->sum-P1extraDim))
i--;
}
return(i);
}
STATIC void eliminate_artificials(lprec *lp, REAL *prow)
{
int i, j, colnr, rownr, P1extraDim = abs(lp->P1extraDim);
for(i = 1; (i <= lp->rows) && (P1extraDim > 0); i++) {
j = lp->var_basic[i];
if(j <= lp->sum-P1extraDim)
continue;
j -= lp->rows;
rownr = get_artificialRow(lp, j);
colnr = find_rowReplacement(lp, rownr, prow, NULL);
#if 0#else
set_basisvar(lp, rownr, colnr);
#endif
del_column(lp, j);
P1extraDim--;
}
lp->P1extraDim = 0;
}
STATIC void clear_artificials(lprec *lp)
{
int i, j, n, P1extraDim;
n = 0;
P1extraDim = abs(lp->P1extraDim);
for(i = 1; (i <= lp->rows) && (n < P1extraDim); i++) {
j = lp->var_basic[i];
if(j <= lp->sum-P1extraDim)
continue;
j = get_artificialRow(lp, j-lp->rows);
set_basisvar(lp, i, j);
n++;
}
#ifdef Paranoia
if(n != lp->P1extraDim)
report(lp, SEVERE, "clear_artificials: Unable to clear all basic artificial variables\n");
#endif
while(P1extraDim > 0) {
i = lp->sum-lp->rows;
del_column(lp, i);
P1extraDim--;
}
lp->P1extraDim = 0;
if(n > 0) {
set_action(&lp->spx_action, ACTION_REINVERT);
lp->basis_valid = TRUE;
}
}
STATIC int primloop(lprec *lp, MYBOOL primalfeasible, REAL primaloffset)
{
MYBOOL primal = TRUE, bfpfinal = FALSE, changedphase = FALSE, forceoutEQ = AUTOMATIC,
primalphase1, pricerCanChange, minit, stallaccept, pendingunbounded;
int i, j, k, colnr = 0, rownr = 0, lastnr = 0,
candidatecount = 0, minitcount = 0, ok = TRUE;
LREAL theta = 0.0;
REAL epsvalue, xviolated = 0.0, cviolated = 0.0,
*prow = NULL, *pcol = NULL,
*drow = lp->drow;
int *workINT = NULL,
*nzdrow = lp->nzdrow;
if(lp->spx_trace)
report(lp, DETAILED, "Entered primal simplex algorithm with feasibility %s\n",
my_boolstr(primalfeasible));
lp->P1extraDim = 0;
if(!primalfeasible) {
lp->simplex_mode = SIMPLEX_Phase1_PRIMAL;
#ifdef Paranoia
if(!verify_basis(lp))
report(lp, SEVERE, "primloop: No valid basis for artificial variables\n");
#endif
#if 0#endif
for(i = 1; i <= lp->rows; i++) {
add_artificial(lp, i, NULL, NULL);
}
if(lp->P1extraDim > 0) {
#if 1
ok = allocREAL(lp, &(lp->drow), lp->sum+1, AUTOMATIC) &&
allocINT(lp, &(lp->nzdrow), lp->sum+1, AUTOMATIC);
if(!ok)
goto Finish;
lp->nzdrow[0] = 0;
drow = lp->drow;
nzdrow = lp->nzdrow;
#endif
mat_validate(lp->matA);
set_OF_p1extra(lp, 0.0);
}
if(lp->spx_trace)
report(lp, DETAILED, "P1extraDim count = %d\n", lp->P1extraDim);
simplexPricer(lp, (MYBOOL)!primal);
invert(lp, INITSOL_USEZERO, TRUE);
}
else {
lp->simplex_mode = SIMPLEX_Phase2_PRIMAL;
restartPricer(lp, (MYBOOL)!primal);
}
ok = allocREAL(lp, &(lp->bsolveVal), lp->rows + 1, FALSE) &&
allocREAL(lp, &prow, lp->sum + 1, TRUE) &&
allocREAL(lp, &pcol, lp->rows + 1, TRUE);
if(is_piv_mode(lp, PRICE_MULTIPLE) && (lp->multiblockdiv > 1)) {
lp->multivars = multi_create(lp, FALSE);
ok &= (lp->multivars != NULL) &&
multi_resize(lp->multivars, lp->sum / lp->multiblockdiv, 2, FALSE, TRUE);
}
if(!ok)
goto Finish;
lp->spx_status = RUNNING;
minit = ITERATE_MAJORMAJOR;
epsvalue = lp->epspivot;
pendingunbounded = FALSE;
ok = stallMonitor_create(lp, FALSE, "primloop");
if(!ok)
goto Finish;
lp->rejectpivot[0] = 0;
while((lp->spx_status == RUNNING) && !userabort(lp, -1)) {
primalphase1 = (MYBOOL) (lp->P1extraDim > 0);
clear_action(&lp->spx_action, ACTION_REINVERT | ACTION_ITERATE);
pricerCanChange = !primalphase1;
stallaccept = stallMonitor_check(lp, rownr, colnr, lastnr, minit, pricerCanChange, &forceoutEQ);
if(!stallaccept)
break;
RetryCol:
#if 0#endif
if(!changedphase) {
i = 0;
do {
i++;
colnr = colprim(lp, drow, nzdrow, (MYBOOL) (minit == ITERATE_MINORRETRY), i, &candidatecount, TRUE, &xviolated);
} while ((colnr == 0) && (i < partial_countBlocks(lp, (MYBOOL) !primal)) &&
partial_blockStep(lp, (MYBOOL) !primal));
if(colnr == 0)
lp->spx_status = OPTIMAL;
if(lp->rejectpivot[0] > 0)
minit = ITERATE_MAJORMAJOR;
if(is_action(lp->spx_action, ACTION_REINVERT))
bfpfinal = TRUE;
}
#ifdef primal_UseRejectionList
if((colnr == 0) && (lp->rejectpivot[0] > 0)) {
lp->spx_status = UNBOUNDED;
if((lp->spx_trace && (lp->bb_totalnodes == 0)) ||
(lp->bb_trace && (lp->bb_totalnodes > 0)))
report(lp, DETAILED, "The model is primal unbounded.\n");
colnr = lp->rejectpivot[1];
rownr = 0;
lp->rejectpivot[0] = 0;
ok = FALSE;
break;
}
#endif
if(colnr > 0) {
changedphase = FALSE;
fsolve(lp, colnr, pcol, NULL, lp->epsmachine, 1.0, TRUE);
if(is_anti_degen(lp, ANTIDEGEN_COLUMNCHECK) && !check_degeneracy(lp, pcol, NULL)) {
if(lp->rejectpivot[0] < DEF_MAXPIVOTRETRY/3) {
i = ++lp->rejectpivot[0];
lp->rejectpivot[i] = colnr;
report(lp, DETAILED, "Entering column %d found to be non-improving due to degeneracy.\n",
colnr);
minit = ITERATE_MINORRETRY;
goto RetryCol;
}
else {
lp->rejectpivot[0] = 0;
report(lp, DETAILED, "Gave up trying to find a strictly improving entering column.\n");
}
}
theta = drow[colnr];
rownr = rowprim(lp, colnr, &theta, pcol, workINT, forceoutEQ, &cviolated);
#ifdef AcceptMarginalAccuracy
if((rownr > 0) && (xviolated+cviolated < lp->epspivot)) {
if(lp->bb_trace || (lp->bb_totalnodes == 0))
report(lp, DETAILED, "primloop: Assuming convergence with reduced accuracy %g.\n",
MAX(xviolated, cviolated));
rownr = 0;
colnr = 0;
goto Optimality;
}
else
#endif
if((lp->P1extraDim != 0) && (rownr == 0) && (colnr <= lp->rows))
rownr = findAnti_artificial(lp, colnr);
if(rownr > 0) {
pendingunbounded = FALSE;
lp->rejectpivot[0] = 0;
set_action(&lp->spx_action, ACTION_ITERATE);
if(!lp->obj_in_basis)
pcol[0] = my_chsign(!lp->is_lower[colnr], drow[colnr]);
lp->bfp_prepareupdate(lp, rownr, colnr, pcol);
}
else {
#ifdef primal_UseRejectionList
if(lp->rejectpivot[0] < DEF_MAXPIVOTRETRY) {
lp->spx_status = RUNNING;
lp->rejectpivot[0]++;
lp->rejectpivot[lp->rejectpivot[0]] = colnr;
report(lp, DETAILED, "...trying to recover via another pivot column.\n");
minit = ITERATE_MINORRETRY;
goto RetryCol;
}
else
#endif
if(!refactRecent(lp) && !pendingunbounded) {
bfpfinal = TRUE;
pendingunbounded = TRUE;
set_action(&lp->spx_action, ACTION_REINVERT);
}
else {
lp->spx_status = UNBOUNDED;
report(lp, DETAILED, "The model is primal unbounded.\n");
break;
}
}
}
else {
Optimality:
if(!primalfeasible || isP1extra(lp)) {
if(feasiblePhase1(lp, epsvalue)) {
lp->spx_status = RUNNING;
if(lp->bb_totalnodes == 0) {
report(lp, NORMAL, "Found feasibility by primal simplex after %10.0f iter.\n",
(double) get_total_iter(lp));
if((lp->usermessage != NULL) && (lp->msgmask & MSG_LPFEASIBLE))
lp->usermessage(lp, lp->msghandle, MSG_LPFEASIBLE);
}
changedphase = FALSE;
primalfeasible = TRUE;
lp->simplex_mode = SIMPLEX_Phase2_PRIMAL;
set_OF_p1extra(lp, 0.0);
if(lp->P1extraDim > 0) {
#ifdef Phase1EliminateRedundant
if((lp->bb_totalnodes == 0) && (MIP_count(lp) == 0)) {
while(lp->P1extraDim > 0) {
i = lp->rows;
while((i > 0) && (lp->var_basic[i] <= lp->sum-lp->P1extraDim))
i--;
#ifdef Paranoia
if(i <= 0) {
report(lp, SEVERE, "primloop: Could not find redundant artificial.\n");
break;
}
#endif
j = lp->var_basic[i]-lp->rows;
k = get_artificialRow(lp, j);
if(lp->is_basic[k]) {
lp->is_basic[lp->rows+j] = FALSE;
del_constraint(lp, k);
}
else
set_basisvar(lp, i, k);
del_column(lp, j);
lp->P1extraDim--;
}
lp->basis_valid = TRUE;
}
else
eliminate_artificials(lp, prow);
#else
lp->P1extraDim = my_flipsign(lp->P1extraDim);
#endif
}
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
else {
lp->spx_status = INFEASIBLE;
minit = ITERATE_MAJORMAJOR;
if(lp->spx_trace)
report(lp, NORMAL, "Model infeasible by primal simplex at iter %10.0f.\n",
(double) get_total_iter(lp));
}
}
else {
}
if((lp->bb_level <= 1) || (lp->improve & IMPROVE_BBSIMPLEX) ) {
set_action(&lp->piv_strategy, PRICE_FORCEFULL);
i = rowdual(lp, lp->rhs, FALSE, FALSE, NULL);
clear_action(&lp->piv_strategy, PRICE_FORCEFULL);
if(i > 0) {
lp->spx_status = LOSTFEAS;
if(lp->total_iter == 0)
report(lp, DETAILED, "primloop: Lost primal feasibility at iter %10.0f: will try to recover.\n",
(double) get_total_iter(lp));
}
}
}
if(is_action(lp->spx_action, ACTION_ITERATE)) {
lastnr = lp->var_basic[rownr];
if(refactRecent(lp) == AUTOMATIC)
minitcount = 0;
else if(minitcount > MAX_MINITUPDATES) {
recompute_solution(lp, INITSOL_USEZERO);
minitcount = 0;
}
minit = performiteration(lp, rownr, colnr, theta, primal,
(MYBOOL) (
(stallaccept != AUTOMATIC)),
NULL, NULL,
pcol, NULL, NULL);
if(minit != ITERATE_MAJORMAJOR)
minitcount++;
if((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT))
break;
else if(minit == ITERATE_MINORMAJOR)
continue;
#ifdef UsePrimalReducedCostUpdate
if(minit == ITERATE_MAJORMAJOR)
update_reducedcosts(lp, primal, lastnr, colnr, pcol, drow);
#endif
if((minit == ITERATE_MAJORMAJOR) && (lastnr > lp->sum - abs(lp->P1extraDim))) {
#ifdef Paranoia
if(lp->is_basic[lastnr] || !lp->is_basic[colnr])
report(lp, SEVERE, "primloop: Invalid basis indicator for variable %d at iter %10.0f.\n",
lastnr, (double) get_total_iter(lp));
#endif
del_column(lp, lastnr-lp->rows);
if(lp->P1extraDim > 0)
lp->P1extraDim--;
else
lp->P1extraDim++;
if(lp->P1extraDim == 0) {
colnr = 0;
changedphase = TRUE;
stallMonitor_reset(lp);
}
}
}
if(lp->spx_status == SWITCH_TO_DUAL)
;
else if(!changedphase && lp->bfp_mustrefactorize(lp)) {
#ifdef ResetMinitOnReinvert
minit = ITERATE_MAJORMAJOR;
#endif
if(!invert(lp, INITSOL_USEZERO, bfpfinal))
lp->spx_status = SINGULAR_BASIS;
bfpfinal = FALSE;
}
}
lp->P1extraDim = abs(lp->P1extraDim);
if(lp->P1extraDim > 0) {
clear_artificials(lp);
if(lp->spx_status != OPTIMAL)
restore_basis(lp);
i = invert(lp, INITSOL_USEZERO, TRUE);
}
#ifdef Paranoia
if(!verify_basis(lp))
report(lp, SEVERE, "primloop: Invalid basis detected due to internal error\n");
#endif
#ifdef ForceDualSimplexInBB
if((lp->bb_totalnodes == 0) && (MIP_count(lp) > 0) &&
((lp->simplex_strategy & SIMPLEX_Phase1_DUAL) == 0)) {
lp->simplex_strategy &= ~SIMPLEX_Phase1_PRIMAL;
lp->simplex_strategy += SIMPLEX_Phase1_DUAL;
}
#endif
Finish:
stallMonitor_finish(lp);
multi_free(&(lp->multivars));
FREE(prow);
FREE(pcol);
FREE(lp->bsolveVal);
return(ok);
}
STATIC int dualloop(lprec *lp, MYBOOL dualfeasible, int dualinfeasibles[], REAL dualoffset)
{
MYBOOL primal = FALSE, inP1extra, dualphase1 = FALSE, changedphase = TRUE,
pricerCanChange, minit, stallaccept, longsteps,
forceoutEQ = FALSE, bfpfinal = FALSE;
int i, colnr = 0, rownr = 0, lastnr = 0,
candidatecount = 0, minitcount = 0,
#ifdef FixInaccurateDualMinit
minitcolnr = 0,
#endif
ok = TRUE;
int *boundswaps = NULL;
LREAL theta = 0.0;
REAL epsvalue, xviolated, cviolated,
*prow = NULL, *pcol = NULL,
*drow = lp->drow;
int *nzprow = NULL, *workINT = NULL,
*nzdrow = lp->nzdrow;
if(lp->spx_trace)
report(lp, DETAILED, "Entered dual simplex algorithm with feasibility %s.\n",
my_boolstr(dualfeasible));
ok = allocREAL(lp, &prow, lp->sum + 1, TRUE) &&
allocINT (lp, &nzprow, lp->sum + 1, FALSE) &&
allocREAL(lp, &pcol, lp->rows + 1, TRUE);
if(!ok)
goto Finish;
inP1extra = (MYBOOL) (dualoffset != 0);
if(inP1extra) {
set_OF_p1extra(lp, dualoffset);
simplexPricer(lp, (MYBOOL)!primal);
invert(lp, INITSOL_USEZERO, TRUE);
}
else
restartPricer(lp, (MYBOOL)!primal);
#if 0#else
longsteps = FALSE;
#endif
#ifdef UseLongStepDualPhase1
longsteps = !dualfeasible && (MYBOOL) (dualinfeasibles != NULL);
#endif
if(longsteps) {
lp->longsteps = multi_create(lp, TRUE);
ok = (lp->longsteps != NULL) &&
multi_resize(lp->longsteps, MIN(lp->boundedvars+2, 11), 1, TRUE, TRUE);
if(!ok)
goto Finish;
#ifdef UseLongStepPruning
lp->longsteps->objcheck = TRUE;
#endif
boundswaps = multi_indexSet(lp->longsteps, FALSE);
}
lp->spx_status = RUNNING;
minit = ITERATE_MAJORMAJOR;
epsvalue = lp->epspivot;
ok = stallMonitor_create(lp, TRUE, "dualloop");
if(!ok)
goto Finish;
lp->rejectpivot[0] = 0;
if(dualfeasible)
lp->simplex_mode = SIMPLEX_Phase2_DUAL;
else
lp->simplex_mode = SIMPLEX_Phase1_DUAL;
if(dualphase1 || inP1extra ||
((lp->fixedvars > 0) && is_anti_degen(lp, ANTIDEGEN_FIXEDVARS))) {
forceoutEQ = AUTOMATIC;
}
#if 1
if(is_anti_degen(lp, ANTIDEGEN_DYNAMIC) && (bin_count(lp, TRUE)*2 > lp->columns)) {
switch (forceoutEQ) {
case FALSE: forceoutEQ = AUTOMATIC;
break;
}
}
#endif
while((lp->spx_status == RUNNING) && !userabort(lp, -1)) {
pricerCanChange = !dualphase1 && !inP1extra;
stallaccept = stallMonitor_check(lp, rownr, colnr, lastnr, minit, pricerCanChange, &forceoutEQ);
if(!stallaccept)
break;
changedphase = FALSE;
dualphase1 &= (MYBOOL) (lp->simplex_mode == SIMPLEX_Phase1_DUAL);
if(longsteps && dualphase1 && !inP1extra) {
obtain_column(lp, dualinfeasibles[1], pcol, NULL, NULL);
i = 2;
for(i = 2; i <= dualinfeasibles[0]; i++)
mat_multadd(lp->matA, pcol, dualinfeasibles[i], 1.0);
ftran(lp, pcol, NULL, lp->epsmachine);
}
#if (defined dual_Phase1PriceEqualities) || (defined dual_UseRejectionList)
RetryRow:
#endif
if(minit != ITERATE_MINORRETRY) {
i = 0;
do {
i++;
rownr = rowdual(lp, my_if(dualphase1, pcol, NULL), forceoutEQ, TRUE, &xviolated);
} while ((rownr == 0) && (i < partial_countBlocks(lp, (MYBOOL) !primal)) &&
partial_blockStep(lp, (MYBOOL) !primal));
}
#ifdef dual_UseRejectionList
if((rownr == 0) && (lp->rejectpivot[0] > 0)) {
lp->spx_status = INFEASIBLE;
if((lp->spx_trace && (lp->bb_totalnodes == 0)) ||
(lp->bb_trace && (lp->bb_totalnodes > 0)))
report(lp, DETAILED, "The model is primal infeasible.\n");
rownr = lp->rejectpivot[1];
colnr = 0;
lp->rejectpivot[0] = 0;
ok = FALSE;
break;
}
#endif
clear_action(&lp->spx_action, ACTION_ITERATE);
if(rownr > 0) {
colnr = coldual(lp, rownr, prow, nzprow, drow, nzdrow,
(MYBOOL) (dualphase1 && !inP1extra),
(MYBOOL) (minit == ITERATE_MINORRETRY), &candidatecount, &cviolated);
if(colnr < 0) {
minit = ITERATE_MAJORMAJOR;
continue;
}
#ifdef AcceptMarginalAccuracy
else if(xviolated+cviolated < lp->epspivot) {
if(lp->bb_trace || (lp->bb_totalnodes == 0))
report(lp, DETAILED, "dualloop: Assuming convergence with reduced accuracy %g.\n",
MAX(xviolated, cviolated));
rownr = 0;
colnr = 0;
}
#endif
if(lp->spx_status == FATHOMED)
break;
}
else
colnr = 0;
if(rownr > 0) {
if(colnr > 0) {
#ifdef Paranoia
if((rownr > lp->rows) || (colnr > lp->sum)) {
report(lp, SEVERE, "dualloop: Invalid row %d(%d) and column %d(%d) pair selected at iteration %.0f\n",
rownr, lp->rows, colnr-lp->columns, lp->columns, (double) get_total_iter(lp));
lp->spx_status = UNKNOWNERROR;
break;
}
#endif
fsolve(lp, colnr, pcol, workINT, lp->epsmachine, 1.0, TRUE);
#ifdef FixInaccurateDualMinit
if(colnr != minitcolnr)
minitcolnr = 0;
#endif
if(pcol[rownr] == 0) {
if(lp->spx_trace)
report(lp, DETAILED, "dualloop: Attempt to divide by zero (pcol[%d])\n", rownr);
if(!refactRecent(lp)) {
report(lp, DETAILED, "...trying to recover by refactorizing basis.\n");
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = FALSE;
}
else {
if(lp->bb_totalnodes == 0)
report(lp, DETAILED, "...cannot recover by refactorizing basis.\n");
lp->spx_status = NUMFAILURE;
ok = FALSE;
}
}
else {
set_action(&lp->spx_action, ACTION_ITERATE);
lp->rejectpivot[0] = 0;
if(!lp->obj_in_basis)
pcol[0] = my_chsign(!lp->is_lower[colnr], drow[colnr]);
theta = lp->bfp_prepareupdate(lp, rownr, colnr, pcol);
if((lp->improve & IMPROVE_THETAGAP) && !refactRecent(lp) &&
(my_reldiff(fabs(theta), fabs(prow[colnr])) >
lp->epspivot*10.0*log(2.0+50.0*lp->rows))) {
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
#ifdef IncreasePivotOnReducedAccuracy
lp->epspivot = MIN(1.0e-4, lp->epspivot*2.0);
#endif
report(lp, DETAILED, "dualloop: Refactorizing at iter %.0f due to loss of accuracy.\n",
(double) get_total_iter(lp));
}
theta = prow[colnr];
compute_theta(lp, rownr, &theta, !lp->is_lower[colnr], 0, primal);
}
}
#ifdef FixInaccurateDualMinit
else if(!refactRecent(lp) && (minit != ITERATE_MAJORMAJOR) && (colnr != minitcolnr)) {
minitcolnr = colnr;
i = invert(lp, INITSOL_USEZERO, TRUE);
if((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT))
break;
else if(!i) {
lp->spx_status = SINGULAR_BASIS;
break;
}
minit = ITERATE_MAJORMAJOR;
continue;
}
#endif
else {
if(!refactRecent(lp)) {
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
#ifdef dual_UseRejectionList
else if(lp->rejectpivot[0] < DEF_MAXPIVOTRETRY) {
lp->spx_status = RUNNING;
lp->rejectpivot[0]++;
lp->rejectpivot[lp->rejectpivot[0]] = rownr;
if(lp->bb_totalnodes == 0)
report(lp, DETAILED, "...trying to find another pivot row!\n");
goto RetryRow;
}
#endif
else if(dualphase1 && (dualoffset != 0)) {
lp->spx_status = LOSTFEAS;
if((lp->spx_trace && (lp->bb_totalnodes == 0)) ||
(lp->bb_trace && (lp->bb_totalnodes > 0)))
report(lp, DETAILED, "dualloop: Model lost dual feasibility.\n");
ok = FALSE;
break;
}
else {
if(lp->spx_status == RUNNING) {
#if 1
if(xviolated < lp->epspivot) {
if(lp->bb_trace || (lp->bb_totalnodes == 0))
report(lp, NORMAL, "The model is primal optimal, but marginally infeasible.\n");
lp->spx_status = OPTIMAL;
break;
}
#endif
lp->spx_status = INFEASIBLE;
if((lp->spx_trace && (lp->bb_totalnodes == 0)) ||
(lp->bb_trace && (lp->bb_totalnodes > 0)))
report(lp, DETAILED, "The model is primal infeasible.\n");
}
ok = FALSE;
break;
}
}
}
else if(inP1extra && !refactRecent(lp) && is_action(lp->improve, IMPROVE_INVERSE)) {
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
else {
bfpfinal = TRUE;
#ifdef dual_RemoveBasicFixedVars
if(inP1extra && (colnr == 0) && (lp->fixedvars > 0) && is_anti_degen(lp, ANTIDEGEN_FIXEDVARS)) {
report(lp, DETAILED, "dualloop: Trying to pivot out %d fixed basic variables at iter %.0f\n",
lp->fixedvars, (double) get_total_iter(lp));
rownr = 0;
while(lp->fixedvars > 0) {
rownr = findBasicFixedvar(lp, rownr, TRUE);
if(rownr == 0) {
colnr = 0;
break;
}
colnr = find_rowReplacement(lp, rownr, prow, nzprow);
if(colnr > 0) {
theta = 0;
performiteration(lp, rownr, colnr, theta, TRUE, FALSE, prow, NULL,
NULL, NULL, NULL);
lp->fixedvars--;
}
}
}
#endif
if(inP1extra && (colnr < 0) && !isPrimalFeasible(lp, lp->epsprimal, NULL, NULL)) {
if(lp->bb_totalnodes == 0) {
if(dualfeasible)
report(lp, DETAILED, "The model is primal infeasible and dual feasible.\n");
else
report(lp, DETAILED, "The model is primal infeasible and dual unbounded.\n");
}
set_OF_p1extra(lp, 0);
inP1extra = FALSE;
set_action(&lp->spx_action, ACTION_REINVERT);
lp->spx_status = INFEASIBLE;
lp->simplex_mode = SIMPLEX_UNDEFINED;
ok = FALSE;
}
else if(inP1extra) {
if(lp->bb_totalnodes == 0) {
report(lp, NORMAL, "Found feasibility by dual simplex after %10.0f iter.\n",
(double) get_total_iter(lp));
if((lp->usermessage != NULL) && (lp->msgmask & MSG_LPFEASIBLE))
lp->usermessage(lp, lp->msghandle, MSG_LPFEASIBLE);
}
set_OF_p1extra(lp, 0);
inP1extra = FALSE;
set_action(&lp->spx_action, ACTION_REINVERT);
#if 1
if((lp->simplex_strategy & SIMPLEX_DUAL_PRIMAL) && (lp->fixedvars == 0))
lp->spx_status = SWITCH_TO_PRIMAL;
#endif
changedphase = TRUE;
}
else {
lp->simplex_mode = SIMPLEX_Phase2_DUAL;
if((lp->fixedvars > 0) && (lp->bb_totalnodes == 0)) {
#ifdef dual_Phase1PriceEqualities
if(forceoutEQ != TRUE) {
forceoutEQ = TRUE;
goto RetryRow;
}
#endif
#ifdef Paranoia
report(lp, NORMAL,
#else
report(lp, DETAILED,
#endif
"Found dual solution with %d fixed slack variables left basic.\n",
lp->fixedvars);
}
colnr = 0;
if((dualoffset != 0) || (lp->bb_level <= 1) || (lp->improve & IMPROVE_BBSIMPLEX) || (lp->bb_rule & NODE_RCOSTFIXING)) {
set_action(&lp->piv_strategy, PRICE_FORCEFULL);
colnr = colprim(lp, drow, nzdrow, FALSE, 1, &candidatecount, FALSE, NULL);
clear_action(&lp->piv_strategy, PRICE_FORCEFULL);
if((dualoffset == 0) && (colnr > 0)) {
lp->spx_status = LOSTFEAS;
if(lp->total_iter == 0)
report(lp, DETAILED, "Recovering lost dual feasibility at iter %10.0f.\n",
(double) get_total_iter(lp));
break;
}
}
if(colnr == 0)
lp->spx_status = OPTIMAL;
else {
lp->spx_status = SWITCH_TO_PRIMAL;
if(lp->total_iter == 0)
report(lp, DETAILED, "Use primal simplex for finalization at iter %10.0f.\n",
(double) get_total_iter(lp));
}
if((lp->total_iter == 0) && (lp->spx_status == OPTIMAL))
report(lp, DETAILED, "Optimal solution with dual simplex at iter %10.0f.\n",
(double) get_total_iter(lp));
}
if(!changedphase)
break;
}
if(is_action(lp->spx_action, ACTION_ITERATE)) {
lastnr = lp->var_basic[rownr];
if(refactRecent(lp) == AUTOMATIC)
minitcount = 0;
else if(minitcount > MAX_MINITUPDATES) {
recompute_solution(lp, INITSOL_USEZERO);
minitcount = 0;
}
minit = performiteration(lp, rownr, colnr, theta, primal,
(MYBOOL) (
(stallaccept != AUTOMATIC)),
prow, nzprow,
pcol, NULL, boundswaps);
if(!lp->is_strongbranch && (lp->solutioncount >= 1) && !lp->spx_perturbed && !inP1extra &&
bb_better(lp, OF_WORKING, OF_TEST_WE)) {
lp->spx_status = FATHOMED;
ok = FALSE;
break;
}
if(minit != ITERATE_MAJORMAJOR)
minitcount++;
if(longsteps && dualphase1 && !inP1extra) {
dualfeasible = isDualFeasible(lp, lp->epsprimal, NULL, dualinfeasibles, NULL);
if(dualfeasible) {
dualphase1 = FALSE;
changedphase = TRUE;
lp->simplex_mode = SIMPLEX_Phase2_DUAL;
}
}
#ifdef UseDualReducedCostUpdate
else if(minit == ITERATE_MAJORMAJOR)
update_reducedcosts(lp, primal, lastnr, colnr, prow, drow);
#endif
if((minit == ITERATE_MAJORMAJOR) && (lastnr <= lp->rows) && is_fixedvar(lp, lastnr))
lp->fixedvars--;
}
if(lp->bfp_mustrefactorize(lp)) {
if(invert(lp, INITSOL_USEZERO, bfpfinal)) {
#if 0#endif
bfpfinal = FALSE;
#ifdef ResetMinitOnReinvert
minit = ITERATE_MAJORMAJOR;
#endif
}
else
lp->spx_status = SINGULAR_BASIS;
}
}
Finish:
stallMonitor_finish(lp);
multi_free(&(lp->longsteps));
FREE(prow);
FREE(nzprow);
FREE(pcol);
return(ok);
}
STATIC int spx_run(lprec *lp, MYBOOL validInvB)
{
int i, j, singular_count, lost_feas_count, *infeasibles = NULL, *boundflip_count;
MYBOOL primalfeasible, dualfeasible, lost_feas_state, isbb;
REAL primaloffset = 0, dualoffset = 0;
lp->current_iter = 0;
lp->current_bswap = 0;
lp->spx_status = RUNNING;
lp->bb_status = lp->spx_status;
lp->P1extraDim = 0;
set_OF_p1extra(lp, 0);
singular_count = 0;
lost_feas_count = 0;
lost_feas_state = FALSE;
lp->simplex_mode = SIMPLEX_DYNAMIC;
lp->fixedvars = 0;
lp->boundedvars = 0;
for(i = 1; i <= lp->rows; i++) {
j = lp->var_basic[i];
if((j <= lp->rows) && is_fixedvar(lp, j))
lp->fixedvars++;
if((lp->upbo[i] < lp->infinite) && (lp->upbo[i] > lp->epsprimal))
lp->boundedvars++;
}
for(; i <= lp->sum; i++){
if((lp->upbo[i] < lp->infinite) && (lp->upbo[i] > lp->epsprimal))
lp->boundedvars++;
}
#ifdef UseLongStepDualPhase1
allocINT(lp, &infeasibles, lp->columns + 1, FALSE);
infeasibles[0] = 0;
#endif
isbb = (MYBOOL) ((MIP_count(lp) > 0) && (lp->bb_level > 1));
if(is_action(lp->spx_action, ACTION_REINVERT)) {
if(isbb && (lp->bb_bounds->nodessolved == 0))
recompute_solution(lp, INITSOL_SHIFTZERO);
else {
i = my_if(is_action(lp->spx_action, ACTION_REBASE), INITSOL_SHIFTZERO, INITSOL_USEZERO);
invert(lp, (MYBOOL) i, TRUE);
}
}
else if(is_action(lp->spx_action, ACTION_REBASE))
recompute_solution(lp, INITSOL_SHIFTZERO);
if(is_action(lp->improve, IMPROVE_DUALFEAS) || (lp->rows == 0))
boundflip_count = &i;
else
boundflip_count = NULL;
while(lp->spx_status == RUNNING) {
dualfeasible = isbb ||
isDualFeasible(lp, lp->epsprimal, boundflip_count, infeasibles, &dualoffset);
if(is_action(lp->spx_action, ACTION_RECOMPUTE))
recompute_solution(lp, INITSOL_USEZERO);
primalfeasible = isPrimalFeasible(lp, lp->epsprimal, NULL, &primaloffset);
if(userabort(lp, -1))
break;
if(lp->spx_trace) {
if(primalfeasible)
report(lp, NORMAL, "Start at primal feasible basis\n");
else if(dualfeasible)
report(lp, NORMAL, "Start at dual feasible basis\n");
else if(lost_feas_count > 0)
report(lp, NORMAL, "Continuing at infeasible basis\n");
else
report(lp, NORMAL, "Start at infeasible basis\n");
}
if(((lp->simplex_strategy & SIMPLEX_Phase1_DUAL) == 0) ||
((MIP_count(lp) > 0) && (lp->total_iter == 0) &&
is_presolve(lp, PRESOLVE_REDUCEMIP))) {
if(!lost_feas_state && primalfeasible && ((lp->simplex_strategy & SIMPLEX_Phase2_DUAL) > 0))
lp->spx_status = SWITCH_TO_DUAL;
else
primloop(lp, primalfeasible, 0.0);
if(lp->spx_status == SWITCH_TO_DUAL)
dualloop(lp, TRUE, NULL, 0.0);
}
else {
if(!lost_feas_state && primalfeasible && ((lp->simplex_strategy & SIMPLEX_Phase2_PRIMAL) > 0))
lp->spx_status = SWITCH_TO_PRIMAL;
else
dualloop(lp, dualfeasible, infeasibles, dualoffset);
if(lp->spx_status == SWITCH_TO_PRIMAL)
primloop(lp, TRUE, 0.0);
}
i = lp->spx_status;
primalfeasible = (MYBOOL) (i == OPTIMAL);
if(primalfeasible || (i == UNBOUNDED))
break;
else if(((i == INFEASIBLE) && is_anti_degen(lp, ANTIDEGEN_INFEASIBLE)) ||
((i == LOSTFEAS) && is_anti_degen(lp, ANTIDEGEN_LOSTFEAS)) ||
((i == NUMFAILURE) && is_anti_degen(lp, ANTIDEGEN_NUMFAILURE)) ||
((i == DEGENERATE) && is_anti_degen(lp, ANTIDEGEN_STALLING))) {
if((lp->bb_level <= 1) || is_anti_degen(lp, ANTIDEGEN_DURINGBB))
break;
#ifdef AssumeHighAccuracyInBB
if((lp->bb_level > 1) && (i == INFEASIBLE))
break;
#endif
}
if(lp->spx_status == SINGULAR_BASIS) {
lost_feas_state = FALSE;
singular_count++;
if(singular_count >= DEF_MAXSINGULARITIES) {
report(lp, IMPORTANT, "spx_run: Failure due to too many singular bases.\n");
lp->spx_status = NUMFAILURE;
break;
}
if(lp->spx_trace || (lp->verbose > DETAILED))
report(lp, NORMAL, "spx_run: Singular basis; attempting to recover.\n");
lp->spx_status = RUNNING;
}
else {
lost_feas_state = (MYBOOL) (lp->spx_status == LOSTFEAS);
#if 0#endif
if(lost_feas_state) {
lost_feas_count++;
if(lost_feas_count < DEF_MAXSINGULARITIES) {
report(lp, DETAILED, "spx_run: Recover lost feasibility at iter %10.0f.\n",
(double) get_total_iter(lp));
lp->spx_status = RUNNING;
}
else {
report(lp, IMPORTANT, "spx_run: Lost feasibility %d times - iter %10.0f and %9.0f nodes.\n",
lost_feas_count, (double) get_total_iter(lp), (double) lp->bb_totalnodes);
lp->spx_status = NUMFAILURE;
}
}
}
}
lp->total_iter += lp->current_iter;
lp->current_iter = 0;
lp->total_bswap += lp->current_bswap;
lp->current_bswap = 0;
FREE(infeasibles);
return(lp->spx_status);
}
lprec *make_lag(lprec *lpserver)
{
int i;
lprec *hlp;
MYBOOL ret;
REAL *duals;
hlp = make_lp(0, lpserver->columns);
if(hlp != NULL) {
set_sense(hlp, is_maxim(lpserver));
hlp->lag_bound = lpserver->bb_limitOF;
for(i = 1; i <= lpserver->columns; i++) {
set_mat(hlp, 0, i, get_mat(lpserver, 0, i));
if(is_binary(lpserver, i))
set_binary(hlp, i, TRUE);
else {
set_int(hlp, i, is_int(lpserver, i));
set_bounds(hlp, i, get_lowbo(lpserver, i), get_upbo(lpserver, i));
}
}
hlp->matL = lpserver->matA;
inc_lag_space(hlp, lpserver->rows, TRUE);
ret = get_ptr_sensitivity_rhs(hlp, &duals, NULL, NULL);
for(i = 1; i <= lpserver->rows; i++) {
hlp->lag_con_type[i] = get_constr_type(lpserver, i);
hlp->lag_rhs[i] = lpserver->orig_rhs[i];
hlp->lambda[i] = (ret) ? duals[i - 1] : 0.0;
}
}
return(hlp);
}
STATIC int heuristics(lprec *lp, int mode)
{
lprec *hlp;
int status = PROCFAIL;
if(lp->bb_level > 1)
return( status );
status = RUNNING;
lp->bb_limitOF = my_chsign(is_maxim(lp), -lp->infinite);
if(FALSE && (lp->int_vars > 0)) {
hlp = make_lag(lp);
status = solve(hlp);
lp->bb_heuristicOF = hlp->best_solution[0];
hlp->matL = NULL;
delete_lp(hlp);
}
lp->timeheuristic = timeNow();
return( status );
}
STATIC int lag_solve(lprec *lp, REAL start_bound, int num_iter)
{
int i, j, citer, nochange, oldpresolve;
MYBOOL LagFeas, AnyFeas, Converged, same_basis;
REAL *OrigObj, *ModObj, *SubGrad, *BestFeasSol;
REAL Zub, Zlb, Znow, Zprev, Zbest, rhsmod, hold;
REAL Phi, StepSize = 0.0, SqrsumSubGrad;
if(lp->spx_status != OPTIMAL) {
lp->lag_status = NOTRUN;
return( lp->lag_status );
}
if(!allocREAL(lp, &OrigObj, lp->columns + 1, FALSE) ||
!allocREAL(lp, &ModObj, lp->columns + 1, TRUE) ||
!allocREAL(lp, &SubGrad, get_Lrows(lp) + 1, TRUE) ||
!allocREAL(lp, &BestFeasSol, lp->sum + 1, TRUE)) {
lp->lag_status = NOMEMORY;
return( lp->lag_status );
}
lp->lag_status = RUNNING;
oldpresolve = lp->do_presolve;
lp->do_presolve = PRESOLVE_NONE;
push_basis(lp, NULL, NULL, NULL);
Zlb = lp->best_solution[0];
Zub = start_bound;
Zbest = Zub;
Znow = Zlb;
Zprev = lp->infinite;
rhsmod = 0;
Phi = DEF_LAGCONTRACT;
LagFeas = FALSE;
Converged= FALSE;
AnyFeas = FALSE;
citer = 0;
nochange = 0;
get_row(lp, 0, OrigObj);
#ifdef DirectOverrideOF
set_OF_override(lp, ModObj);
#endif
OrigObj[0] = get_rh(lp, 0);
for(i = 1 ; i <= get_Lrows(lp); i++)
lp->lambda[i] = 0;
while((lp->lag_status == RUNNING) && (citer < num_iter)) {
citer++;
LagFeas = TRUE;
Converged = TRUE;
SqrsumSubGrad = 0;
for(i = 1; i <= get_Lrows(lp); i++) {
hold = lp->lag_rhs[i];
for(j = 1; j <= lp->columns; j++)
hold -= mat_getitem(lp->matL, i, j) * lp->best_solution[lp->rows + j];
if(LagFeas) {
if(lp->lag_con_type[i] == EQ) {
if(fabs(hold) > lp->epsprimal)
LagFeas = FALSE;
}
else if(hold < -lp->epsprimal)
LagFeas = FALSE;
}
if(Converged && (fabs(my_reldiff(hold , SubGrad[i])) > lp->lag_accept))
Converged = FALSE;
SubGrad[i] = hold;
SqrsumSubGrad += hold * hold;
}
SqrsumSubGrad = sqrt(SqrsumSubGrad);
#if 1
Converged &= LagFeas;
#endif
if(Converged)
break;
Znow = lp->best_solution[0] - rhsmod;
if(Znow > Zub) {
Phi *= DEF_LAGCONTRACT;
StepSize *= (Zub-Zlb) / (Znow-Zlb);
}
else
#define LagBasisContract
#ifdef LagBasisContract
StepSize = Phi * (2-DEF_LAGCONTRACT) * (Zub - Znow) / SqrsumSubGrad;
#else
StepSize = Phi * (Zub - Znow) / SqrsumSubGrad;
#endif
for(i = 1; i <= get_Lrows(lp); i++) {
lp->lambda[i] += StepSize * SubGrad[i];
if((lp->lag_con_type[i] != EQ) && (lp->lambda[i] > 0)) {
if(Znow < Zub)
lp->lambda[i] = 0;
}
}
if(LagFeas && (Znow < Zbest)) {
MEMCOPY(BestFeasSol, lp->best_solution, lp->sum+1);
hold = OrigObj[0];
for(i = 1; i <= lp->columns; i++)
hold += lp->best_solution[lp->rows + i] * OrigObj[i];
BestFeasSol[0] = hold;
if(lp->lag_trace)
report(lp, NORMAL, "lag_solve: Improved feasible solution at iteration %d of %g\n",
citer, hold);
Zbest = Znow;
AnyFeas = TRUE;
nochange = 0;
}
else if(Znow == Zprev) {
nochange++;
if(nochange > LAG_SINGULARLIMIT) {
Phi *= 0.5;
nochange = 0;
}
}
Zprev = Znow;
for(j = 1; j <= lp->columns; j++) {
hold = 0;
for(i = 1; i <= get_Lrows(lp); i++)
hold += lp->lambda[i] * mat_getitem(lp->matL, i, j);
ModObj[j] = OrigObj[j] - my_chsign(is_maxim(lp), hold);
#ifndef DirectOverrideOF
set_mat(lp, 0, j, ModObj[j]);
#endif
}
rhsmod = my_chsign(is_maxim(lp), get_rh(lp, 0));
for(i = 1; i <= get_Lrows(lp); i++)
rhsmod += lp->lambda[i] * lp->lag_rhs[i];
if(lp->lag_trace) {
report(lp, IMPORTANT, "Zub: %10g Zlb: %10g Stepsize: %10g Phi: %10g Feas %d\n",
(double) Zub, (double) Zlb, (double) StepSize, (double) Phi, LagFeas);
for(i = 1; i <= get_Lrows(lp); i++)
report(lp, IMPORTANT, "%3d SubGrad %10g lambda %10g\n",
i, (double) SubGrad[i], (double) lp->lambda[i]);
if(lp->sum < 20)
print_lp(lp);
}
i = spx_solve(lp);
if(lp->spx_status == UNBOUNDED) {
if(lp->lag_trace) {
report(lp, NORMAL, "lag_solve: Unbounded solution encountered with this OF:\n");
for(i = 1; i <= lp->columns; i++)
report(lp, NORMAL, RESULTVALUEMASK " ", (double) ModObj[i]);
}
goto Leave;
}
else if((lp->spx_status == NUMFAILURE) || (lp->spx_status == PROCFAIL) ||
(lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) ||
(lp->spx_status == INFEASIBLE)) {
lp->lag_status = lp->spx_status;
}
#ifdef LagBasisContract
same_basis = compare_basis(lp);
if(LagFeas &&
!same_basis) {
pop_basis(lp, FALSE);
push_basis(lp, NULL, NULL, NULL);
Phi *= DEF_LAGCONTRACT;
}
if(lp->lag_trace) {
report(lp, DETAILED, "lag_solve: Simplex status code %d, same basis %s\n",
lp->spx_status, my_boolstr(same_basis));
print_solution(lp, 1);
}
#endif
}
if(AnyFeas) {
lp->lag_bound = my_chsign(is_maxim(lp), Zbest);
for(i = 0; i <= lp->sum; i++)
lp->solution[i] = BestFeasSol[i];
transfer_solution(lp, TRUE);
if(!is_maxim(lp))
for(i = 1; i <= get_Lrows(lp); i++)
lp->lambda[i] = my_flipsign(lp->lambda[i]);
}
Leave:
if(citer >= num_iter) {
if(AnyFeas)
lp->lag_status = FEASFOUND;
else
lp->lag_status = NOFEASFOUND;
}
else
lp->lag_status = lp->spx_status;
if(lp->lag_status == OPTIMAL) {
report(lp, NORMAL, "\nLagrangean convergence achieved in %d iterations\n", citer);
i = check_solution(lp, lp->columns,
lp->best_solution, lp->orig_upbo, lp->orig_lowbo, lp->epssolution);
}
else {
report(lp, NORMAL, "\nUnsatisfactory convergence achieved over %d Lagrangean iterations.\n",
citer);
if(AnyFeas)
report(lp, NORMAL, "The best feasible Lagrangean objective function value was %g\n",
lp->best_solution[0]);
}
#ifdef DirectOverrideOF
set_OF_override(lp, NULL);
#else
for(i = 1; i <= lp->columns; i++)
set_mat(lp, 0, i, OrigObj[i]);
#endif
FREE(BestFeasSol);
FREE(SubGrad);
FREE(OrigObj);
FREE(ModObj);
pop_basis(lp, FALSE);
lp->do_presolve = oldpresolve;
return( lp->lag_status );
}
STATIC int spx_solve(lprec *lp)
{
int status;
MYBOOL iprocessed;
lp->total_iter = 0;
lp->total_bswap = 0;
lp->perturb_count = 0;
lp->bb_maxlevel = 1;
lp->bb_totalnodes = 0;
lp->bb_improvements = 0;
lp->bb_strongbranches= 0;
lp->is_strongbranch = FALSE;
lp->bb_level = 0;
lp->bb_solutionlevel = 0;
lp->best_solution[0] = my_chsign(is_maxim(lp), lp->infinite);
if(lp->invB != NULL)
lp->bfp_restart(lp);
lp->spx_status = presolve(lp);
if(lp->spx_status == PRESOLVED) {
status = lp->spx_status;
goto Reconstruct;
}
else if(lp->spx_status != RUNNING)
goto Leave;
iprocessed = !lp->wasPreprocessed;
if(!preprocess(lp) || userabort(lp, -1))
goto Leave;
if(mat_validate(lp->matA)) {
lp->solutioncount = 0;
lp->real_solution = lp->infinite;
set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT);
lp->bb_break = FALSE;
status = run_BB(lp);
if(iprocessed)
postprocess(lp);
Reconstruct:
if(!postsolve(lp, status))
report(lp, SEVERE, "spx_solve: Failure during postsolve.\n");
goto Leave;
}
if(lp->bb_trace || lp->spx_trace)
report(lp, CRITICAL, "spx_solve: The current LP seems to be invalid\n");
lp->spx_status = NUMFAILURE;
Leave:
lp->timeend = timeNow();
if((lp->lag_status != RUNNING) && (lp->invB != NULL)) {
int itemp;
REAL test;
itemp = lp->bfp_nonzeros(lp, TRUE);
test = 100;
if(lp->total_iter > 0)
test *= (REAL) lp->total_bswap/lp->total_iter;
report(lp, NORMAL, "\n ");
report(lp, NORMAL, "MEMO: lp_solve version %d.%d.%d.%d for %d bit OS, with %d bit REAL variables.\n",
MAJORVERSION, MINORVERSION, RELEASE, BUILD, 8*sizeof(void *), 8*sizeof(REAL));
report(lp, NORMAL, " In the total iteration count %.0f, %.0f (%.1f%%) were bound flips.\n",
(double) lp->total_iter, (double) lp->total_bswap, test);
report(lp, NORMAL, " There were %d refactorizations, %d triggered by time and %d by density.\n",
lp->bfp_refactcount(lp, BFP_STAT_REFACT_TOTAL),
lp->bfp_refactcount(lp, BFP_STAT_REFACT_TIMED),
lp->bfp_refactcount(lp, BFP_STAT_REFACT_DENSE));
report(lp, NORMAL, " ... on average %.1f major pivots per refactorization.\n",
get_refactfrequency(lp, TRUE));
report(lp, NORMAL, " The largest [%s] fact(B) had %d NZ entries, %.1fx largest basis.\n",
lp->bfp_name(), itemp, lp->bfp_efficiency(lp));
if(lp->perturb_count > 0)
report(lp, NORMAL, " The bounds were relaxed via perturbations %d times.\n",
lp->perturb_count);
if(MIP_count(lp) > 0) {
if(lp->bb_solutionlevel > 0)
report(lp, NORMAL, " The maximum B&B level was %d, %.1fx MIP order, %d at the optimal solution.\n",
lp->bb_maxlevel, (double) lp->bb_maxlevel / (MIP_count(lp)+lp->int_vars), lp->bb_solutionlevel);
else
report(lp, NORMAL, " The maximum B&B level was %d, %.1fx MIP order, with %.0f nodes explored.\n",
lp->bb_maxlevel, (double) lp->bb_maxlevel / (MIP_count(lp)+lp->int_vars), (double) get_total_nodes(lp));
if(GUB_count(lp) > 0)
report(lp, NORMAL, " %d general upper-bounded (GUB) structures were employed during B&B.\n",
GUB_count(lp));
}
report(lp, NORMAL, " The constraint matrix inf-norm is %g, with a dynamic range of %g.\n",
lp->matA->infnorm, lp->matA->dynrange);
report(lp, NORMAL, " Time to load data was %.3f seconds, presolve used %.3f seconds,\n",
lp->timestart-lp->timecreate, lp->timepresolved-lp->timestart);
report(lp, NORMAL, " ... %.3f seconds in simplex solver, in total %.3f seconds.\n",
lp->timeend-lp->timepresolved, lp->timeend-lp->timecreate);
}
return( lp->spx_status );
}
int lin_solve(lprec *lp)
{
int status = NOTRUN;
lp->lag_status = NOTRUN;
if(lp->columns == 0) {
default_basis(lp);
lp->spx_status = NOTRUN;
return( lp->spx_status);
}
unset_OF_p1extra(lp);
free_duals(lp);
FREE(lp->drow);
FREE(lp->nzdrow);
if(lp->bb_cuttype != NULL)
freecuts_BB(lp);
lp->timestart = timeNow();
lp->timeheuristic = 0;
lp->timepresolved = 0;
lp->timeend = 0;
if(heuristics(lp, AUTOMATIC) != RUNNING)
return( INFEASIBLE );
status = spx_solve(lp);
if((get_Lrows(lp) > 0) && (lp->lag_status == NOTRUN)) {
if(status == OPTIMAL)
status = lag_solve(lp, lp->bb_heuristicOF, DEF_LAGMAXITERATIONS);
else
report(lp, IMPORTANT, "\nCannot do Lagrangean optimization since root model was not solved.\n");
}
lp->bb_heuristicOF = my_chsign(is_maxim(lp), lp->infinite);
if((lp->spx_status == OPTIMAL) && (lp->bb_totalnodes > 0)) {
if((lp->bb_break && !bb_better(lp, OF_DUALLIMIT, OF_TEST_BE)) )
status = lp->spx_status = SUBOPTIMAL;
}
return( status );
}