s14h
/
prism-accumulation
1
0
Fork 0
Code Releases Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3738 Commits
5 Branches
0 Tags
31 MiB
 
 
 
 
 
 
Tree: bc8f512cf7
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'bc8f512cf7'
${ noResults }
prism-accumulation/prism/ext/lpsolve55/src/lp_solve_5.5/lp_simplex.c

2206 lines
75 KiB
Raw Blame History

/*
Core optimization drivers for lp_solve v5.0+
----------------------------------------------------------------------------------
Author: Michel Berkelaar (to lp_solve v3.2),
Kjell Eikland (v4.0 and forward)
Contact:
License terms: LGPL.
Requires: lp_lib.h, lp_simplex.h, lp_presolve.h, lp_pricerPSE.h
Release notes:
v5.0.0 1 January 2004 New unit applying stacked basis and bounds storage.
v5.0.1 31 January 2004 Moved B&B routines to separate file and implemented
a new runsolver() general purpose call method.
v5.0.2 1 May 2004 Changed routine names to be more intuitive.
v5.1.0 10 January 2005 Created modular stalling/cycling functions.
Rewrote dualloop() to optimize long dual and
also streamlined primloop() correspondingly.
v5.2.0 20 March 2005 Reimplemented primal phase 1 logic.
Made multiple pricing finally work (primal simplex).
----------------------------------------------------------------------------------
*/
#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; /* Expand degeneracy/stalling tolerance range */
#endif
monitor->limitstall[TRUE] = monitor->limitstall[FALSE];
if(monitor->oldpivrule == PRICER_DEVEX) /* Increase tolerance since primal Steepest Edge is expensive */
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->epsvalue; */
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;
/* Accept unconditionally if this is the first or second call */
monitor->active = FALSE;
if(monitor->Icount <= 1) {
if(monitor->Icount == 1) {
monitor->prevobj = lp->rhs[0];
monitor->previnfeas = deltaobj;
}
monitor->Icount++;
return( acceptance );
}
/* Define progress as primal objective less sum of (primal/dual) infeasibilities */
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);
/* Check if we have a stationary solution at selected tolerance level;
allow some difference in case we just refactorized the basis. */
deltaobj = my_reldiff(monitor->thisobj, monitor->prevobj);
deltaobj = fabs(deltaobj); /* Pre v5.2 version */
isStalled = (MYBOOL) (deltaobj < monitor->epsvalue);
/* Also require that we have a measure of infeasibility-stalling */
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
refvalue *= 1000*log10(9.0+lp->sum);
#endif
testvalue = my_reldiff(monitor->thisinfeas, monitor->previnfeas);
isStalled &= (fabs(testvalue) < refvalue);
/* Check if we should force "major" pivoting, i.e. no bound flips;
this is activated when we see the feasibility deteriorate */
/* if(!isStalled && (testvalue > 0) && (TRUE || is_action(lp->anti_degen, ANTIDEGEN_BOUNDFLIP))) */
#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
isCreeping |= stallMonitor_creepingObj(lp);
/* isCreeping |= stallMonitor_shortSteps(lp); */
#endif
if(isStalled || isCreeping) {
/* Update counters along with specific tolerance for bound flips */
#if 1
if(minit != ITERATE_MAJORMAJOR) {
if(++monitor->Mcycle > 2) {
monitor->Mcycle = 0;
monitor->Ncycle++;
}
}
else
#endif
monitor->Ncycle++;
/* Start to monitor for variable cycling if this is the initial stationarity */
if(monitor->Ncycle <= 1) {
monitor->Ccycle = colnr;
monitor->Rcycle = rownr;
}
/* Check if we should change pivoting strategy */
else if(isCreeping || /* We have OF creep */
(monitor->Ncycle > monitor->limitstall[monitor->isdual]) || /* KE empirical value */
((monitor->Ccycle == rownr) && (monitor->Rcycle == colnr))) { /* Obvious cycling */
monitor->active = TRUE;
/* Try to force out equality slacks to combat degeneracy */
if((lp->fixedvars > 0) && (*forceoutEQ != TRUE)) {
*forceoutEQ = TRUE;
goto Proceed;
}
/* Our options are now to select an alternative rule or to do bound perturbation;
check if these options are available to us or if we must signal failure and break out. */
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 );
}
/* See if we can do the appropriate alternative rule. */
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)) {
/* Switch rule to combat degeneracy. */
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);
}
/* If not, code for bound relaxation/perturbation */
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 );
}
}
}
/* Otherwise change back to original selection strategy as soon as possible */
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;
}
/* Update objective progress tracker */
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)
/* This routine is called for each constraint at the start of
primloop and the primal problem is infeasible. Its
purpose is to add artificial variables and associated
objective function values to populate primal phase 1. */
{
MYBOOL add;
/* Make sure we don't add unnecessary artificials, i.e. avoid
cases where the slack variable is enough */
add = !isBasisVarFeasible(lp, lp->epspivot, forrownr);
if(add) {
int *rownr = NULL, i, bvar, ii;
REAL *avalue = NULL, rhscoef, acoef;
MATrec *mat = lp->matA;
/* Check the simple case where a slack is basic */
for(i = 1; i <= lp->rows; i++) {
if(lp->var_basic[i] == forrownr)
break;
}
acoef = 1;
/* If not, look for any basic user variable that has a
non-zero coefficient in the current constraint row */
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 no candidate was found above, gamble on using the densest column available */
#if 0
if(i > lp->rows) {
int len = 0;
bvar = 0;
for(i = 1; i <= lp->rows; i++) {
ii = lp->var_basic[i] - lp->rows;
if((ii <= 0) || (ii > (lp->columns-lp->P1extraDim)))
continue;
if(mat_collength(mat, ii) > len) {
len = mat_collength(mat, ii);
bvar = i;
}
}
i = bvar;
acoef = 1;
}
#endif
bvar = i;
add = (MYBOOL) (bvar <= lp->rows);
if(add) {
rhscoef = lp->rhs[forrownr];
/* Create temporary sparse array storage */
if(nzarray == NULL)
allocREAL(lp, &avalue, 2, FALSE);
else
avalue = nzarray;
if(idxarray == NULL)
allocINT(lp, &rownr, 2, FALSE);
else
rownr = idxarray;
/* Set the objective coefficient */
rownr[0] = 0;
avalue[0] = my_chsign(is_chsign(lp, 0), 1);
/* Set the constraint row coefficient */
rownr[1] = forrownr;
avalue[1] = my_chsign(is_chsign(lp, forrownr), my_sign(rhscoef/acoef));
/* Add the column of artificial variable data to the user data matrix */
add_columnex(lp, 2, avalue, rownr);
/* Free the temporary sparse array storage */
if(idxarray == NULL)
FREE(rownr);
if(nzarray == NULL)
FREE(avalue);
/* Now set the artificial variable to be basic */
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
/* Return the row index of the singleton */
colnr = mat->col_end[colnr-1];
colnr = COL_MAT_ROWNR(colnr);
return( colnr );
}
STATIC int findAnti_artificial(lprec *lp, int colnr)
/* Primal simplex: Find a basic artificial variable to swap
against the non-basic slack variable, if possible */
{
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);
/* Should we find the artificial's slack direct "antibody"? */
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
performiteration(lp, rownr, colnr, 0.0, TRUE, FALSE, prow, NULL,
NULL, NULL, NULL);
#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;
/* Substitute any basic artificial variable for its slack counterpart */
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
/* Delete any remaining non-basic artificial variables */
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));
/* Add sufficent number of artificial variables to make the problem feasible
through the first phase; delete when primal feasibility has been achieved */
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
/* First check if we can get away with a single artificial variable */
if(lp->equalities == 0) {
i = (int) feasibilityOffset(lp, !primal);
add_artificial(lp, i, prow, (int *) pcol);
}
else
#endif
/* Otherwise add as many artificial variables as is necessary
to force primal feasibility. */
for(i = 1; i <= lp->rows; i++) {
add_artificial(lp, i, NULL, NULL);
}
/* Make sure we update the working objective */
if(lp->P1extraDim > 0) {
#if 1 /* v5.1 code: Not really necessary since we do not price the artificial
variables (stored at the end of the column list, they are initially
basic and are never allowed to enter the basis, once they exit) */
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);
}
/* Create work arrays and optionally the multiple pricing structure */
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;
/* Initialize regular primal simplex algorithm variables */
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;
/* Iterate while we are successful; exit when the model is infeasible/unbounded,
or we must terminate due to numeric instability or user-determined reasons */
while((lp->spx_status == RUNNING) && !userabort(lp, -1)) {
primalphase1 = (MYBOOL) (lp->P1extraDim > 0);
clear_action(&lp->spx_action, ACTION_REINVERT | ACTION_ITERATE);
/* Check if we have stalling (from numerics or degenerate cycling) */
pricerCanChange = !primalphase1;
stallaccept = stallMonitor_check(lp, rownr, colnr, lastnr, minit, pricerCanChange, &forceoutEQ);
if(!stallaccept)
break;
/* Find best column to enter the basis */
RetryCol:
#if 0
if(verify_solution(lp, FALSE, "spx_loop") > 0)
i = 1; /* This is just a debug trap */
#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));
/* Handle direct outcomes */
if(colnr == 0)
lp->spx_status = OPTIMAL;
if(lp->rejectpivot[0] > 0)
minit = ITERATE_MAJORMAJOR;
/* See if accuracy check during compute_reducedcosts flagged refactorization */
if(is_action(lp->spx_action, ACTION_REINVERT))
bfpfinal = TRUE;
}
/* Make sure that we do not erroneously conclude that an unbounded model is optimal */
#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
/* Check if we found an entering variable (indicating that we are still dual infeasible) */
if(colnr > 0) {
changedphase = FALSE;
fsolve(lp, colnr, pcol, NULL, lp->epsmachine, 1.0, TRUE); /* Solve entering column for Pi */
/* Do special anti-degeneracy column selection, if specified */
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");
}
}
/* Find the leaving variable that gives the most stringent bound on the entering variable */
theta = drow[colnr];
rownr = rowprim(lp, colnr, &theta, pcol, workINT, forceoutEQ, &cviolated);
#ifdef AcceptMarginalAccuracy
/* Check for marginal accuracy */
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
/* See if we can do a straight artificial<->slack replacement (when "colnr" is a slack) */
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) /* We must manually copy the reduced cost for RHS update */
pcol[0] = my_chsign(!lp->is_lower[colnr], drow[colnr]);
lp->bfp_prepareupdate(lp, rownr, colnr, pcol);
}
/* We may be unbounded... */
else {
/* First make sure that we are not suffering from precision loss */
#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
/* Check that we are not having numerical problems */
if(!refactRecent(lp) && !pendingunbounded) {
bfpfinal = TRUE;
pendingunbounded = TRUE;
set_action(&lp->spx_action, ACTION_REINVERT);
}
/* Conclude that the model is unbounded */
else {
lp->spx_status = UNBOUNDED;
report(lp, DETAILED, "The model is primal unbounded.\n");
break;
}
}
}
/* We handle optimality and phase 1 infeasibility ... */
else {
Optimality:
/* Handle possible transition from phase 1 to phase 2 */
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);
/* We can do two things now;
1) delete the rows belonging to those variables, since they are redundant, OR
2) drive out the existing artificial variables via pivoting. */
if(lp->P1extraDim > 0) {
#ifdef Phase1EliminateRedundant
/* If it is not a MIP model we can try to delete redundant rows */
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
/* Obtain column and row indeces */
j = lp->var_basic[i]-lp->rows;
k = get_artificialRow(lp, j);
/* Delete row before column due to basis "compensation logic" */
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;
}
/* Otherwise we drive out the artificials by elimination pivoting */
else
eliminate_artificials(lp, prow);
#else
/* Indicate phase 2 with artificial variables by negating P1extraDim */
lp->P1extraDim = my_flipsign(lp->P1extraDim);
#endif
}
/* We must refactorize since the OF changes from phase 1 to phase 2 */
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
/* We are infeasible in phase 1 */
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));
}
}
/* Handle phase 1 optimality */
else {
/* (Do nothing special) */
}
/* Check if we are still primal feasible; the default assumes that this check
is not necessary after the relaxed problem has been solved satisfactorily. */
if((lp->bb_level <= 1) || (lp->improve & IMPROVE_BBSIMPLEX) /* || (lp->bb_rule & NODE_RCOSTFIXING) */) { /* NODE_RCOSTFIXING fix */
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));
}
}
}
/* Pivot row/col and update the inverse */
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) (/*(candidatecount > 1) && */
(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
/* Do a fast update of the reduced costs in preparation for the next iteration */
if(minit == ITERATE_MAJORMAJOR)
update_reducedcosts(lp, primal, lastnr, colnr, pcol, drow);
#endif
/* Detect if an auxiliary variable has left the basis and delete it; if
the non-basic variable only changed bound (a "minor iteration"), the
basic artificial variable did not leave and there is nothing to do */
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;
}
}
/* Remove any remaining artificial variables (feasible or infeasible model) */
lp->P1extraDim = abs(lp->P1extraDim);
/* if((lp->P1extraDim > 0) && (lp->spx_status != DEGENERATE)) { */
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
/* Switch to dual phase 1 simplex for MIP models during
B&B phases, since this is typically far more efficient */
#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);
} /* primloop */
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));
/* Allocate work arrays */
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;
/* Set non-zero P1extraVal value to force dual feasibility when the dual
simplex is used as a phase 1 algorithm for the primal simplex.
The value will be reset when primal feasibility has been achieved, or
a dual non-feasibility has been encountered (no candidate for a first
leaving variable) */
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);
/* Prepare dual long-step structures */
#if 0
longsteps = TRUE;
#elif 0
longsteps = (MYBOOL) ((MIP_count(lp) > 0) && (lp->bb_level > 1));
#elif 0
longsteps = (MYBOOL) ((MIP_count(lp) > 0) && (lp->solutioncount >= 1));
#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);
}
/* Do regular dual simplex variable initializations */
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;
/* Check if we have equality slacks in the basis and we should try to
drive them out in order to reduce chance of degeneracy in Phase 1.
forceoutEQ = FALSE : Only eliminate assured "good" violated
equality constraint slacks
AUTOMATIC: Seek more elimination of equality constraint
slacks (but not as aggressive as the rule
used in lp_solve v4.0 and earlier)
TRUE: Force remaining equality slacks out of the
basis */
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;
/* case AUTOMATIC: forceoutEQ = TRUE;
break;
default: forceoutEQ = TRUE; */
}
}
#endif
while((lp->spx_status == RUNNING) && !userabort(lp, -1)) {
/* Check if we have stalling (from numerics or degenerate cycling) */
pricerCanChange = !dualphase1 && !inP1extra;
stallaccept = stallMonitor_check(lp, rownr, colnr, lastnr, minit, pricerCanChange, &forceoutEQ);
if(!stallaccept)
break;
/* Store current LP index for reference at next iteration */
changedphase = FALSE;
/* Compute (pure) dual phase1 offsets / reduced costs if appropriate */
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);
/* Solve (note that solved pcol will be used instead of lp->rhs) */
ftran(lp, pcol, NULL, lp->epsmachine);
}
/* Do minor iterations (non-basic variable bound flips) for as
long as possible since this is a cheap way of iterating */
#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));
}
/* Make sure that we do not erroneously conclude that an infeasible model is optimal */
#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
/* If we found a leaving variable, find a matching entering one */
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
/* Check if the long-dual found reason to prune the B&B tree */
if(lp->spx_status == FATHOMED)
break;
}
else
colnr = 0;
/* Process primal-infeasible row */
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
/* Prevent bound flip-flops during minor iterations; used to detect
infeasibility after triggering of minor iteration accuracy management */
if(colnr != minitcolnr)
minitcolnr = 0;
#endif
/* Getting division by zero here; catch it and try to recover */
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) /* We must manually copy the reduced cost for RHS update */
pcol[0] = my_chsign(!lp->is_lower[colnr], drow[colnr]);
theta = lp->bfp_prepareupdate(lp, rownr, colnr, pcol);
/* Verify numeric accuracy of the basis factorization and change to
the "theoretically" correct version of the theta */
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))) { /* This is my kludge - KE */
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
/* Force reinvertion and try another row if we did not find a bound-violated leaving column */
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
/* We may be infeasible, have lost dual feasibility, or simply have no valid entering
variable for the selected row. The strategy is to refactorize if we suspect numerical
problems and loss of dual feasibility; this is done if it has been a while since
refactorization. If not, first try to select a different row/leaving variable to
see if a valid entering variable can be found. Otherwise, determine this model
as infeasible. */
else {
/* As a first option, try to recover from any numerical trouble by refactorizing */
if(!refactRecent(lp)) {
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
#ifdef dual_UseRejectionList
/* Check for pivot size issues */
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
/* Check if we may have lost dual feasibility if we also did phase 1 here */
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;
}
/* Otherwise just determine that we are infeasible */
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;
}
}
}
/* Make sure that we enter the primal simplex with a high quality solution */
else if(inP1extra && !refactRecent(lp) && is_action(lp->improve, IMPROVE_INVERSE)) {
set_action(&lp->spx_action, ACTION_REINVERT);
bfpfinal = TRUE;
}
/* High quality solution with no leaving candidates available ... */
else {
bfpfinal = TRUE;
#ifdef dual_RemoveBasicFixedVars
/* See if we should try to eliminate basic fixed variables;
can be time-consuming for some models */
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
/* Check if we are INFEASIBLE for the case that the dual is used
as phase 1 before the primal simplex phase 2 */
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;
}
/* Check if we are FEASIBLE (and possibly also optimal) for the case that the
dual is used as phase 1 before the primal simplex phase 2 */
else if(inP1extra) {
/* Set default action; force an update of the rhs vector, adjusted for
the new P1extraVal=0 (set here so that usermessage() behaves properly) */
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
/* Optionally try another dual loop, if so selected by the user */
if((lp->simplex_strategy & SIMPLEX_DUAL_PRIMAL) && (lp->fixedvars == 0))
lp->spx_status = SWITCH_TO_PRIMAL;
#endif
changedphase = TRUE;
}
/* We are primal feasible and also optimal if we were in phase 2 */
else {
lp->simplex_mode = SIMPLEX_Phase2_DUAL;
/* Check if we still have equality slacks stuck in the basis; drive them out? */
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);
}
/* Check if we are still dual feasible; the default assumes that this check
is not necessary after the relaxed problem has been solved satisfactorily. */
colnr = 0;
if((dualoffset != 0) || (lp->bb_level <= 1) || (lp->improve & IMPROVE_BBSIMPLEX) || (lp->bb_rule & NODE_RCOSTFIXING)) { /* NODE_RCOSTFIXING fix */
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));
}
/* Determine if we are ready to break out of the loop */
if(!changedphase)
break;
}
/* Check if we are allowed to iterate on the chosen column and row */
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) (/*(candidatecount > 1) && */
(stallaccept != AUTOMATIC)),
prow, nzprow,
pcol, NULL, boundswaps);
/* Check if we should abandon iterations on finding that there is no
hope that this branch can improve on the incumbent B&B solution */
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++;
/* Update reduced costs for (pure) dual long-step phase 1 */
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
/* Do a fast update of reduced costs in preparation for the next iteration */
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--;
}
/* Refactorize if required to */
if(lp->bfp_mustrefactorize(lp)) {
if(invert(lp, INITSOL_USEZERO, bfpfinal)) {
#if 0
/* Verify dual feasibility in case we are attempting the extra dual loop */
if(changedphase && (dualoffset != 0) && !inP1extra && (lp->spx_status != SWITCH_TO_PRIMAL)) {
#if 1
if(!isDualFeasible(lp, lp->epsdual, &colnr, NULL, NULL)) {
#else
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(colnr > 0) {
#endif
lp->spx_status = SWITCH_TO_PRIMAL;
colnr = 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;
/* Compute the number of fixed basic and bounded variables (used in long duals) */
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
/* Reinvert for initialization, if necessary */
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))
/* if(isbb && (lp->bb_basis->pivots == 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);
/* Optionally try to do bound flips to obtain dual feasibility */
if(is_action(lp->improve, IMPROVE_DUALFEAS) || (lp->rows == 0))
boundflip_count = &i;
else
boundflip_count = NULL;
/* Loop for as long as is needed */
while(lp->spx_status == RUNNING) {
/* Check for dual and primal feasibility */
dualfeasible = isbb ||
isDualFeasible(lp, lp->epsprimal, boundflip_count, infeasibles, &dualoffset);
/* Recompute if the dual feasibility check included bound flips */
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");
}
/* Now do the simplex magic */
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);
}
/* Check for simplex outcomes that always involve breaking out of the loop;
this includes optimality, unboundedness, pure infeasibility (i.e. not
loss of feasibility), numerical failure and perturbation-based degeneracy
handling */
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))) {
/* Check if we should not loop here, but do perturbations */
if((lp->bb_level <= 1) || is_anti_degen(lp, ANTIDEGEN_DURINGBB))
break;
/* Assume that accuracy during B&B is high and that infeasibility is "real" */
#ifdef AssumeHighAccuracyInBB
if((lp->bb_level > 1) && (i == INFEASIBLE))
break;
#endif
}
/* Check for outcomes that may involve trying another simplex loop */
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;
/* Singular pivots are simply skipped by the inversion, leaving a row's
slack variable in the basis instead of the singular user variable. */
}
else {
lost_feas_state = (MYBOOL) (lp->spx_status == LOSTFEAS);
#if 0
/* Optionally handle loss of numerical accuracy as loss of feasibility,
but only attempt a single loop to try to recover from this. */
lost_feas_state |= (MYBOOL) ((lp->spx_status == NUMFAILURE) && (lost_feas_count < 1));
#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;
}
}
}
}
/* Update iteration tallies before returning */
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);
} /* spx_run */
lprec *make_lag(lprec *lpserver)
{
int i;
lprec *hlp;
MYBOOL ret;
REAL *duals;
/* Create a Lagrangean solver instance */
hlp = make_lp(0, lpserver->columns);
if(hlp != NULL) {
/* First create and core variable data */
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));
}
}
/* Then fill data for the Lagrangean constraints */
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)
/* Initialize / bound a MIP problem */
{
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)) {
/* 1. Copy the problem into a new relaxed instance, extracting Lagrangean constraints */
hlp = make_lag(lp);
/* 2. Run the Lagrangean relaxation */
status = solve(hlp);
/* 3. Copy the key results (bound) into the original problem */
lp->bb_heuristicOF = hlp->best_solution[0];
/* 4. Delete the helper heuristic */
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;
/* Make sure we have something to work with */
if(lp->spx_status != OPTIMAL) {
lp->lag_status = NOTRUN;
return( lp->lag_status );
}
/* Allocate iteration arrays */
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;
/* Prepare for Lagrangean iterations using results from relaxed problem */
oldpresolve = lp->do_presolve;
lp->do_presolve = PRESOLVE_NONE;
push_basis(lp, NULL, NULL, NULL);
/* Initialize variables (assume minimization problem in overall structure) */
Zlb = lp->best_solution[0];
Zub = start_bound;
Zbest = Zub;
Znow = Zlb;
Zprev = lp->infinite;
rhsmod = 0;
Phi = DEF_LAGCONTRACT; /* In the range 0-2.0 to guarantee convergence */
/* Phi = 0.15; */
LagFeas = FALSE;
Converged= FALSE;
AnyFeas = FALSE;
citer = 0;
nochange = 0;
/* Initialize reference and solution vectors; don't bother about the
original OF offset since we are maintaining an offset locally. */
/* #define DirectOverrideOF */
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;
/* Iterate to convergence, failure or user-specified termination */
while((lp->lag_status == RUNNING) && (citer < num_iter)) {
citer++;
/* Compute constraint feasibility gaps and associated sum of squares,
and determine feasibility over the Lagrangean constraints;
SubGrad is the subgradient, which here is identical to the slack. */
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;
}
/* Test for convergence and update */
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;
/* Modify step parameters and initialize ahead of next iteration */
Znow = lp->best_solution[0] - rhsmod;
if(Znow > Zub) {
/* Handle exceptional case where we overshoot */
Phi *= DEF_LAGCONTRACT;
StepSize *= (Zub-Zlb) / (Znow-Zlb);
}
else
#define LagBasisContract
#ifdef LagBasisContract
/* StepSize = Phi * (Zub - Znow) / SqrsumSubGrad; */
StepSize = Phi * (2-DEF_LAGCONTRACT) * (Zub - Znow) / SqrsumSubGrad;
#else
StepSize = Phi * (Zub - Znow) / SqrsumSubGrad;
#endif
/* Compute the new dual price vector (Lagrangean multipliers, lambda) */
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)) {
/* Handle case where we overshoot and need to correct (see above) */
if(Znow < Zub)
lp->lambda[i] = 0;
}
}
/* normalizeVector(lp->lambda, get_Lrows(lp)); */
/* Save the current vector if it is better */
if(LagFeas && (Znow < Zbest)) {
/* Recompute the objective function value in terms of the original values */
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);
/* Reset variables */
Zbest = Znow;
AnyFeas = TRUE;
nochange = 0;
}
else if(Znow == Zprev) {
nochange++;
if(nochange > LAG_SINGULARLIMIT) {
Phi *= 0.5;
nochange = 0;
}
}
Zprev = Znow;
/* Recompute the objective function values for the next iteration */
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
}
/* Recompute the fixed part of the new objective function */
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];
/* Print trace/debugging information, if specified */
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);
}
/* Solve the Lagrangean relaxation, handle failures and compute
the Lagrangean objective value, if successful */
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;
}
/* Compare optimal bases and contract if we have basis stationarity */
#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
}
/* Transfer solution values */
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]);
}
/* Do standard postprocessing */
Leave:
/* Set status variables and report */
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]);
}
/* Restore the original objective function */
#ifdef DirectOverrideOF
set_OF_override(lp, NULL);
#else
for(i = 1; i <= lp->columns; i++)
set_mat(lp, 0, i, OrigObj[i]);
#endif
/* ... and then free memory */
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)) {
/* Do standard initializations */
lp->solutioncount = 0;
lp->real_solution = lp->infinite;
set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT);
lp->bb_break = FALSE;
/* Do the call to the real underlying solver (note that
run_BB is replaceable with any compatible MIP solver) */
status = run_BB(lp);
/* Restore modified problem */
if(iprocessed)
postprocess(lp);
/* Restore data related to presolve (mainly a placeholder as of v5.1) */
Reconstruct:
if(!postsolve(lp, status))
report(lp, SEVERE, "spx_solve: Failure during postsolve.\n");
goto Leave;
}
/* If we get here, mat_validate(lp) failed. */
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 );
} /* spx_solve */
int lin_solve(lprec *lp)
{
int status = NOTRUN;
/* Don't do anything in case of an empty model */
lp->lag_status = NOTRUN;
/* if(get_nonzeros(lp) == 0) { */
if(lp->columns == 0) {
default_basis(lp);
lp->spx_status = NOTRUN;
return( /* OPTIMAL */ lp->spx_status);
}
/* Otherwise reset selected arrays before solving */
unset_OF_p1extra(lp);
free_duals(lp);
FREE(lp->drow);
FREE(lp->nzdrow);
if(lp->bb_cuttype != NULL)
freecuts_BB(lp);
/* Reset/initialize timers */
lp->timestart = timeNow();
lp->timeheuristic = 0;
lp->timepresolved = 0;
lp->timeend = 0;
/* Do heuristics ahead of solving the model */
if(heuristics(lp, AUTOMATIC) != RUNNING)
return( INFEASIBLE );
/* Solve the full, prepared model */
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");
}
/* Reset heuristic in preparation for next run (if any) */
lp->bb_heuristicOF = my_chsign(is_maxim(lp), lp->infinite);
/* Check that correct status code is returned */
/*
peno 26.12.07
status was not set to SUBOPTIMAL, only lp->spx_status
Bug occured by a change in 5.5.0.10 when && (lp->bb_totalnodes > 0) was added
added status =
See UnitTest3
*/
/*
peno 12.01.08
If an integer solution is found with the same objective value as the relaxed solution then
searching is stopped. This by setting lp->bb_break. However this resulted in a report of SUBOPTIMAL
solution. For this, && !bb_better(lp, OF_DUALLIMIT, OF_TEST_BE) is added in the test.
See UnitTest20
*/
if((lp->spx_status == OPTIMAL) && (lp->bb_totalnodes > 0)) {
if((lp->bb_break && !bb_better(lp, OF_DUALLIMIT, OF_TEST_BE)) /* ||
ISMASKSET(lp->trace, TRACE_NOBBLIMIT) */)
status = lp->spx_status = SUBOPTIMAL;
}
return( status );
}