limp-cbc-0.3.2.0: cbits/coin/IdiSolve.cpp
/* $Id: IdiSolve.cpp 1878 2012-08-30 15:43:19Z forrest $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#include "CoinPragma.hpp"
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "Idiot.hpp"
#define FIT
#ifdef FIT
#define HISTORY 8
#else
#define HISTORY 7
#endif
#define NSOLVE HISTORY-1
static void solveSmall(int nsolve, double **aIn, double **a, double * b)
{
int i, j;
/* copy */
for (i = 0; i < nsolve; i++) {
for (j = 0; j < nsolve; j++) {
a[i][j] = aIn[i][j];
}
}
for (i = 0; i < nsolve; i++) {
/* update using all previous */
double diagonal;
int j;
for (j = i; j < nsolve; j++) {
int k;
double value = a[i][j];
for (k = 0; k < i; k++) {
value -= a[k][i] * a[k][j];
}
a[i][j] = value;
}
diagonal = a[i][i];
if (diagonal < 1.0e-20) {
diagonal = 0.0;
} else {
diagonal = 1.0 / sqrt(diagonal);
}
a[i][i] = diagonal;
for (j = i + 1; j < nsolve; j++) {
a[i][j] *= diagonal;
}
}
/* update */
for (i = 0; i < nsolve; i++) {
int j;
double value = b[i];
for (j = 0; j < i; j++) {
value -= b[j] * a[j][i];
}
value *= a[i][i];
b[i] = value;
}
for (i = nsolve - 1; i >= 0; i--) {
int j;
double value = b[i];
for (j = i + 1; j < nsolve; j++) {
value -= b[j] * a[i][j];
}
value *= a[i][i];
b[i] = value;
}
}
IdiotResult
Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol,
double * pi, double * /*djs*/, const double * cost ,
const double * /*rowlower*/,
const double * rowupper, const double * /*lower*/,
const double * /*upper*/, const double * elemnt,
const int * row, const CoinBigIndex * columnStart,
const int * length, int extraBlock, int * rowExtra,
double * solExtra, double * elemExtra, double * /*upperExtra*/,
double * costExtra, double weight)
{
IdiotResult result;
double objvalue = 0.0;
double sum1 = 0.0, sum2 = 0.0;
int i;
for (i = 0; i < nrows; i++) {
rowsol[i] = -rowupper[i];
}
for (i = 0; i < ncols; i++) {
CoinBigIndex j;
double value = colsol[i];
if (value) {
objvalue += value * cost[i];
if (elemnt) {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += elemnt[j] * value;
}
} else {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += value;
}
}
}
}
/* adjust to make as feasible as possible */
/* no */
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
double element = elemExtra[i];
int irow = rowExtra[i];
objvalue += solExtra[i] * costExtra[i];
rowsol[irow] += solExtra[i] * element;
}
}
for (i = 0; i < nrows; i++) {
double value = rowsol[i];
sum1 += fabs(value);
sum2 += value * value;
pi[i] = -2.0 * weight * value;
}
result.infeas = sum1;
result.objval = objvalue;
result.weighted = objvalue + weight * sum2;
result.sumSquared = sum2;
return result;
}
IdiotResult
Idiot::IdiSolve(
int nrows, int ncols, double * COIN_RESTRICT rowsol , double * COIN_RESTRICT colsol,
double * COIN_RESTRICT pi, double * COIN_RESTRICT djs, const double * COIN_RESTRICT origcost , double * COIN_RESTRICT rowlower,
double * COIN_RESTRICT rowupper, const double * COIN_RESTRICT lower,
const double * COIN_RESTRICT upper, const double * COIN_RESTRICT elemnt,
const int * row, const CoinBigIndex * columnStart,
const int * length, double * COIN_RESTRICT lambda,
int maxIts, double mu, double drop,
double maxmin, double offset,
int strategy, double djTol, double djExit, double djFlag,
CoinThreadRandom * randomNumberGenerator)
{
IdiotResult result;
int i, j, k, iter;
double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
double tolerance = 1.0e-8;
double * history[HISTORY+1];
int ncolx;
int nChange;
int extraBlock = 0;
int * rowExtra = NULL;
double * COIN_RESTRICT solExtra = NULL;
double * COIN_RESTRICT elemExtra = NULL;
double * COIN_RESTRICT upperExtra = NULL;
double * COIN_RESTRICT costExtra = NULL;
double * COIN_RESTRICT useCostExtra = NULL;
double * COIN_RESTRICT saveExtra = NULL;
double * COIN_RESTRICT cost = NULL;
double saveValue = 1.0e30;
double saveOffset = offset;
double useOffset = offset;
/*#define NULLVECTOR*/
#ifndef NULLVECTOR
int nsolve = NSOLVE;
#else
int nsolve = NSOLVE + 1; /* allow for null vector */
#endif
int nflagged;
double * COIN_RESTRICT thetaX;
double * COIN_RESTRICT djX;
double * COIN_RESTRICT bX;
double * COIN_RESTRICT vX;
double ** aX;
double **aworkX;
double ** allsum;
double * COIN_RESTRICT saveSol = 0;
const double * COIN_RESTRICT useCost = cost;
double bestSol = 1.0e60;
double weight = 0.5 / mu;
char * statusSave = new char[2*ncols];
char * statusWork = statusSave + ncols;
#define DJTEST 5
double djSave[DJTEST];
double largestDj = 0.0;
double smallestDj = 1.0e60;
double maxDj = 0.0;
int doFull = 0;
#define SAVEHISTORY 10
#define EVERY (2*SAVEHISTORY)
#define AFTER SAVEHISTORY*(HISTORY+1)
#define DROP 5
double after = AFTER;
double obj[DROP];
double kbad = 0, kgood = 0;
if (strategy & 128) after = 999999; /* no acceleration at all */
for (i = 0; i < DROP; i++) {
obj[i] = 1.0e70;
}
//#define FOUR_GOES 2
#ifdef FOUR_GOES
double * COIN_RESTRICT pi2 = new double [3*nrows];
double * COIN_RESTRICT rowsol2 = new double [3*nrows];
double * COIN_RESTRICT piX[4];
double * COIN_RESTRICT rowsolX[4];
int startsX[2][5];
int nChangeX[4];
double maxDjX[4];
double objvalueX[4];
int nflaggedX[4];
piX[0]=pi;
piX[1]=pi2;
piX[2]=pi2+nrows;
piX[3]=piX[2]+nrows;
rowsolX[0]=rowsol;
rowsolX[1]=rowsol2;
rowsolX[2]=rowsol2+nrows;
rowsolX[3]=rowsolX[2]+nrows;
#endif
allsum = new double * [nsolve];
aX = new double * [nsolve];
aworkX = new double * [nsolve];
thetaX = new double[nsolve];
vX = new double[nsolve];
bX = new double[nsolve];
djX = new double[nsolve];
allsum[0] = pi;
for (i = 0; i < nsolve; i++) {
if (i) allsum[i] = new double[nrows];
aX[i] = new double[nsolve];
aworkX[i] = new double[nsolve];
}
/* check = rows */
for (i = 0; i < nrows; i++) {
if (rowupper[i] - rowlower[i] > tolerance) {
extraBlock++;
}
}
cost = new double[ncols];
memset(rowsol, 0, nrows * sizeof(double));
for (i = 0; i < ncols; i++) {
CoinBigIndex j;
double value = origcost[i] * maxmin;
double value2 = colsol[i];
if (elemnt) {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
value += elemnt[j] * lambda[irow];
rowsol[irow] += elemnt[j] * value2;
}
} else {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
value += lambda[irow];
rowsol[irow] += value2;
}
}
cost[i] = value;
}
if (extraBlock) {
rowExtra = new int[extraBlock];
solExtra = new double[extraBlock];
elemExtra = new double[extraBlock];
upperExtra = new double[extraBlock];
costExtra = new double[extraBlock];
saveExtra = new double[extraBlock];
extraBlock = 0;
int nbad = 0;
for (i = 0; i < nrows; i++) {
if (rowupper[i] - rowlower[i] > tolerance) {
double smaller, difference;
double value;
saveExtra[extraBlock] = rowupper[i];
if (fabs(rowupper[i]) > fabs(rowlower[i])) {
smaller = rowlower[i];
value = -1.0;
} else {
smaller = rowupper[i];
value = 1.0;
}
if (fabs(smaller) > 1.0e10) {
if (!nbad)
COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
i, smaller));
nbad++;
if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
abort();
if (fabs(rowupper[i]) > fabs(rowlower[i])) {
rowlower[i] = -0.9e10;
smaller = rowlower[i];
value = -1.0;
} else {
rowupper[i] = 0.9e10;
saveExtra[extraBlock] = rowupper[i];
smaller = rowupper[i];
value = 1.0;
}
}
difference = rowupper[i] - rowlower[i];
difference = CoinMin(difference, 1.0e31);
rowupper[i] = smaller;
elemExtra[extraBlock] = value;
solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
if (solExtra[extraBlock] < 0.0) solExtra[extraBlock] = 0.0;
if (solExtra[extraBlock] > difference) solExtra[extraBlock] = difference;
costExtra[extraBlock] = lambda[i] * value;
upperExtra[extraBlock] = difference;
rowsol[i] += value * solExtra[extraBlock];
rowExtra[extraBlock++] = i;
}
}
if (nbad)
COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad));
}
for (i = 0; i < nrows; i++) {
offset += lambda[i] * rowsol[i];
}
if ((strategy & 256) != 0) {
/* save best solution */
saveSol = new double[ncols];
CoinMemcpyN(colsol, ncols, saveSol);
if (extraBlock) {
useCostExtra = new double[extraBlock];
memset(useCostExtra, 0, extraBlock * sizeof(double));
}
useCost = origcost;
useOffset = saveOffset;
} else {
useCostExtra = costExtra;
useCost = cost;
useOffset = offset;
}
ncolx = ncols + extraBlock;
for (i = 0; i < HISTORY + 1; i++) {
history[i] = new double[ncolx];
}
for (i = 0; i < DJTEST; i++) {
djSave[i] = 1.0e30;
}
for (i = 0; i < ncols; i++) {
if (upper[i] - lower[i]) {
statusSave[i] = 0;
} else {
statusSave[i] = 1;
}
}
// for two pass method
int start[2];
int stop[2];
int direction = -1;
start[0] = 0;
stop[0] = ncols;
start[1] = 0;
stop[1] = 0;
iter = 0;
for (; iter < maxIts; iter++) {
double sum1 = 0.0, sum2 = 0.0;
double lastObj = 1.0e70;
int good = 0, doScale = 0;
if (strategy & 16) {
int ii = iter / EVERY + 1;
ii = ii * EVERY;
if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
double * COIN_RESTRICT x = history[HISTORY-1];
for (i = HISTORY - 1; i > 0; i--) {
history[i] = history[i-1];
}
history[0] = x;
CoinMemcpyN(colsol, ncols, history[0]);
CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
}
}
if ((iter % SAVEHISTORY) == 0 || doFull) {
if ((strategy & 16) == 0) {
double * COIN_RESTRICT x = history[HISTORY-1];
for (i = HISTORY - 1; i > 0; i--) {
history[i] = history[i-1];
}
history[0] = x;
CoinMemcpyN(colsol, ncols, history[0]);
CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
}
}
/* start full try */
if ((iter % EVERY) == 0 || doFull) {
// for next pass
direction = - direction;
// randomize.
// The cast is to avoid gcc compiler warning
int kcol = static_cast<int>(ncols * randomNumberGenerator->randomDouble());
if (kcol == ncols)
kcol = ncols - 1;
if (direction > 0) {
start[0] = kcol;
stop[0] = ncols;
start[1] = 0;
stop[1] = kcol;
#ifdef FOUR_GOES
for (int itry=0;itry<2;itry++) {
int chunk=(stop[itry]-start[itry]+FOUR_GOES-1)/FOUR_GOES;
startsX[itry][0]=start[itry];
for (int i=1;i<5;i++)
startsX[itry][i]=CoinMin(stop[itry],startsX[itry][i-1]+chunk);
}
#endif
} else {
start[0] = kcol;
stop[0] = -1;
start[1] = ncols - 1;
stop[1] = kcol;
#ifdef FOUR_GOES
for (int itry=0;itry<2;itry++) {
int chunk=(start[itry]-stop[itry]+FOUR_GOES-1)/FOUR_GOES;
startsX[itry][0]=start[itry];
for (int i=1;i<5;i++)
startsX[itry][i]=CoinMax(stop[itry],startsX[itry][i-1]-chunk);
}
#endif
}
int itry = 0;
/*if ((strategy&16)==0) {
double * COIN_RESTRICT x=history[HISTORY-1];
for (i=HISTORY-1;i>0;i--) {
history[i]=history[i-1];
}
history[0]=x;
CoinMemcpyN(colsol,ncols,history[0]);
CoinMemcpyN(solExtra,extraBlock,history[0]+ncols);
}*/
while (!good) {
itry++;
#define MAXTRY 5
if (iter > after && doScale < 2 && itry < MAXTRY) {
/* now full one */
for (i = 0; i < nrows; i++) {
rowsol[i] = -rowupper[i];
}
sum2 = 0.0;
objvalue = 0.0;
memset(pi, 0, nrows * sizeof(double));
{
double * COIN_RESTRICT theta = thetaX;
double * COIN_RESTRICT dj = djX;
double * COIN_RESTRICT b = bX;
double ** a = aX;
double ** awork = aworkX;
double * COIN_RESTRICT v = vX;
double c;
#ifdef FIT
int ntot = 0, nsign = 0, ngood = 0, mgood[4] = {0, 0, 0, 0};
double diff1, diff2, val0, val1, val2, newValue;
CoinMemcpyN(colsol, ncols, history[HISTORY-1]);
CoinMemcpyN(solExtra, extraBlock, history[HISTORY-1] + ncols);
#endif
dj[0] = 0.0;
for (i = 1; i < nsolve; i++) {
dj[i] = 0.0;
memset(allsum[i], 0, nrows * sizeof(double));
}
for (i = 0; i < ncols; i++) {
double value2 = colsol[i];
if (value2 > lower[i] + tolerance) {
if(value2 < (upper[i] - tolerance)) {
int k;
objvalue += value2 * cost[i];
#ifdef FIT
ntot++;
val0 = history[0][i];
val1 = history[1][i];
val2 = history[2][i];
diff1 = val0 - val1;
diff2 = val1 - val2;
if (diff1*diff2 >= 0.0) {
nsign++;
if (fabs(diff1) < fabs(diff2)) {
int ii = static_cast<int>(fabs(4.0 * diff1 / diff2));
if (ii == 4) ii = 3;
mgood[ii]++;
ngood++;
}
if (fabs(diff1) < 0.75 * fabs(diff2)) {
newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
} else {
newValue = val1 + 4.0 * diff1;
}
} else {
newValue = 0.333333333 * (val0 + val1 + val2);
}
if (newValue > upper[i] - tolerance) {
newValue = upper[i];
} else if (newValue < lower[i] + tolerance) {
newValue = lower[i];
}
history[HISTORY-1][i] = newValue;
#endif
for (k = 0; k < HISTORY - 1; k++) {
value = history[k][i] - history[k+1][i];
dj[k] += value * cost[i];
v[k] = value;
}
if (elemnt) {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
for (k = 0; k < HISTORY - 1; k++) {
allsum[k][irow] += elemnt[j] * v[k];
}
rowsol[irow] += elemnt[j] * value2;
}
} else {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
for (k = 0; k < HISTORY - 1; k++) {
allsum[k][irow] += v[k];
}
rowsol[irow] += value2;
}
}
} else {
/* at ub */
colsol[i] = upper[i];
value2 = colsol[i];
objvalue += value2 * cost[i];
if (elemnt) {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += elemnt[j] * value2;
}
} else {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += value2;
}
}
}
} else {
/* at lb */
if (value2) {
objvalue += value2 * cost[i];
if (elemnt) {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += elemnt[j] * value2;
}
} else {
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
rowsol[irow] += value2;
}
}
}
}
}
#ifdef FIT
/*printf("total %d, same sign %d, good %d %d %d %d %d\n",
ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
#endif
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
double element = elemExtra[i];
int irow = rowExtra[i];
objvalue += solExtra[i] * costExtra[i];
if (solExtra[i] > tolerance
&& solExtra[i] < (upperExtra[i] - tolerance)) {
double value2 = solExtra[i];
int k;
for (k = 0; k < HISTORY - 1; k++) {
value = history[k][i+ncols] - history[k+1][i+ncols];
dj[k] += value * costExtra[i];
allsum[k][irow] += element * value;
}
rowsol[irow] += element * value2;
} else {
double value2 = solExtra[i];
double element = elemExtra[i];
int irow = rowExtra[i];
rowsol[irow] += element * value2;
}
}
}
#ifdef NULLVECTOR
if ((strategy & 64)) {
double djVal = dj[0];
for (i = 0; i < ncols - nrows; i++) {
double value2 = colsol[i];
if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
value = history[0][i] - history[1][i];
} else {
value = 0.0;
}
history[HISTORY][i] = value;
}
for (; i < ncols; i++) {
double value2 = colsol[i];
double delta;
int irow = i - (ncols - nrows);
double oldSum = allsum[0][irow];;
if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
delta = history[0][i] - history[1][i];
} else {
delta = 0.0;
}
djVal -= delta * cost[i];
oldSum -= delta;
delta = - oldSum;
djVal += delta * cost[i];
history[HISTORY][i] = delta;
}
dj[HISTORY-1] = djVal;
djVal = 0.0;
for (i = 0; i < ncols; i++) {
double value2 = colsol[i];
if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance ||
i >= ncols - nrows) {
int k;
value = history[HISTORY][i];
djVal += value * cost[i];
for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
int irow = row[j];
allsum[nsolve-1][irow] += value;
}
}
}
printf("djs %g %g\n", dj[HISTORY-1], djVal);
}
#endif
for (i = 0; i < nsolve; i++) {
int j;
b[i] = 0.0;
for (j = 0; j < nsolve; j++) {
a[i][j] = 0.0;
}
}
c = 0.0;
for (i = 0; i < nrows; i++) {
double value = rowsol[i];
for (k = 0; k < nsolve; k++) {
v[k] = allsum[k][i];
b[k] += v[k] * value;
}
c += value * value;
for (k = 0; k < nsolve; k++) {
for (j = k; j < nsolve; j++) {
a[k][j] += v[k] * v[j];
}
}
}
sum2 = c;
if (itry == 1) {
lastObj = objvalue + weight * sum2;
}
for (k = 0; k < nsolve; k++) {
b[k] = - (weight * b[k] + 0.5 * dj[k]);
for (j = k; j < nsolve; j++) {
a[k][j] *= weight;
a[j][k] = a[k][j];
}
}
c *= weight;
for (k = 0; k < nsolve; k++) {
theta[k] = b[k];
}
solveSmall(nsolve, a, awork, theta);
if ((strategy & 64) != 0) {
value = 10.0;
for (k = 0; k < nsolve; k++) {
value = CoinMax(value, fabs(theta[k]));
}
if (value > 10.0 && ((logLevel_ & 4) != 0)) {
printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
}
value = 10.0 / value;
for (k = 0; k < nsolve; k++) {
theta[k] *= value;
}
}
for (i = 0; i < ncolx; i++) {
double valueh = 0.0;
for (k = 0; k < HISTORY - 1; k++) {
value = history[k][i] - history[k+1][i];
valueh += value * theta[k];
}
#ifdef NULLVECTOR
value = history[HISTORY][i];
valueh += value * theta[HISTORY-1];
#endif
history[HISTORY][i] = valueh;
}
}
#ifdef NULLVECTOR
if ((strategy & 64)) {
for (i = 0; i < ncols - nrows; i++) {
if (colsol[i] <= lower[i] + tolerance
|| colsol[i] >= (upper[i] - tolerance)) {
history[HISTORY][i] = 0.0;;
}
}
tolerance = -tolerance; /* switch off test */
}
#endif
if (!doScale) {
for (i = 0; i < ncols; i++) {
if (colsol[i] > lower[i] + tolerance
&& colsol[i] < (upper[i] - tolerance)) {
value = history[HISTORY][i];
colsol[i] += value;
if (colsol[i] < lower[i] + tolerance) {
colsol[i] = lower[i];
} else if (colsol[i] > upper[i] - tolerance) {
colsol[i] = upper[i];
}
}
}
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
if (solExtra[i] > tolerance
&& solExtra[i] < (upperExtra[i] - tolerance)) {
value = history[HISTORY][i+ncols];
solExtra[i] += value;
if (solExtra[i] < 0.0) {
solExtra[i] = 0.0;
} else if (solExtra[i] > upperExtra[i]) {
solExtra[i] = upperExtra[i];
}
}
}
}
} else {
double theta = 1.0;
double saveTheta = theta;
for (i = 0; i < ncols; i++) {
if (colsol[i] > lower[i] + tolerance
&& colsol[i] < (upper[i] - tolerance)) {
value = history[HISTORY][i];
if (value > 0) {
if (theta * value + colsol[i] > upper[i]) {
theta = (upper[i] - colsol[i]) / value;
}
} else if (value < 0) {
if (colsol[i] + theta * value < lower[i]) {
theta = (lower[i] - colsol[i]) / value;
}
}
}
}
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
if (solExtra[i] > tolerance
&& solExtra[i] < (upperExtra[i] - tolerance)) {
value = history[HISTORY][i+ncols];
if (value > 0) {
if (theta * value + solExtra[i] > upperExtra[i]) {
theta = (upperExtra[i] - solExtra[i]) / value;
}
} else if (value < 0) {
if (solExtra[i] + theta * value < 0.0) {
theta = -solExtra[i] / value;
}
}
}
}
}
if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
if (theta < saveTheta) {
printf(" - modified theta %g\n", theta);
}
}
for (i = 0; i < ncols; i++) {
if (colsol[i] > lower[i] + tolerance
&& colsol[i] < (upper[i] - tolerance)) {
value = history[HISTORY][i];
colsol[i] += value * theta;
}
}
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
if (solExtra[i] > tolerance
&& solExtra[i] < (upperExtra[i] - tolerance)) {
value = history[HISTORY][i+ncols];
solExtra[i] += value * theta;
}
}
}
}
#ifdef NULLVECTOR
tolerance = fabs(tolerance); /* switch back on */
#endif
if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
printf("\n");
}
}
good = 1;
result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
rowlower, rowupper, lower, upper,
elemnt, row, columnStart, length, extraBlock, rowExtra,
solExtra, elemExtra, upperExtra, useCostExtra,
weight);
weightedObj = result.weighted;
if (!iter) saveValue = weightedObj;
objvalue = result.objval;
sum1 = result.infeas;
if (saveSol) {
if (result.weighted < bestSol) {
COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter,
result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
bestSol = result.weighted;
CoinMemcpyN(colsol, ncols, saveSol);
}
}
#ifdef FITz
if (iter > after) {
IdiotResult result2;
double ww, oo, ss;
if (extraBlock) abort();
result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
rowlower, rowupper, lower, upper,
elemnt, row, columnStart, extraBlock, rowExtra,
solExtra, elemExtra, upperExtra, costExtra,
weight);
ww = result2.weighted;
oo = result2.objval;
ss = result2.infeas;
printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
if (ww < weightedObj && ww < lastObj) {
printf(" taken");
ntaken++;
saving += weightedObj - ww;
weightedObj = ww;
objvalue = oo;
sum1 = ss;
CoinMemcpyN(row2, nrows, rowsol);
CoinMemcpyN(pi2, nrows, pi);
CoinMemcpyN(sol2, ncols, colsol);
result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
rowlower, rowupper, lower, upper,
elemnt, row, columnStart, extraBlock, rowExtra,
solExtra, elemExtra, upperExtra, costExtra,
weight);
weightedObj = result.weighted;
objvalue = result.objval;
sum1 = result.infeas;
if (ww < weightedObj) abort();
} else {
printf(" not taken");
nottaken++;
}
}
#endif
/*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
if((logLevel_ & 16) != 0 && doScale) {
printf("Weighted objective from %g to %g **** bad move\n",
lastObj, weightedObj);
}
if (doScale) {
good = 1;
}
if ((strategy & 3) == 1) {
good = 0;
if (weightedObj > lastObj + djExit) {
if ((logLevel_ & 16) != 0) {
printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
}
CoinMemcpyN(history[0], ncols, colsol);
CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
good = 1;
}
} else if ((strategy & 3) == 2) {
if (weightedObj > lastObj + 0.1 * maxDj) {
CoinMemcpyN(history[0], ncols, colsol);
CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
doScale++;
good = 0;
}
} else if ((strategy & 3) == 3) {
if (weightedObj > lastObj + 0.001 * maxDj) {
/*doScale++;*/
good = 0;
}
}
}
}
if ((iter % checkFrequency_) == 0) {
double best = weightedObj;
double test = obj[0];
for (i = 1; i < DROP; i++) {
obj[i-1] = obj[i];
if (best > obj[i]) best = obj[i];
}
obj[DROP-1] = best;
if (test - best < drop && (strategy & 8) == 0) {
if ((logLevel_ & 8) != 0) {
printf("Exiting as drop in %d its is %g after %d iterations\n",
DROP * checkFrequency_, test - best, iter);
}
goto RETURN;
}
}
if ((iter % logFreq_) == 0) {
double piSum = 0.0;
for (i = 0; i < nrows; i++) {
piSum += (rowsol[i] + rowupper[i]) * pi[i];
}
if ((logLevel_ & 2) != 0) {
printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
piSum * maxmin - useOffset, maxDj);
}
}
CoinMemcpyN(statusSave, ncols, statusWork);
nflagged = 0;
}
nChange = 0;
doFull = 0;
maxDj = 0.0;
// go through forwards or backwards and starting at odd places
#ifdef FOUR_GOES
for (int i=1;i<FOUR_GOES;i++) {
cilk_spawn memcpy(piX[i], pi, nrows * sizeof(double));
cilk_spawn memcpy(rowsolX[i], rowsol, nrows * sizeof(double));
}
for (int i=0;i<FOUR_GOES;i++) {
nChangeX[i]=0;
maxDjX[i]=0.0;
objvalueX[i]=0.0;
nflaggedX[i]=0;
}
cilk_sync;
#endif
//printf("PASS\n");
#ifdef FOUR_GOES
cilk_for (int iPar = 0; iPar < FOUR_GOES; iPar++) {
double * COIN_RESTRICT pi=piX[iPar];
double * COIN_RESTRICT rowsol = rowsolX[iPar];
for (int itry = 0; itry < 2; itry++) {
int istop;
int istart;
#if 0
int chunk = (start[itry]+stop[itry]+FOUR_GOES-1)/FOUR_GOES;
if (iPar == 0) {
istart=start[itry];
istop=(start[itry]+stop[itry])>>1;
} else {
istart=(start[itry]+stop[itry])>>1;
istop = stop[itry];
}
#endif
#if 0
printf("istart %d istop %d direction %d array %d %d new %d %d\n",
istart,istop,direction,start[itry],stop[itry],
startsX[itry][iPar],startsX[itry][iPar+1]);
#endif
istart=startsX[itry][iPar];
istop=startsX[itry][iPar+1];
#else
for (int itry = 0; itry < 2; itry++) {
int istart = start[itry];
int istop = stop[itry];
#endif
for (int icol=istart; icol != istop; icol += direction) {
if (!statusWork[icol]) {
CoinBigIndex j;
double value = colsol[icol];
double djval = cost[icol];
double djval2, value2;
double theta, a, b, c;
if (elemnt) {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
djval -= elemnt[j] * pi[irow];
}
} else {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
djval -= pi[irow];
}
}
/*printf("xx iter %d seq %d djval %g value %g\n",
iter,i,djval,value);*/
if (djval > 1.0e-5) {
value2 = (lower[icol] - value);
} else {
value2 = (upper[icol] - value);
}
djval2 = djval * value2;
djval = fabs(djval);
if (djval > djTol) {
if (djval2 < -1.0e-4) {
#ifndef FOUR_GOES
nChange++;
if (djval > maxDj) maxDj = djval;
#else
nChangeX[iPar]++;
if (djval > maxDjX[iPar]) maxDjX[iPar] = djval;
#endif
/*if (djval>3.55e6) {
printf("big\n");
}*/
a = 0.0;
b = 0.0;
c = 0.0;
djval2 = cost[icol];
if (elemnt) {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
double value = rowsol[irow];
c += value * value;
a += elemnt[j] * elemnt[j];
b += value * elemnt[j];
}
} else {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
double value = rowsol[irow];
c += value * value;
a += 1.0;
b += value;
}
}
a *= weight;
b = b * weight + 0.5 * djval2;
c *= weight;
/* solve */
theta = -b / a;
#ifndef FOUR_GOES
if ((strategy & 4) != 0) {
double valuep, thetap;
value2 = a * theta * theta + 2.0 * b * theta;
thetap = 2.0 * theta;
valuep = a * thetap * thetap + 2.0 * b * thetap;
if (valuep < value2 + djTol) {
theta = thetap;
kgood++;
} else {
kbad++;
}
}
#endif
if (theta > 0.0) {
if (theta < upper[icol] - colsol[icol]) {
value2 = theta;
} else {
value2 = upper[icol] - colsol[icol];
}
} else {
if (theta > lower[icol] - colsol[icol]) {
value2 = theta;
} else {
value2 = lower[icol] - colsol[icol];
}
}
colsol[icol] += value2;
#ifndef FOUR_GOES
objvalue += cost[icol] * value2;
#else
objvalueX[iPar] += cost[icol] * value2;
#endif
if (elemnt) {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
double value;
rowsol[irow] += elemnt[j] * value2;
value = rowsol[irow];
pi[irow] = -2.0 * weight * value;
}
} else {
for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
int irow = row[j];
double value;
rowsol[irow] += value2;
value = rowsol[irow];
pi[irow] = -2.0 * weight * value;
}
}
} else {
/* dj but at bound */
if (djval > djFlag) {
statusWork[icol] = 1;
#ifndef FOUR_GOES
nflagged++;
#else
nflaggedX[iPar]++;
#endif
}
}
}
}
}
#ifdef FOUR_GOES
}
#endif
}
#ifdef FOUR_GOES
for (int i=0;i<FOUR_GOES;i++) {
nChange += nChangeX[i];
maxDj = CoinMax(maxDj,maxDjX[i]);
objvalue += objvalueX[i];
nflagged += nflaggedX[i];
}
cilk_for (int i = 0; i < nrows; i++) {
#if FOUR_GOES==2
rowsol[i] = 0.5 * (rowsolX[0][i] + rowsolX[1][i]);
pi[i] = 0.5 * (piX[0][i] + piX[1][i]);
#elif FOUR_GOES==3
pi[i] = 0.33333333333333 * (piX[0][i] + piX[1][i]+piX[2][i]);
rowsol[i] = 0.3333333333333 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]);
#else
pi[i] = 0.25 * (piX[0][i] + piX[1][i]+piX[2][i]+piX[3][i]);
rowsol[i] = 0.25 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]+rowsolX[3][i]);
#endif
}
#endif
if (extraBlock) {
for (int i = 0; i < extraBlock; i++) {
double value = solExtra[i];
double djval = costExtra[i];
double djval2, value2;
double element = elemExtra[i];
double theta, a, b, c;
int irow = rowExtra[i];
djval -= element * pi[irow];
/*printf("xxx iter %d extra %d djval %g value %g\n",
iter,irow,djval,value);*/
if (djval > 1.0e-5) {
value2 = -value;
} else {
value2 = (upperExtra[i] - value);
}
djval2 = djval * value2;
if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
nChange++;
a = 0.0;
b = 0.0;
c = 0.0;
djval2 = costExtra[i];
value = rowsol[irow];
c += value * value;
a += element * element;
b += element * value;
a *= weight;
b = b * weight + 0.5 * djval2;
c *= weight;
/* solve */
theta = -b / a;
if (theta > 0.0) {
value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
} else {
value2 = CoinMax(theta, -solExtra[i]);
}
solExtra[i] += value2;
rowsol[irow] += element * value2;
value = rowsol[irow];
pi[irow] = -2.0 * weight * value;
}
}
}
if ((iter % 10) == 2) {
for (int i = DJTEST - 1; i > 0; i--) {
djSave[i] = djSave[i-1];
}
djSave[0] = maxDj;
largestDj = CoinMax(largestDj, maxDj);
smallestDj = CoinMin(smallestDj, maxDj);
for (int i = DJTEST - 1; i > 0; i--) {
maxDj += djSave[i];
}
maxDj = maxDj / static_cast<double> (DJTEST);
if (maxDj < djExit && iter > 50) {
//printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
break;
}
if (nChange < 100) {
djTol *= 0.5;
}
}
}
RETURN:
if (kgood || kbad) {
COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad));
}
result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
rowlower, rowupper, lower, upper,
elemnt, row, columnStart, length, extraBlock, rowExtra,
solExtra, elemExtra, upperExtra, useCostExtra,
weight);
result.djAtBeginning = largestDj;
result.djAtEnd = smallestDj;
result.dropThis = saveValue - result.weighted;
if (saveSol) {
if (result.weighted < bestSol) {
bestSol = result.weighted;
CoinMemcpyN(colsol, ncols, saveSol);
} else {
COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n",
result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
}
}
if (saveSol) {
if (extraBlock) {
delete [] useCostExtra;
}
CoinMemcpyN(saveSol, ncols, colsol);
delete [] saveSol;
}
for (i = 0; i < nsolve; i++) {
if (i) delete [] allsum[i];
delete [] aX[i];
delete [] aworkX[i];
}
delete [] thetaX;
delete [] djX;
delete [] bX;
delete [] vX;
delete [] aX;
delete [] aworkX;
delete [] allsum;
delete [] cost;
#ifdef FOUR_GOES
delete [] pi2 ;
delete [] rowsol2 ;
#endif
for (i = 0; i < HISTORY + 1; i++) {
delete [] history[i];
}
delete [] statusSave;
/* do original costs objvalue*/
result.objval = 0.0;
for (i = 0; i < ncols; i++) {
result.objval += colsol[i] * origcost[i];
}
if (extraBlock) {
for (i = 0; i < extraBlock; i++) {
int irow = rowExtra[i];
rowupper[irow] = saveExtra[i];
}
delete [] rowExtra;
delete [] solExtra;
delete [] elemExtra;
delete [] upperExtra;
delete [] costExtra;
delete [] saveExtra;
}
result.iteration = iter;
result.objval -= saveOffset;
result.weighted = result.objval + weight * result.sumSquared;
return result;
}