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.
2889 Commits
5 Branches
0 Tags
31 MiB
 
 
 
 
 
 
Tree: 432d017bbf
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '432d017bbf'
${ noResults }
prism-accumulation/prism/ext/lpsolve55/src/lp_solve_5.5/lp_price.c

2107 lines
67 KiB
Raw Blame History

#include <string.h>
#include "commonlib.h"
#include "lp_lib.h"
#include "lp_report.h"
#include "lp_pricePSE.h"
#include "lp_price.h"
#if libBLAS > 0
#include "myblas.h"
#endif
#ifdef FORTIFY
# include "lp_fortify.h"
#endif
/* Simplex pricing utility module - w/interface for lp_solve v5.0+
-------------------------------------------------------------------------
Author: Kjell Eikland
Contact: kjell.eikland@broadpark.no
License terms: LGPL.
Requires: lp_lib.h, commonlib.h
Release notes:
v1.0.0 1 July 2004 Routines extracted from lp_lib.
v1.0.1 10 July 2004 Added comparison operators for determination
of entering and leaving variables.
Added routines for multiple and partial
pricing and made corresponding changes to
colprim and rowdual.
v1.0.2 20 August 2004 Implemented relative pivot size control in
rowprim and rowdual.
v1.1.0 15 October 2004 Added dual long step logic.
v1.1.1 22 October 2004 Added bound sort order to variable selections.
v1.2.0 24 March 2005 Completed multiple pricing logic.
------------------------------------------------------------------------- */
/* Comparison operators for entering and leaving variables for both the primal and
dual simplexes. The functions compare a candidate variable with an incumbent. */
int CMP_CALLMODEL compareImprovementVar(const pricerec *current, const pricerec *candidate)
{
register int result = COMP_PREFERNONE;
register lprec *lp = current->lp;
register REAL testvalue, margin = PREC_IMPROVEGAP;
int currentcolno, currentvarno = current->varno,
candidatecolno, candidatevarno = candidate->varno;
MYBOOL isdual = candidate->isdual;
if(isdual) {
candidatevarno = lp->var_basic[candidatevarno];
currentvarno = lp->var_basic[currentvarno];
}
candidatecolno = candidatevarno - lp->rows;
currentcolno = currentvarno - lp->rows;
/* Do pivot-based selection unless Bland's (first index) rule is active */
if(lp->_piv_rule_ != PRICER_FIRSTINDEX) {
MYBOOL candbetter;
/* Find the largest value - normalize in case of the dual, since
constraint violation is expressed as a negative number. */
/* Use absolute test for "small numbers", relative otherwise */
testvalue = candidate->pivot;
if(fabs(testvalue) < LIMIT_ABS_REL)
testvalue -= current->pivot;
else
testvalue = my_reldiff(testvalue, current->pivot);
if(isdual)
testvalue = -testvalue;
candbetter = (MYBOOL) (testvalue > 0);
if(candbetter) {
if(testvalue > margin)
result = COMP_PREFERCANDIDATE;
}
#if 0 /* Give more opportunity to optimize on non-primary criteria */
else if (testvalue < -margin)
#else /* Give reduced opportunity to optimize on non-primary criteria */
else if (testvalue < -lp->epsvalue)
#endif
result = COMP_PREFERINCUMBENT;
#ifdef UseSortOnBound
/* Extra selection criterion based on the variable's range;
variable with - DUAL: small bound out; PRIMAL: large bound in */
if(result == COMP_PREFERNONE) {
testvalue = lp->upbo[candidatevarno] - lp->upbo[currentvarno];
if(testvalue < -margin)
result = COMP_PREFERINCUMBENT;
else if(testvalue > margin)
result = COMP_PREFERCANDIDATE;
result = my_chsign(isdual, result);
}
#endif
#ifdef UseSortOnColumnLength
/* Prevent long columns from entering the basis */
if(result == COMP_PREFERNONE) {
if(candidatecolno > 0)
testvalue = mat_collength(lp->matA, candidatecolno) +
(is_obj_in_basis(lp) && (lp->obj[candidatecolno] != 0) ? 1 : 0);
else
testvalue = 1;
if(currentcolno > 0)
testvalue -= mat_collength(lp->matA, currentcolno) +
(is_obj_in_basis(lp) && (lp->obj[currentcolno] != 0) ? 1 : 0);
else
testvalue -= 1;
if(testvalue > 0)
result = COMP_PREFERINCUMBENT;
else if(testvalue < 0)
result = COMP_PREFERCANDIDATE;
result = my_chsign(isdual, result);
}
#endif
/* Select absolute best if the non-primary criteria failed to separate */
if((result == COMP_PREFERNONE) && candbetter) {
result = COMP_PREFERCANDIDATE;
goto Finish;
}
}
/* Final tie-breakers */
if(result == COMP_PREFERNONE) {
/* Add randomization tie-braker */
if(lp->piv_strategy & PRICE_RANDOMIZE) {
result = my_sign(PRICER_RANDFACT - rand_uniform(lp, 1.0));
if(candidatevarno < currentvarno)
result = -result;
}
/* Resolve ties via index ordinal */
if(result == COMP_PREFERNONE) {
if(candidatevarno < currentvarno)
result = COMP_PREFERCANDIDATE;
else /* if(candidatevarno > currentvarno) */
result = COMP_PREFERINCUMBENT;
if(lp->_piv_left_)
result = -result;
}
}
Finish:
return( result );
}
int CMP_CALLMODEL compareSubstitutionVar(const pricerec *current, const pricerec *candidate)
{
register int result = COMP_PREFERNONE;
register lprec *lp = current->lp;
register REAL testvalue = candidate->theta,
margin = current->theta;
MYBOOL isdual = candidate->isdual, candbetter;
int currentcolno, currentvarno = current->varno,
candidatecolno, candidatevarno = candidate->varno;
if(!isdual) {
candidatevarno = lp->var_basic[candidatevarno];
currentvarno = lp->var_basic[currentvarno];
}
candidatecolno = candidatevarno - lp->rows;
currentcolno = currentvarno - lp->rows;
/* Compute the ranking test metric. */
if(isdual) {
testvalue = fabs(testvalue);
margin = fabs(margin);
}
/* Use absolute test for "small numbers", relative otherwise */
if(fabs(testvalue) < LIMIT_ABS_REL)
testvalue -= margin;
else
testvalue = my_reldiff(testvalue, margin);
/* Find if the new Theta is smaller or near equal (i.e. testvalue <= eps)
compared to the previous best; ties will be broken by pivot size or index
NB! The margin below is essential in maintaining primal/dual feasibility
during the primal/dual simplex, respectively. Sometimes a small
value prevents the selection of a suitable pivot, thereby weakening
the numerical stability of some models */
margin = PREC_SUBSTFEASGAP;
candbetter = (MYBOOL) (testvalue < 0);
if(candbetter) {
if(testvalue < -margin)
result = COMP_PREFERCANDIDATE;
}
else if(testvalue > margin)
result = COMP_PREFERINCUMBENT;
/* Resolve a tie */
if(result == COMP_PREFERNONE) {
REAL currentpivot = fabs(current->pivot),
candidatepivot = fabs(candidate->pivot);
/* Handle first index / Bland's rule specially */
if(lp->_piv_rule_ == PRICER_FIRSTINDEX) {
#if 1
/* Special secondary selection by pivot size (limited stability protection) */
margin = candidate->epspivot;
if((candidatepivot >= margin) && (currentpivot < margin))
result = COMP_PREFERCANDIDATE;
#endif
}
else {
/* General secondary selection based on pivot size */
#if 0
if(candidatepivot > MIN_STABLEPIVOT)
testvalue = my_reldiff(testvalue, currentpivot);
else
#endif
testvalue = candidatepivot - currentpivot;
if(testvalue > margin)
result = COMP_PREFERCANDIDATE;
else if(testvalue < -margin)
result = COMP_PREFERINCUMBENT;
#ifdef UseSortOnBound
/* Extra selection criterion based on the variable's range;
variable with - PRIMAL: small bound out; DUAL: large bound in */
if(result == COMP_PREFERNONE) {
testvalue = lp->upbo[candidatevarno] - lp->upbo[currentvarno];
if(testvalue < -margin)
result = COMP_PREFERCANDIDATE;
else if(testvalue > margin)
result = COMP_PREFERINCUMBENT;
result = my_chsign(isdual, result);
}
#endif
#ifdef UseSortOnColumnLength
/* Prevent long columns from entering the basis */
if(result == COMP_PREFERNONE) {
if(candidatecolno > 0)
testvalue = mat_collength(lp->matA, candidatecolno) +
(is_obj_in_basis(lp) && (lp->obj[candidatecolno] != 0) ? 1 : 0);
else
testvalue = 1;
if(currentcolno > 0)
testvalue -= mat_collength(lp->matA, currentcolno) +
(is_obj_in_basis(lp) && (lp->obj[currentcolno] != 0) ? 1 : 0);
else
testvalue -= 1;
if(testvalue > 0)
result = COMP_PREFERCANDIDATE;
else if(testvalue < 0)
result = COMP_PREFERINCUMBENT;
result = my_chsign(isdual, result);
}
#endif
}
}
/* Select absolute best if the non-primary criteria failed to separate */
if((result == COMP_PREFERNONE) && candbetter) {
result = COMP_PREFERCANDIDATE;
goto Finish;
}
/* Final tie-breakers */
if(result == COMP_PREFERNONE) {
/* Add randomization tie-braker */
if(lp->piv_strategy & PRICE_RANDOMIZE) {
result = my_sign(PRICER_RANDFACT - rand_uniform(lp, 1.0));
if(candidatevarno < currentvarno)
result = -result;
}
/* Resolve ties via index ordinal (also prefers slacks over user variables) */
if(result == COMP_PREFERNONE) {
if(candidatevarno < currentvarno)
result = COMP_PREFERCANDIDATE;
else /* if(candidatevarno > currentvarno) */
result = COMP_PREFERINCUMBENT;
if(lp->_piv_left_)
result = -result;
}
}
Finish:
return( result );
}
int CMP_CALLMODEL compareBoundFlipVar(const pricerec *current, const pricerec *candidate)
{
register REAL testvalue, margin;
register int result = COMP_PREFERNONE;
register lprec *lp = current->lp;
MYBOOL candbetter;
int currentvarno = current->varno,
candidatevarno = candidate->varno;
if(!current->isdual) {
candidatevarno = lp->var_basic[candidatevarno];
currentvarno = lp->var_basic[currentvarno];
}
/* Compute the ranking test metric. */
testvalue = candidate->theta;
margin = current->theta;
if(candidate->isdual) {
testvalue = fabs(testvalue);
margin = fabs(margin);
}
if(fabs(margin) < LIMIT_ABS_REL)
testvalue -= margin;
else
testvalue = my_reldiff(testvalue, margin);
/* Find if the new Theta is smaller or near equal (i.e. testvalue <= eps)
compared to the previous best; ties will be broken by pivot size or index */
margin = PREC_SUBSTFEASGAP;
candbetter = (MYBOOL) (testvalue < 0);
if(candbetter) {
if(testvalue < -margin)
result = COMP_PREFERCANDIDATE;
}
else if(testvalue > margin)
result = COMP_PREFERINCUMBENT;
/* Resolve a tie */
if(result == COMP_PREFERNONE) {
/* Tertiary selection based on priority for large pivot sizes */
if(result == COMP_PREFERNONE) {
REAL currentpivot = fabs(current->pivot),
candidatepivot = fabs(candidate->pivot);
if(candidatepivot > currentpivot+margin)
result = COMP_PREFERCANDIDATE;
else if(candidatepivot < currentpivot-margin)
result = COMP_PREFERINCUMBENT;
}
/* Secondary selection based on priority for narrow-bounded variables */
if(result == COMP_PREFERNONE)
result = compareREAL(&(lp->upbo[currentvarno]),
&(lp->upbo[candidatevarno]));
}
/* Select absolute best if the non-primary criteria failed to separate */
if((result == COMP_PREFERNONE) && candbetter) {
result = COMP_PREFERCANDIDATE;
goto Finish;
}
/* Quaternary selection by index value */
if(result == COMP_PREFERNONE) {
if(candidatevarno < currentvarno)
result = COMP_PREFERCANDIDATE;
else
result = COMP_PREFERINCUMBENT;
if(lp->_piv_left_)
result = -result;
}
Finish:
return( result );
}
/* Validity operators for entering and leaving columns for both the primal and dual
simplex. All candidates must satisfy these tests to qualify to be allowed to be
a subject for the comparison functions/operators. */
STATIC MYBOOL validImprovementVar(pricerec *candidate)
{
register REAL candidatepivot = fabs(candidate->pivot);
#ifdef Paranoia
return( (MYBOOL) ((candidate->varno > 0) && (candidatepivot > candidate->lp->epsvalue)) );
#else
return( (MYBOOL) (candidatepivot > candidate->lp->epsvalue) );
#endif
}
STATIC MYBOOL validSubstitutionVar(pricerec *candidate)
{
register lprec *lp = candidate->lp;
register REAL theta = (candidate->isdual ? fabs(candidate->theta) : candidate->theta);
#ifdef Paranoia
if(candidate->varno <= 0)
return( FALSE );
else
#endif
if(fabs(candidate->pivot) >= lp->infinite)
return( (MYBOOL) (theta < lp->infinite) );
else
return( (MYBOOL) ((theta < lp->infinite) &&
(fabs(candidate->pivot) >= candidate->epspivot)) );
}
int CMP_CALLMODEL compareImprovementQS(const UNIONTYPE QSORTrec *current, const UNIONTYPE QSORTrec *candidate)
{
return( compareImprovementVar((pricerec *) current->pvoidint2.ptr, (pricerec *) candidate->pvoidint2.ptr) );
}
int CMP_CALLMODEL compareSubstitutionQS(const UNIONTYPE QSORTrec *current, const UNIONTYPE QSORTrec *candidate)
{
return( compareBoundFlipVar((pricerec *) current->pvoidint2.ptr, (pricerec *) candidate->pvoidint2.ptr) );
/* return( compareSubstitutionVar((pricerec *) current->self, (pricerec *) candidate->self) ); */
}
/* Function to add a valid pivot candidate into the specified list */
STATIC int addCandidateVar(pricerec *candidate, multirec *multi, findCompare_func findCompare, MYBOOL allowSortedExpand)
{
int insertpos, delta = 1;
pricerec *targetrec;
/* Find the insertion point (if any) */
if((multi->freeList[0] == 0) ||
(multi->sorted && allowSortedExpand) ||
(candidate->isdual && (multi->used == 1) && ((multi->step_last >= multi->epszero) ||
multi_truncatingvar(multi, ((pricerec *) (multi->sortedList[0].pvoidreal.ptr))->varno)))
) {
UNIONTYPE QSORTrec searchTarget;
/* Make sure that the list is sorted before the search for an insertion point */
if((multi->freeList[0] == 0) && !multi->sorted) {
multi->sorted = QS_execute(multi->sortedList, multi->used, findCompare, &insertpos);
multi->dirty = (MYBOOL) (insertpos > 0);
}
/* Perform the search */
searchTarget.pvoidint2.ptr = (void *) candidate;
insertpos = sizeof(searchTarget);
insertpos = findIndexEx(&searchTarget, multi->sortedList-delta, multi->used, delta, insertpos, findCompare, TRUE);
if(insertpos > 0)
return( -1 );
insertpos = -insertpos - delta;
/* Check if the candidate is worse than the worst of the list */
if(((insertpos >= multi->size) && (multi->freeList[0] == 0)) ||
((insertpos == multi->used) && (!allowSortedExpand ||
(multi->step_last >= multi->epszero))))
return( -1 );
#ifdef Paranoia
/* Do validation */
if((insertpos < 0) || (insertpos > multi->used))
return( -1 );
#endif
/* Define the target for storing the candidate;
Case 1: List is full and we must discard the previously worst candidate
Case 2: List is not full and we simply use the next free position */
if(multi->freeList[0] == 0)
targetrec = (pricerec *) multi->sortedList[multi->used-1].pvoidreal.ptr;
else {
delta = multi->freeList[0]--;
delta = multi->freeList[delta];
targetrec = &(multi->items[delta]);
}
}
else {
delta = multi->freeList[0]--;
delta = multi->freeList[delta];
targetrec = &(multi->items[delta]);
insertpos = multi->used;
}
/* Insert the new candidate record in the data store */
MEMCOPY(targetrec, candidate, 1);
/* Store the pointer data and handle tree cases:
Case 1: The list is unsorted and not full; simply append pointer to list,
Case 2: The list is sorted and full; insert the pointer by discarding previous last,
Case 3: The list is sorted and not full; shift the inferior items down, and increment count */
if((multi->used < multi->size) && (insertpos >= multi->used)) {
QS_append(multi->sortedList, insertpos, targetrec);
multi->used++;
}
else {
if(multi->used == multi->size)
QS_insert(multi->sortedList, insertpos, targetrec, multi->size-1); /* Discard previous last */
else {
QS_insert(multi->sortedList, insertpos, targetrec, multi->used); /* Keep previous last */
multi->used++;
}
}
multi->active = insertpos;
#ifdef Paranoia
if((insertpos >= multi->size) || (insertpos >= multi->used))
report(multi->lp, SEVERE, "addCandidateVar: Insertion point beyond limit!\n");
#endif
return( insertpos );
}
STATIC MYBOOL findImprovementVar(pricerec *current, pricerec *candidate, MYBOOL collectMP, int *candidatecount)
/* PRIMAL: Find a variable to enter the basis
DUAL: Find a variable to leave the basis
Allowed variable set: Any pivot PRIMAL:larger or DUAL:smaller than threshold value of 0 */
{
MYBOOL Action = FALSE,
#ifdef ExtractedValidityTest
Accept = TRUE;
#else /* Check for validity and compare result with previous best */
Accept = validImprovementVar(candidate);
#endif
if(Accept) {
if(candidatecount != NULL)
(*candidatecount)++;
if(collectMP) {
if(addCandidateVar(candidate, current->lp->multivars, (findCompare_func *) compareImprovementQS, FALSE) < 0)
return(Action);
}
if(current->varno > 0)
Accept = (MYBOOL) (compareImprovementVar(current, candidate) > 0);
}
/* Apply candidate if accepted */
if(Accept) {
(*current) = *candidate;
/* Force immediate acceptance for Bland's rule using the primal simplex */
if(!candidate->isdual)
Action = (MYBOOL) (candidate->lp->_piv_rule_ == PRICER_FIRSTINDEX);
}
return(Action);
}
/* Bound flip variable accumulation routine */
STATIC MYBOOL collectMinorVar(pricerec *candidate, multirec *longsteps, MYBOOL isphase2, MYBOOL isbatch)
{
int inspos;
/* 1. Check for ratio and pivot validity (to have the extra flexibility that all
bound-flip candidates are also possible as basis-entering variables */
if(!validSubstitutionVar(candidate))
return( FALSE );
/* 2. If the free-list is empty we need to see if we have a better candidate,
and for this the candidate list has to be sorted by merit */
if(!isbatch &&
!longsteps->sorted && (longsteps->used > 1) &&
((longsteps->freeList[0] == 0) ||
multi_truncatingvar(longsteps, candidate->varno) ||
(longsteps->step_last >= longsteps->epszero) )) {
longsteps->sorted = QS_execute(longsteps->sortedList, longsteps->used,
(findCompare_func *) compareSubstitutionQS, &inspos);
longsteps->dirty = (MYBOOL) (inspos > 0);
if(longsteps->dirty)
multi_recompute(longsteps, 0, isphase2, TRUE);
}
/* 3. Now handle three cases...
- Add to the list when the list is not full and there is opportunity for improvement,
- Check if we should replace an incumbent when the list is full,
- Check if we should replace an incumbent when the list is not full, there is no room
for improvement, but the current candidate is better than an incumbent. */
inspos = addCandidateVar(candidate, longsteps, (findCompare_func *) compareSubstitutionQS, TRUE);
/* 4. Recompute steps and objective, and (if relevant) determine if we
may be suboptimal in relation to an incumbent MILP solution. */
return( (MYBOOL) (inspos >= 0) &&
((isbatch == TRUE) || multi_recompute(longsteps, inspos, isphase2, TRUE)) );
}
STATIC MYBOOL findSubstitutionVar(pricerec *current, pricerec *candidate, int *candidatecount)
/* PRIMAL: Find a variable to leave the basis
DUAL: Find a variable to enter the basis
Allowed variable set: "Equal-valued" smallest thetas! */
{
MYBOOL Action = FALSE,
#ifdef ExtractedValidityTest
Accept = TRUE;
#else /* Check for validity and comparison result with previous best */
Accept = validSubstitutionVar(candidate);
#endif
if(Accept) {
if(candidatecount != NULL)
(*candidatecount)++;
if(current->varno != 0)
Accept = (MYBOOL) (compareSubstitutionVar(current, candidate) > 0);
}
/* Apply candidate if accepted */
if(Accept) {
(*current) = *candidate;
/* Force immediate acceptance for Bland's rule using the dual simplex */
#ifdef ForceEarlyBlandRule
if(candidate->isdual)
Action = (MYBOOL) (candidate->lp->_piv_rule_ == PRICER_FIRSTINDEX);
#endif
}
return(Action);
}
/* Partial pricing management routines */
STATIC partialrec *partial_createBlocks(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = (partialrec *) calloc(1, sizeof(*blockdata));
blockdata->lp = lp;
blockdata->blockcount = 1;
blockdata->blocknow = 1;
blockdata->isrow = isrow;
return(blockdata);
}
STATIC int partial_countBlocks(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( 1 );
else
return( blockdata->blockcount );
}
STATIC int partial_activeBlocks(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( 1 );
else
return( blockdata->blocknow );
}
STATIC void partial_freeBlocks(partialrec **blockdata)
{
if((blockdata == NULL) || (*blockdata == NULL))
return;
FREE((*blockdata)->blockend);
FREE((*blockdata)->blockpos);
FREE(*blockdata);
}
/* Function to provide for left-right or right-left scanning of entering/leaving
variables; note that *end must have been initialized by the calling routine! */
STATIC void makePriceLoop(lprec *lp, int *start, int *end, int *delta)
{
int offset = is_piv_mode(lp, PRICE_LOOPLEFT);
if((offset) ||
(((lp->total_iter+offset) % 2 == 0) && is_piv_mode(lp, PRICE_LOOPALTERNATE))) {
*delta = -1; /* Step backwards - "left" */
swapINT(start, end);
lp->_piv_left_ = TRUE;
}
else {
*delta = 1; /* Step forwards - "right" */
lp->_piv_left_ = FALSE;
}
}
/* Routine to verify accuracy of the current basis factorization */
STATIC MYBOOL serious_facterror(lprec *lp, REAL *bvector, int maxcols, REAL tolerance)
{
int i, j, ib, ie, nz, nc;
REAL sum, tsum = 0, err = 0;
MATrec *mat = lp->matA;
if(bvector == 0)
bvector = lp->bsolveVal;
nc =0;
nz = 0;
for(i = 1; (i <= lp->rows) && (nc <= maxcols); i++) {
/* Do we have a non-slack variable? (we choose to skip slacks,
since they have "natural" good accuracy properties) */
j = lp->var_basic[i] - lp->rows;
if(j <= 0)
continue;
nc++;
/* Compute cross product for basic, non-slack column */
ib = mat->col_end[j-1];
ie = mat->col_end[j];
nz += ie - ib;
sum = get_OF_active(lp, j+lp->rows, bvector[0]);
for(; ib < ie; ib++)
sum += COL_MAT_VALUE(ib)*bvector[COL_MAT_ROWNR(ib)];
/* Catch high precision early, so we don't to uneccessary work */
tsum += sum;
SETMAX(err, fabs(sum));
if((tsum / nc > tolerance / 100) && (err < tolerance / 100))
break;
}
err /= mat->infnorm;
return( (MYBOOL) (err >= tolerance) );
}
/* Computation of reduced costs */
STATIC void update_reducedcosts(lprec *lp, MYBOOL isdual, int leave_nr, int enter_nr, REAL *prow, REAL *drow)
{
/* "Fast" update of the dual reduced cost vector; note that it must be called
after the pivot operation and only applies to a major "true" iteration */
int i;
REAL hold;
if(isdual) {
hold = -drow[enter_nr]/prow[enter_nr];
for(i=1; i <= lp->sum; i++)
if(!lp->is_basic[i]) {
if(i == leave_nr)
drow[i] = hold;
else {
drow[i] += hold*prow[i];
my_roundzero(drow[i], lp->epsmachine);
}
}
}
else
report(lp, SEVERE, "update_reducedcosts: Cannot update primal reduced costs!\n");
}
STATIC void compute_reducedcosts(lprec *lp, MYBOOL isdual, int row_nr, int *coltarget, MYBOOL dosolve,
REAL *prow, int *nzprow,
REAL *drow, int *nzdrow,
int roundmode)
{
REAL epsvalue = lp->epsvalue; /* Any larger value can produce a suboptimal result */
roundmode |= MAT_ROUNDRC;
if(isdual) {
bsolve_xA2(lp, coltarget,
row_nr, prow, epsvalue, nzprow, /* Calculate net sensitivity given a leaving variable */
0, drow, epsvalue, nzdrow, /* Calculate the net objective function values */
roundmode);
}
else {
REAL *bVector;
#if 1 /* Legacy mode, that is possibly a little faster */
if((lp->multivars == NULL) && (lp->P1extraDim == 0))
bVector = drow;
else
#endif
bVector = lp->bsolveVal;
if(dosolve) {
bsolve(lp, 0, bVector, lp->bsolveIdx, epsvalue*DOUBLEROUND, 1.0);
if(!isdual && (row_nr == 0) && (lp->improve & IMPROVE_SOLUTION) && !refactRecent(lp) &&
serious_facterror(lp, bVector, lp->rows, lp->epsvalue))
set_action(&lp->spx_action, ACTION_REINVERT);
}
prod_xA(lp, coltarget,
bVector, lp->bsolveIdx, epsvalue, 1.0,
drow, nzdrow, roundmode);
}
}
/* Primal: Prevent acceptance of an entering variable when the magnitude of
other candidates is also very small.
Dual: Prevent acceptance of a leaving variable when the magnitude of
other candidates is also very small.
Both of these cases are associated with numerical stalling, which we could
argue should be detected and handled by the stalling monitor routine. */
STATIC MYBOOL verify_stability(lprec *lp, MYBOOL isprimal, REAL xfeas, REAL sfeas, int nfeas)
{
MYBOOL testOK = TRUE;
return( testOK );
#if 1
/* Try to make dual feasibility as tight as possible */
if(!isprimal)
/* if(lp->P1extraVal == 0) */
{
xfeas /= (1+lp->rhsmax);
sfeas /= (1+lp->rhsmax);
}
#endif
xfeas = fabs(xfeas); /* Maximum (positive) infeasibility */
/* if(xfeas < lp->epspivot) { */
if(xfeas < lp->epssolution) {
REAL f;
sfeas = fabs(sfeas); /* Make sum of infeasibilities positive */
xfeas = (sfeas-xfeas)/nfeas; /* Average "residual" feasibility */
f = 1 + log10((REAL) nfeas); /* Some numerical complexity scalar */
/* Numerical errors can interact to cause non-convergence, and the
idea is to relax the tolerance to account for this and only
marginally weakening the (user-specified) tolerance. */
if((sfeas-xfeas) < f*lp->epsprimal)
testOK = FALSE;
}
return( testOK );
}
/* Find an entering column for the case that the specified basic variable
is fixed or zero - typically used for artificial variable elimination */
STATIC int find_rowReplacement(lprec *lp, int rownr, REAL *prow, int *nzprow)
/* The logic in this section generally follows Chvatal: Linear Programming, p. 130
Basically, the function is a specialized coldual(). */
{
int i, bestindex;
REAL bestvalue;
/* Solve for "local reduced cost" */
set_action(&lp->piv_strategy, PRICE_FORCEFULL);
compute_reducedcosts(lp, TRUE, rownr, NULL, TRUE,
prow, nzprow, NULL, NULL, MAT_ROUNDDEFAULT);
clear_action(&lp->piv_strategy, PRICE_FORCEFULL);
/* Find a suitably non-singular variable to enter ("most orthogonal") */
bestindex = 0;
bestvalue = 0;
for(i = 1; i <= lp->sum-abs(lp->P1extraDim); i++) {
if(!lp->is_basic[i] && !is_fixedvar(lp, i) &&
(fabs(prow[i]) > bestvalue)) {
bestindex = i;
bestvalue = fabs(prow[i]);
}
}
/* Prepare to update inverse and pivot/iterate (compute Bw=a) */
if(i > lp->sum-abs(lp->P1extraDim))
bestindex = 0;
else
fsolve(lp, bestindex, prow, nzprow, lp->epsmachine, 1.0, TRUE);
return( bestindex );
}
/* Find the primal simplex entering non-basic column variable */
STATIC int colprim(lprec *lp, REAL *drow, int *nzdrow, MYBOOL skipupdate, int partialloop, int *candidatecount, MYBOOL updateinfeas, REAL *xviol)
{
int i, ix, iy, iz, ninfeas, nloop = 0;
REAL f, sinfeas, xinfeas, epsvalue = lp->epsdual;
pricerec current, candidate;
MYBOOL collectMP = FALSE;
int *coltarget = NULL;
/* Identify pivot column according to pricing strategy; set
entering variable initial threshold reduced cost value to "0" */
current.pivot = lp->epsprimal; /* Minimum acceptable improvement */
current.varno = 0;
current.lp = lp;
current.isdual = FALSE;
candidate.lp = lp;
candidate.isdual = FALSE;
*candidatecount = 0;
/* Update local value of pivot setting and determine active multiple pricing set */
lp->_piv_rule_ = get_piv_rule(lp);
doLoop:
nloop++;
if((lp->multivars != NULL) && ((lp->simplex_mode & SIMPLEX_PRIMAL_PRIMAL) != 0)) {
collectMP = multi_mustupdate(lp->multivars);
if(collectMP) {
multi_restart(lp->multivars);
coltarget = NULL;
}
else
coltarget = multi_indexSet(lp->multivars, FALSE);
}
/* Compute reduced costs c - c*Inv(B), if necessary
(i.e. the previous iteration was not a "minor" iteration/bound flip) */
if(!skipupdate) {
#ifdef UsePrimalReducedCostUpdate
/* Recompute from scratch only at the beginning, otherwise update */
if((lp->current_iter > 0) && (refactRecent(lp) == AUTOMATIC))
#endif
compute_reducedcosts(lp, FALSE, 0, coltarget, (MYBOOL) ((nloop <= 1) || (partialloop > 1)),
NULL, NULL,
drow, nzdrow,
MAT_ROUNDDEFAULT);
}
/* Loop over active partial column set; we presume that reduced costs
have only been updated for columns in the active partial range. */
ix = 1;
iy = nzdrow[0];
ninfeas = 0;
xinfeas = 0;
sinfeas = 0;
makePriceLoop(lp, &ix, &iy, &iz);
iy *= iz;
for(; ix*iz <= iy; ix += iz) {
i = nzdrow[ix];
#if 0 /* Not necessary since we masked them out in compute_reducedcosts() */
if(i > lp->sum-abs(lp->P1extraDim))
continue;
#endif
/* Check if the pivot candidate is on the block-list */
if(lp->rejectpivot[0] > 0) {
int kk;
for(kk = 1; (kk <= lp->rejectpivot[0]) && (i != lp->rejectpivot[kk]); kk++);
if(kk <= lp->rejectpivot[0])
continue;
}
/* Retrieve the applicable reduced cost - threshold should not be smaller than 0 */
f = my_chsign(lp->is_lower[i], drow[i]);
if(f <= epsvalue)
continue;
/* Find entering variable according to strategy (largest positive f) */
ninfeas++;
SETMAX(xinfeas, f);
sinfeas += f;
candidate.pivot = normalizeEdge(lp, i, f, FALSE);
candidate.varno = i;
if(findImprovementVar(&current, &candidate, collectMP, candidatecount))
break;
}
/* Check if we should loop again after a multiple pricing update */
if(lp->multivars != NULL) {
if(collectMP) {
if(!lp->multivars->sorted)
lp->multivars->sorted = QS_execute(lp->multivars->sortedList, lp->multivars->used,
(findCompare_func *) compareImprovementQS, NULL);
coltarget = multi_indexSet(lp->multivars, TRUE);
}
else if((current.varno == 0) && (lp->multivars->retries == 0)) {
ix = partial_blockStart(lp, FALSE);
iy = partial_blockEnd(lp, FALSE);
lp->multivars->used = 0;
lp->multivars->retries++;
goto doLoop;
}
/* Shrink the candidate list */
lp->multivars->retries = 0;
if(current.varno != 0)
multi_removevar(lp->multivars, current.varno);
}
/* Check for optimality */
if(xviol != NULL)
*xviol = xinfeas;
if(updateinfeas)
lp->suminfeas = fabs(sinfeas);
if((lp->multivars == NULL) && (current.varno > 0) &&
!verify_stability(lp, TRUE, xinfeas, sinfeas, ninfeas))
current.varno = 0;
/* Produce statistics */
if(lp->spx_trace) {
if(current.varno > 0)
report(lp, DETAILED, "colprim: Column %d reduced cost = " RESULTVALUEMASK "\n",
current.varno, current.pivot);
else
report(lp, DETAILED, "colprim: No positive reduced costs found, optimality!\n");
}
return( current.varno );
} /* colprim */
/* Find the primal simplex leaving basic column variable */
STATIC int rowprim(lprec *lp, int colnr, LREAL *theta, REAL *pcol, int *nzpcol, MYBOOL forceoutEQ, REAL *xviol)
{
int i, ii, iy, iz, Hpass, k, *nzlist;
LREAL f, savef;
REAL Heps, Htheta, Hlimit, epsvalue, epspivot, p = 0.0;
pricerec current, candidate;
MYBOOL isupper = !lp->is_lower[colnr], HarrisTwoPass = FALSE;
/* Update local value of pivot setting */
lp->_piv_rule_ = get_piv_rule(lp);
if(nzpcol == NULL)
nzlist = (int *) mempool_obtainVector(lp->workarrays, lp->rows+1, sizeof(*nzlist));
else
nzlist = nzpcol;
/* Find unconditional non-zeros and optionally compute relative size of epspivot */
epspivot = lp->epspivot;
epsvalue = lp->epsvalue;
Hlimit = 0;
Htheta = 0;
k = 0;
for(i = 1; i <= lp->rows; i++) {
p = fabs(pcol[i]);
if(p > Hlimit)
Hlimit = p;
if(p > epsvalue) {
k++;
nzlist[k] = i;
SETMAX(Htheta, p);
}
#ifdef Paranoia
else {
if(lp->spx_trace)
report(lp, FULL, "rowprim: Row %d with pivot " RESULTVALUEMASK " rejected as too small\n",
i, p);
}
#endif
}
if(xviol != NULL)
*xviol = Htheta;
Htheta = 0;
/* Update non-zero list based on the new pivot threshold */
#ifdef UseRelativePivot_Primal
/* epspivot *= sqrt(lp->matA->dynrange) / lp->matA->infnorm; */
epspivot /= MAX(1, sqrt(lp->matA->colmax[colnr]));
iy = k;
k = 0;
p = 0;
for(ii = 1; ii <= iy; ii++) {
i = nzlist[ii];
p = fabs(pcol[i]);
/* Compress the list of valid alternatives, if appropriate */
if(p > epspivot) {
k++;
nzlist[k] = i;
}
#ifdef Paranoia
else {
if(lp->spx_trace)
report(lp, FULL, "rowprim: Row %d with pivot " RESULTVALUEMASK " rejected as too small\n",
i, p);
}
#endif
}
#endif
/* Initialize counters */
nzlist[0] = k;
k = 0;
Retry:
k++;
HarrisTwoPass = is_piv_mode(lp, PRICE_HARRISTWOPASS);
if(HarrisTwoPass)
Hpass = 1;
else
Hpass = 2;
current.theta = lp->infinite;
current.pivot = 0;
current.varno = 0;
current.isdual = FALSE;
current.epspivot = epspivot;
current.lp = lp;
candidate.epspivot = epspivot;
candidate.isdual = FALSE;
candidate.lp = lp;
savef = 0;
for(; Hpass <= 2; Hpass++) {
Htheta = lp->infinite;
if(Hpass == 1) {
Hlimit = lp->infinite; /* Don't apply any limit in the first pass */
Heps = epspivot/lp->epsprimal; /* Scaled to lp->epsprimal used in compute_theta() */
}
else {
Hlimit = current.theta; /* This is the smallest Theta of the first pass */
Heps = 0.0;
}
current.theta = lp->infinite;
current.pivot = 0;
current.varno = 0;
savef = 0;
ii = 1;
iy = nzlist[0];
makePriceLoop(lp, &ii, &iy, &iz);
iy *= iz;
for(; ii*iz <= iy; ii += iz) {
i = nzlist[ii];
f = pcol[i];
candidate.theta = f;
candidate.pivot = f;
candidate.varno = i;
/*i =*/ compute_theta(lp, i, &candidate.theta, isupper,
my_if(lp->upbo[lp->var_basic[i]] < lp->epsprimal, Heps/10, Heps), TRUE);
if(fabs(candidate.theta) >= lp->infinite) {
savef = f;
candidate.theta = 2*lp->infinite;
continue;
}
/* Find the candidate leaving variable according to strategy (smallest theta) */
if((Hpass == 2) && (candidate.theta > Hlimit))
continue;
/* Give a slight preference to fixed variables (mainly equality slacks) */
if(forceoutEQ) {
p = candidate.pivot;
if(lp->upbo[lp->var_basic[i]] < lp->epsprimal) {
/* Give an extra early boost to equality slack elimination, if specified */
if(forceoutEQ == AUTOMATIC)
candidate.pivot *= 1.0+lp->epspivot;
else
candidate.pivot *= 10.0;
}
}
if(HarrisTwoPass) {
f = candidate.theta;
if(Hpass == 2)
candidate.theta = 1;
if(findSubstitutionVar(&current, &candidate, NULL))
break;
if((Hpass == 2) && (current.varno == candidate.varno))
Htheta = f;
}
else
if(findSubstitutionVar(&current, &candidate, NULL))
break;
/* Restore temporarily modified pivot */
if(forceoutEQ && (current.varno == candidate.varno))
current.pivot = p;
}
}
if(HarrisTwoPass)
current.theta = Htheta;
/* Handle case of no available leaving variable */
if(current.varno == 0) {
if(lp->upbo[colnr] >= lp->infinite) {
/* Optionally try again with reduced pivot threshold level */
if(k < 2) {
epspivot = epspivot / 10;
goto Retry;
}
}
else {
#if 1
i = 1;
while((pcol[i] >= 0) && (i <= lp->rows))
i++;
if(i > lp->rows) { /* Empty column with upper bound! */
lp->is_lower[colnr] = !lp->is_lower[colnr];
/* lp->is_lower[colnr] = FALSE; */
lp->rhs[0] += lp->upbo[colnr]*pcol[0];
}
else /* if(pcol[i]<0) */
{
current.varno = i;
}
#endif
}
}
else if(current.theta >= lp->infinite) {
report(lp, IMPORTANT, "rowprim: Numeric instability pcol[%d] = %g, rhs[%d] = %g, upbo = %g\n",
current.varno, savef, current.varno, lp->rhs[current.varno],
lp->upbo[lp->var_basic[current.varno]]);
}
/* Return working array to pool */
if(nzpcol == NULL)
mempool_releaseVector(lp->workarrays, (char *) nzlist, FALSE);
if(lp->spx_trace)
report(lp, DETAILED, "row_prim: %d, pivot size = " RESULTVALUEMASK "\n",
current.varno, current.pivot);
/* *theta = current.theta; */
*theta = fabs(current.theta);
return(current.varno);
} /* rowprim */
/* Find the dual simplex leaving basic variable */
STATIC int rowdual(lprec *lp, REAL *rhvec, MYBOOL forceoutEQ, MYBOOL updateinfeas, REAL *xviol)
{
int k, i, iy, iz, ii, ninfeas;
register REAL rh;
REAL up, lo = 0,
epsvalue, sinfeas, xinfeas;
pricerec current, candidate;
MYBOOL collectMP = FALSE;
/* Initialize */
if(rhvec == NULL)
rhvec = lp->rhs;
epsvalue = lp->epsdual;
current.pivot = -epsvalue; /* Initialize leaving variable threshold; "less than 0" */
current.theta = 0;
current.varno = 0;
current.isdual = TRUE;
current.lp = lp;
candidate.isdual = TRUE;
candidate.lp = lp;
/* Loop over active partial row set */
if(is_action(lp->piv_strategy, PRICE_FORCEFULL)) {
k = 1;
iy = lp->rows;
}
else {
k = partial_blockStart(lp, TRUE);
iy = partial_blockEnd(lp, TRUE);
}
ninfeas = 0;
xinfeas = 0;
sinfeas = 0;
makePriceLoop(lp, &k, &iy, &iz);
iy *= iz;
for(; k*iz <= iy; k += iz) {
/* Map loop variable to target */
i = k;
/* Check if the pivot candidate is on the block-list */
if(lp->rejectpivot[0] > 0) {
int kk;
for(kk = 1; (kk <= lp->rejectpivot[0]) && (i != lp->rejectpivot[kk]); kk++);
if(kk <= lp->rejectpivot[0])
continue;
}
/* Set local variables - express violation as a negative number */
ii = lp->var_basic[i];
up = lp->upbo[ii];
lo = 0;
rh = rhvec[i];
if(rh > up)
rh = up - rh;
else
rh -= lo;
up -= lo;
/* Analyze relevant constraints ...
KE version skips uninteresting alternatives and gives a noticeable speedup */
/* if((rh < -epsvalue*sqrt(lp->matA->rowmax[i])) || */
if((rh < -epsvalue) ||
((forceoutEQ == TRUE) && (up < epsvalue))) { /* It causes instability to remove the "TRUE" test */
/* Accumulate stats */
ninfeas++;
SETMIN(xinfeas, rh);
sinfeas += rh;
/* Give a slight preference to fixed variables (mainly equality slacks) */
if(up < epsvalue) {
/* Break out immediately if we are directed to force slacks out of the basis */
if(forceoutEQ == TRUE) {
current.varno = i;
current.pivot = -1;
break;
}
/* Give an extra early boost to equality slack elimination, if specified */
if(forceoutEQ == AUTOMATIC)
rh *= 10.0;
else /* .. or just the normal. marginal boost */
rh *= 1.0+lp->epspivot;
}
/* Select leaving variable according to strategy (the most negative/largest violation) */
candidate.pivot = normalizeEdge(lp, i, rh, TRUE);
candidate.varno = i;
if(findImprovementVar(&current, &candidate, collectMP, NULL))
break;
}
}
/* Verify infeasibility */
if(updateinfeas)
lp->suminfeas = fabs(sinfeas);
if((ninfeas > 1) &&
!verify_stability(lp, FALSE, xinfeas, sinfeas, ninfeas)) {
report(lp, IMPORTANT, "rowdual: Check for reduced accuracy and tolerance settings.\n");
current.varno = 0;
}
/* Produce statistics */
if(lp->spx_trace) {
report(lp, NORMAL, "rowdual: Infeasibility sum " RESULTVALUEMASK " in %7d constraints.\n",
sinfeas, ninfeas);
if(current.varno > 0) {
report(lp, DETAILED, "rowdual: rhs[%d] = " RESULTVALUEMASK "\n",
current.varno, lp->rhs[current.varno]);
}
else
report(lp, FULL, "rowdual: Optimality - No primal infeasibilities found\n");
}
if(xviol != NULL)
*xviol = fabs(xinfeas);
return(current.varno);
} /* rowdual */
STATIC void longdual_testset(lprec *lp, int which, int rownr, REAL *prow, int *nzprow,
REAL *drow, int *nzdrow)
{
int i,j;
REAL F = lp->infinite;
if(which == 0) { /* Maros Example-1 - raw data */
j = 1; i = lp->rows+j; lp->upbo[i] = 0; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 2; drow[i] = -1;
j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 5;
j = 4; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 3; drow[i] = -6;
j = 5; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -4; drow[i] = -2;
j = 6; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 0;
j = 7; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
j = 8; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
j = 9; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 4;
j = 10; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 10;
nzprow[0] = i-lp->rows;
lp->rhs[rownr] = -11;
lp->upbo[lp->var_basic[rownr]] = F;
lp->rhs[0] = 1;
}
else if(which == 1) { /* Maros Example-1 - presorted in correct order */
j = 1; i = lp->rows+j; lp->upbo[i] = 0; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 2; drow[i] = -1;
j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 5;
j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -4; drow[i] = -2;
j = 4; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
j = 5; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 0;
j = 6; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
j = 7; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
j = 8; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 3; drow[i] = -6;
j = 9; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 4;
j = 10; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 10;
nzprow[0] = i-lp->rows;
lp->rhs[rownr] = -11;
lp->upbo[lp->var_basic[rownr]] = F;
lp->rhs[0] = 1;
}
else if(which == 10) { /* Maros Example-2 - raw data */
j = 1; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 3; drow[i] = 3;
j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
j = 4; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -1; drow[i] = -2;
j = 5; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
j = 6; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 3; drow[i] = 9;
nzprow[0] = i-lp->rows;
lp->rhs[rownr] = 14;
lp->upbo[lp->var_basic[rownr]] = 2;
lp->rhs[0] = 6;
}
}
/* Find the dual simplex entering non-basic variable */
STATIC int coldual(lprec *lp, int row_nr, REAL *prow, int *nzprow,
REAL *drow, int *nzdrow,
MYBOOL dualphase1, MYBOOL skipupdate,
int *candidatecount, REAL *xviol)
{
int i, iy, iz, ix, k, nbound;
LREAL w, g, quot;
REAL viol, p, epspivot = lp->epspivot;
#ifdef MachinePrecRoundRHS
REAL epsvalue = lp->epsmachine;
#else
REAL epsvalue = lp->epsvalue;
#endif
pricerec current, candidate;
MYBOOL isbatch = FALSE, /* Requires that lp->longsteps->size > lp->sum */
dolongsteps = (MYBOOL) (lp->longsteps != NULL);
/* Initialize */
if(xviol != NULL)
*xviol = lp->infinite;
if(dolongsteps && !dualphase1)
dolongsteps = AUTOMATIC; /* Sets Phase1 = TRUE, Phase2 = AUTOMATIC */
current.theta = lp->infinite;
current.pivot = 0;
current.varno = 0;
current.epspivot = epspivot;
current.isdual = TRUE;
current.lp = lp;
candidate.epspivot = epspivot;
candidate.isdual = TRUE;
candidate.lp = lp;
*candidatecount = 0;
/* Compute reduced costs */
if(!skipupdate) {
#ifdef UseDualReducedCostUpdate
/* Recompute from scratch only at the beginning, otherwise update */
if((lp->current_iter > 0) && (refactRecent(lp) < AUTOMATIC))
compute_reducedcosts(lp, TRUE, row_nr, NULL, TRUE,
prow, nzprow,
NULL, NULL,
MAT_ROUNDDEFAULT);
else
#endif
compute_reducedcosts(lp, TRUE, row_nr, NULL, TRUE,
prow, nzprow,
drow, nzdrow,
MAT_ROUNDDEFAULT);
}
#if 0
/* Override all above to do in-line testing with fixed test set */
if(lp->rows > 1 && lp->columns > 10)
longdual_testset(lp, 10, row_nr, prow, nzprow, drow, nzdrow);
#endif
/* Compute the current violation of the bounds of the outgoing variable,
negative for violation of lower bound, positive for upper bound violation.
(Basic variables are always lower-bounded, by lp_solve convention) */
g = 1;
viol = lp->rhs[row_nr];
if(viol > 0) { /* Check if the leaving variable is >= its upper bound */
p = lp->upbo[lp->var_basic[row_nr]];
if(p < lp->infinite) {
viol -= p;
my_roundzero(viol, epsvalue);
if(viol > 0)
g = -1;
}
/* Do validation of numerics */
if(g == 1) {
if(viol >= lp->infinite) {
report(lp, IMPORTANT, "coldual: Large basic solution value %g at iter %.0f indicates numerical instability\n",
lp->rhs[row_nr], (double) get_total_iter(lp));
lp->spx_status = NUMFAILURE;
return( 0 );
}
if(skipupdate)
report(lp, DETAILED, "coldual: Inaccurate bound-flip accuracy at iter %.0f\n",
(double) get_total_iter(lp));
else
report(lp, SEVERE, "coldual: Leaving variable %d does not violate bounds at iter %.0f\n",
row_nr, (double) get_total_iter(lp));
return( -1 );
}
}
/* Update local value of pivot setting */
lp->_piv_rule_ = get_piv_rule(lp);
/* Condense list of relevant targets */
p = 0;
k = 0;
nbound = 0;
ix = 1;
iy = nzprow[0];
for(ix = 1; ix <= iy; ix++) {
i = nzprow[ix];
w = prow[i] * g; /* Change sign if upper bound of the leaving variable is violated */
/* Change sign if the non-basic variable is currently upper-bounded */
/* w *= 2*lp->is_lower[i] - 1; */ /* fails on AIX!!! */
w = my_chsign(!lp->is_lower[i], w);
/* Check if the candidate is worth using for anything */
if(w < -epsvalue) {
/* Tally bounded variables */
if(lp->upbo[i] < lp->infinite)
nbound++;
/* Update the nz-index */
k++;
nzprow[k] = i;
SETMAX(p, -w);
}
#ifdef Paranoia
else {
if(lp->spx_trace) {
report(lp, FULL, "coldual: Candidate variable prow[%d] rejected with %g too small\n",
i, w);
}
}
#endif
}
nzprow[0] = k;
if(xviol != NULL)
*xviol = p;
#ifdef UseRelativePivot_Dual
/* epspivot *= sqrt(lp->matA->dynrange) / lp->matA->infnorm; */
epspivot /= MAX(1, sqrt(lp->matA->rowmax[row_nr]));
#endif
current.epspivot = epspivot;
candidate.epspivot = epspivot;
/* Initialize the long-step structures if indicated */
if(dolongsteps) {
if((nzprow[0] <= 1) || (nbound == 0)) { /* Don't bother */
dolongsteps = FALSE;
lp->longsteps->indexSet[0] = 0;
}
else {
multi_restart(lp->longsteps);
multi_valueInit(lp->longsteps, g*viol, lp->rhs[0]);
}
}
/* Loop over all entering column candidates */
ix = 1;
iy = nzprow[0];
makePriceLoop(lp, &ix, &iy, &iz);
iy *= iz;
for(; ix*iz <= iy; ix += iz) {
i = nzprow[ix];
/* Compute the dual ratio (prow = w and drow = cbar in Chvatal's "nomenclatura") */
w = prow[i] * g; /* Change sign if upper bound of the leaving variable is violated */
quot = -drow[i] / w; /* Remember this sign-reversal in multi_recompute! */
/* Apply the selected pivot strategy (smallest theta) */
candidate.theta = quot; /* Note that abs() is applied in findSubstitutionVar */
candidate.pivot = w;
candidate.varno = i;
/* Collect candidates for minor iterations/bound flips */
if(dolongsteps) {
if(isbatch && (ix == iy))
isbatch = AUTOMATIC;
if(collectMinorVar(&candidate, lp->longsteps, (MYBOOL) (dolongsteps == AUTOMATIC), isbatch) &&
lp->spx_trace)
report(lp, DETAILED, "coldual: Long-dual break point with %d bound-flip variables\n",
lp->longsteps->used);
if(lp->spx_status == FATHOMED)
return( 0 );
}
/* We have a candidate for entering the basis; check if it is better than the incumbent */
else if(findSubstitutionVar(&current, &candidate, candidatecount))
break;
}
/* Set entering variable and long-step bound swap variables */
if(dolongsteps) {
*candidatecount = lp->longsteps->used;
i = multi_enteringvar(lp->longsteps, NULL, 3);
}
else
i = current.varno;
if(lp->spx_trace)
report(lp, NORMAL, "coldual: Entering column %d, reduced cost %g, pivot value %g, bound swaps %d\n",
i, drow[i], prow[i], multi_used(lp->longsteps));
return( i );
} /* coldual */
INLINE REAL normalizeEdge(lprec *lp, int item, REAL edge, MYBOOL isdual)
{
#if 1
/* Don't use the pricer "close to home", since this can possibly
worsen the final feasibility picture (mainly a Devex issue?) */
if(fabs(edge) > lp->epssolution)
#endif
edge /= getPricer(lp, item, isdual);
if((lp->piv_strategy & PRICE_RANDOMIZE) != 0)
edge *= (1.0-PRICER_RANDFACT) + PRICER_RANDFACT*rand_uniform(lp, 1.0);
return( edge );
}
/* Support routines for block detection and partial pricing */
STATIC int partial_findBlocks(lprec *lp, MYBOOL autodefine, MYBOOL isrow)
{
int i, jj, n, nb, ne, items;
REAL hold, biggest, *sum = NULL;
MATrec *mat = lp->matA;
partialrec *blockdata;
if(!mat_validate(mat))
return( 1 );
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
items = IF(isrow, lp->rows, lp->columns);
allocREAL(lp, &sum, items+1, FALSE);
/* Loop over items and compute the average column index for each */
sum[0] = 0;
for(i = 1; i <= items; i++) {
n = 0;
if(isrow) {
nb = mat->row_end[i-1];
ne = mat->row_end[i];
}
else {
nb = mat->col_end[i-1];
ne = mat->col_end[i];
}
n = ne-nb;
sum[i] = 0;
if(n > 0) {
if(isrow)
for(jj = nb; jj < ne; jj++)
sum[i] += ROW_MAT_COLNR(jj);
else
for(jj = nb; jj < ne; jj++)
sum[i] += COL_MAT_ROWNR(jj);
sum[i] /= n;
}
else
sum[i] = sum[i-1];
}
/* Loop over items again, find largest difference and make monotone */
hold = 0;
biggest = 0;
for(i = 2; i <= items; i++) {
hold = sum[i] - sum[i-1];
if(hold > 0) {
if(hold > biggest)
biggest = hold;
}
else
hold = 0;
sum[i-1] = hold;
}
/* Loop over items again and find differences exceeding threshold;
the discriminatory power of this routine depends strongly on the
magnitude of the scaling factor - from empirical evidence > 0.9 */
biggest = MAX(1, 0.9*biggest);
n = 0;
nb = 0;
ne = 0;
for(i = 1; i < items; i++)
if(sum[i] > biggest) {
ne += i-nb; /* Compute sum of index gaps between maxima */
nb = i;
n++; /* Increment count */
}
/* Clean up */
FREE(sum);
/* Require that the maxima are spread "nicely" across the columns,
otherwise return that there is only one monolithic block.
(This is probably an area for improvement in the logic!) */
if(n > 0) {
ne /= n; /* Average index gap between maxima */
i = IF(isrow, lp->columns, lp->rows);
nb = i / ne; /* Another estimated block count */
if(abs(nb - n) > 2) /* Probably Ok to require equality (nb==n)*/
n = 1;
else if(autodefine) /* Generate row/column break-indeces for partial pricing */
set_partialprice(lp, nb, NULL, isrow);
}
else
n = 1;
return( n );
}
STATIC int partial_blockStart(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( 1 );
else {
if((blockdata->blocknow < 1) || (blockdata->blocknow > blockdata->blockcount))
blockdata->blocknow = 1;
return( blockdata->blockend[blockdata->blocknow-1] );
}
}
STATIC int partial_blockEnd(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( IF(isrow, lp->rows, lp->sum) );
else {
if((blockdata->blocknow < 1) || (blockdata->blocknow > blockdata->blockcount))
blockdata->blocknow = 1;
return( blockdata->blockend[blockdata->blocknow]-1 );
}
}
STATIC int partial_blockNextPos(lprec *lp, int block, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
#ifdef Paranoia
if((blockdata == NULL) || (block <= 1) || (block > blockdata->blockcount)) {
report(lp, SEVERE, "partial_blockNextPos: Invalid block %d specified.\n",
block);
return( -1 );
}
#endif
block--;
if(blockdata->blockpos[block] == blockdata->blockend[block+1])
blockdata->blockpos[block] = blockdata->blockend[block];
else
blockdata->blockpos[block]++;
return( blockdata->blockpos[block] );
}
STATIC MYBOOL partial_blockStep(lprec *lp, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( FALSE );
else if(blockdata->blocknow < blockdata->blockcount) {
blockdata->blocknow++;
return( TRUE);
}
else {
blockdata->blocknow = 1;
return( TRUE );
}
}
STATIC MYBOOL partial_isVarActive(lprec *lp, int varno, MYBOOL isrow)
{
partialrec *blockdata;
blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
if(blockdata == NULL)
return( TRUE );
else {
return( (MYBOOL) ((varno >= blockdata->blockend[blockdata->blocknow-1]) &&
(varno < blockdata->blockend[blockdata->blocknow])) );
}
}
/* Multiple pricing routines */
STATIC multirec *multi_create(lprec *lp, MYBOOL truncinf)
{
multirec *multi;
multi = (multirec *) calloc(1, sizeof(*multi));
if(multi != NULL) {
multi->active = 1;
multi->lp = lp;
multi->epszero = lp->epsprimal;
multi->truncinf = truncinf;
}
return(multi);
}
STATIC void multi_free(multirec **multi)
{
if((multi == NULL) || (*multi == NULL))
return;
FREE((*multi)->items);
FREE((*multi)->valueList);
FREE((*multi)->indexSet);
FREE((*multi)->freeList);
FREE((*multi)->sortedList);
FREE(*multi);
}
STATIC MYBOOL multi_mustupdate(multirec *multi)
{
return( (MYBOOL) ((multi != NULL) &&
(multi->used < multi->limit)) );
}
STATIC MYBOOL multi_resize(multirec *multi, int blocksize, int blockdiv, MYBOOL doVlist, MYBOOL doIset)
{
MYBOOL ok = TRUE;
if((blocksize > 1) && (blockdiv > 0)) {
int oldsize = multi->size;
multi->size = blocksize;
if(blockdiv > 1)
multi->limit += (multi->size-oldsize) / blockdiv;
multi->items = (pricerec *) realloc(multi->items, (multi->size+1)*sizeof(*(multi->items)));
multi->sortedList = (UNIONTYPE QSORTrec *) realloc(multi->sortedList, (multi->size+1)*sizeof(*(multi->sortedList)));
ok = (multi->items != NULL) && (multi->sortedList != NULL) &&
allocINT(multi->lp, &(multi->freeList), multi->size+1, AUTOMATIC);
if(ok) {
int i, n;
if(oldsize == 0)
i = 0;
else
i = multi->freeList[0];
multi->freeList[0] = i + (multi->size-oldsize);
for(n = multi->size - 1, i++; i <= multi->freeList[0]; i++, n--)
multi->freeList[i] = n;
}
if(doVlist)
ok &= allocREAL(multi->lp, &(multi->valueList), multi->size+1, AUTOMATIC);
if(doIset) {
ok &= allocINT(multi->lp, &(multi->indexSet), multi->size+1, AUTOMATIC);
if(ok && (oldsize == 0))
multi->indexSet[0] = 0;
}
if(!ok)
goto Undo;
}
else {
Undo:
multi->size = 0;
FREE(multi->items);
FREE(multi->valueList);
FREE(multi->indexSet);
FREE(multi->freeList);
FREE(multi->sortedList);
}
multi->active = 1;
return( ok );
}
STATIC int multi_size(multirec *multi)
{
if(multi == NULL)
return( 0 );
else
return( multi->size );
}
STATIC int multi_used(multirec *multi)
{
if(multi == NULL)
return( 0 );
else
return( multi->used );
}
STATIC int multi_restart(multirec *multi)
{
int i, n = multi->used;
multi->used = 0;
multi->sorted = FALSE;
multi->dirty = FALSE;
if(multi->freeList != NULL) {
for(i = 1; i <= multi->size; i++)
multi->freeList[i] = multi->size - i;
multi->freeList[0] = multi->size;
}
#if 0
if(multi->indexSet != NULL)
multi->indexSet[0] = 0;
#endif
return( n );
}
STATIC void multi_valueInit(multirec *multi, REAL step_base, REAL obj_base)
{
multi->step_base = multi->step_last = step_base;
multi->obj_base = multi->obj_last = obj_base;
#ifdef Paranoia
if(step_base > 0)
report(multi->lp, SEVERE, "multi_valueInit: Positive constraint violation %g provided at iteration %6.0f\n",
step_base, (double) get_total_iter(multi->lp));
#endif
}
STATIC REAL *multi_valueList(multirec *multi)
{
return(multi->valueList);
}
STATIC int *multi_indexSet(multirec *multi, MYBOOL regenerate)
{
if(regenerate)
multi_populateSet(multi, NULL, -1);
return(multi->indexSet);
}
STATIC int multi_getvar(multirec *multi, int item)
{
#ifdef Paranoia
if((item < 1) || (item >= multi->size))
return(-1);
#endif
return( ((pricerec *) &(multi->sortedList[item].pvoidreal.ptr))->varno );
}
STATIC MYBOOL multi_recompute(multirec *multi, int index, MYBOOL isphase2, MYBOOL fullupdate)
{
int i, n;
REAL lB, uB, Alpha, this_theta, prev_theta;
lprec *lp = multi->lp;
pricerec *thisprice;
/* Define target update window */
if(multi->dirty) {
index = 0;
n = multi->used - 1;
}
else if(fullupdate)
n = multi->used - 1;
else
n = index;
/* Initialize accumulators from the specified update index */
if(index == 0) {
multi->maxpivot = 0;
multi->maxbound = 0;
multi->step_last = multi->step_base;
multi->obj_last = multi->obj_base;
thisprice = NULL;
this_theta = 0;
}
else {
multi->obj_last = multi->valueList[index-1];
multi->step_last = multi->sortedList[index-1].pvoidreal.realval;
thisprice = (pricerec *) (multi->sortedList[index-1].pvoidreal.ptr);
this_theta = thisprice->theta;
}
/* Update step lengths and objective values */
while((index <= n) && (multi->step_last < multi->epszero)) {
/* Update parameters for this loop */
prev_theta = this_theta;
thisprice = (pricerec *) (multi->sortedList[index].pvoidreal.ptr);
this_theta = thisprice->theta;
Alpha = fabs(thisprice->pivot);
uB = lp->upbo[thisprice->varno];
lB = 0;
SETMAX(multi->maxpivot, Alpha);
SETMAX(multi->maxbound, uB);
/* Do the value updates */
if(isphase2) {
multi->obj_last += (this_theta - prev_theta) * multi->step_last; /* Sign-readjusted from coldual()/Maros */
if(uB >= lp->infinite)
multi->step_last = lp->infinite;
else
multi->step_last += Alpha*(uB-lB);
}
else {
multi->obj_last += (this_theta - prev_theta) * multi->step_last; /* Sign-readjusted from coldual()/Maros */
multi->step_last += Alpha;
}
/* Store updated values at the indexed locations */
multi->sortedList[index].pvoidreal.realval = multi->step_last;
multi->valueList[index] = multi->obj_last;
#ifdef Paranoia
if(lp->spx_trace &&
(multi->step_last > lp->infinite))
report(lp, SEVERE, "multi_recompute: A very large step-size %g was generated at iteration %6.0f\n",
multi->step_last, (double) get_total_iter(lp));
#endif
index++;
}
/* Discard candidates entered earlier that now make the OF worsen, and
make sure that the released positions are added to the free list. */
n = index;
while(n < multi->used) {
i = ++multi->freeList[0];
multi->freeList[i] = (int) (((pricerec *) multi->sortedList[n].pvoidreal.ptr) - multi->items);
n++;
}
multi->used = index;
if(multi->sorted && (index == 1))
multi->sorted = FALSE;
multi->dirty = FALSE;
/* Return TRUE if the step is now positive */
return( (MYBOOL) (multi->step_last >= multi->epszero) );
}
STATIC MYBOOL multi_truncatingvar(multirec *multi, int varnr)
{
return( multi->truncinf && is_infinite(multi->lp, multi->lp->upbo[varnr]) );
}
STATIC MYBOOL multi_removevar(multirec *multi, int varnr)
{
int i = 1;
int *coltarget = multi->indexSet;
if(coltarget == NULL)
return( FALSE );
while((i <= multi->used) && (coltarget[i] != varnr))
i++;
if(i > multi->used)
return( FALSE );
for(; i < multi->used; i++)
coltarget[i] = coltarget[i+1];
coltarget[0]--;
multi->used--;
multi->dirty = TRUE;
return( TRUE );
}
STATIC int multi_enteringvar(multirec *multi, pricerec *current, int priority)
{
lprec *lp = multi->lp;
int i = 0, bestindex, colnr;
REAL bound, score, bestscore = -lp->infinite;
REAL b1, b2, b3;
pricerec *candidate, *bestcand;
/* Check that we have a candidate */
multi->active = bestindex = 0;
if((multi == NULL) || (multi->used == 0))
return( bestindex );
/* Check for pruning possibility of the B&B tree */
if(multi->objcheck && (lp->solutioncount > 0) &&
bb_better(lp, OF_WORKING | OF_PROJECTED, OF_TEST_WE)) {
lp->spx_status = FATHOMED;
return( bestindex );
}
/* Check the trivial case */
if(multi->used == 1) {
bestcand = (pricerec *) (multi->sortedList[bestindex].pvoidreal.ptr);
goto Finish;
}
/* Set priority weights */
Redo:
switch(priority) {
case 0: b1 = 0.0, b2 = 0.0, b3 = 1.0; /* Only OF */
bestindex = multi->used - 2; break;
case 1: b1 = 0.2, b2 = 0.3, b3 = 0.5; break; /* Emphasize OF */
case 2: b1 = 0.3, b2 = 0.5, b3 = 0.2; break; /* Emphasize bound */
case 3: b1 = 0.6, b2 = 0.2, b3 = 0.2; break; /* Emphasize pivot */
case 4: b1 = 1.0, b2 = 0.0, b3 = 0.0; break; /* Only pivot */
default: b1 = 0.4, b2 = 0.2, b3 = 0.4; /* Balanced default */
}
bestcand = (pricerec *) (multi->sortedList[bestindex].pvoidreal.ptr);
/* Loop over all candidates to get the best entering candidate;
start at the end to try to maximize the chain length */
for(i = multi->used - 1; i >= 0; i--) {
candidate = (pricerec *) (multi->sortedList[i].pvoidreal.ptr);
colnr = candidate->varno;
bound = lp->upbo[colnr];
score = fabs(candidate->pivot) / multi->maxpivot;
score = pow(1.0 + score , b1) *
pow(1.0 + log(bound / multi->maxbound + 1), b2) *
pow(1.0 + (REAL) i / multi->used , b3);
if(score > bestscore) {
bestscore = score;
bestindex = i;
bestcand = candidate;
}
}
/* Do pivot protection */
if((priority < 4) && (fabs(bestcand->pivot) < lp->epssolution)) {
bestindex = 0;
priority++;
goto Redo;
}
Finish:
/* Make sure we shrink the list and update */
multi->active = colnr = bestcand->varno;
if(bestindex < multi->used - 1) {
#if 0
/* if(lp->upbo[colnr] >= lp->infinite) */
QS_swap(multi->sortedList, bestindex, multi->used-1);
multi_recompute(multi, bestindex, (bestcand->isdual == AUTOMATIC), TRUE);
#else
multi->used = i + 1;
#endif
}
multi_populateSet(multi, NULL, multi->active);
/* Compute the entering theta and update parameters */
score = (multi->used == 1 ? multi->step_base : multi->sortedList[multi->used-2].pvoidreal.realval);
score /= bestcand->pivot;
score = my_chsign(!lp->is_lower[multi->active], score);
if(lp->spx_trace &&
(fabs(score) > 1/lp->epsprimal))
report(lp, IMPORTANT, "multi_enteringvar: A very large Theta %g was generated (pivot %g)\n",
score, bestcand->pivot);
multi->step_base = score;
if(current != NULL)
*current = *bestcand;
return( multi->active );
}
STATIC REAL multi_enteringtheta(multirec *multi)
{
return( multi->step_base );
}
STATIC int multi_populateSet(multirec *multi, int **list, int excludenr)
{
int n = 0;
if(list == NULL)
list = &(multi->indexSet);
if((multi->used > 0) &&
((*list != NULL) || allocINT(multi->lp, list, multi->size+1, FALSE))) {
int i, colnr;
for(i = 0; i < multi->used; i++) {
colnr = ((pricerec *) (multi->sortedList[i].pvoidreal.ptr))->varno;
if((colnr != excludenr) &&
/* Prevent an unbounded variable from "bound-flip"; this could
actually indicate that we should let the entering variable be
bound-swapped (in the case that it is bounded), but we
disregard this possibility here, since it brings with it
issues of pivot size, etc. */
((excludenr > 0) && (multi->lp->upbo[colnr] < multi->lp->infinite))) {
n++;
(*list)[n] = colnr;
}
}
(*list)[0] = n;
}
return( n );
}