IdiSolve.cpp

// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

#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 * rowsol , double * colsol,
              double * pi, double * djs, const double * origcost , const double * rowlower,
              double * rowupper, const double * lower,
              const double * upper, const double * elemnt, 
                    const int * row, const CoinBigIndex * columnStart,
                    const int * length, double * lambda,
                    int maxIts,double mu,double drop,
                    double maxmin, double offset,
                    int strategy,double djTol,double djExit,double djFlag)
{
  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 * solExtra=NULL;
  double * elemExtra=NULL;
  double * upperExtra=NULL;
  double * costExtra=NULL;
  double * useCostExtra=NULL;
  double * saveExtra=NULL;
  double * 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 * thetaX;
  double * djX;
  double * bX;
  double * vX;
  double ** aX;
  double **aworkX;
  double ** allsum;
  double * saveSol=0;
  const double * 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;
  }
  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;
    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) {
          printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
                 i,smaller);
          abort();
        }
        difference=rowupper[i]-rowlower[i];
        difference = min(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;
      }
    }
  }
  for (i=0;i<nrows;i++) {
    offset+=lambda[i]*rowsol[i];
  }
  if ((strategy&256)!=0) {
    /* save best solution */
    saveSol=new double[ncols];
    memcpy(saveSol,colsol,ncols*sizeof(double));
    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,djtot;
    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 * x=history[HISTORY-1];
        for (i=HISTORY-1;i>0;i--) {
          history[i]=history[i-1];
        }
        history[0]=x;
        memcpy(history[0],colsol,ncols*sizeof(double));
        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
      }
    }
    if ((iter % SAVEHISTORY)==0||doFull) {
      if ((strategy&16)==0) {
        double * x=history[HISTORY-1];
        for (i=HISTORY-1;i>0;i--) {
          history[i]=history[i-1];
        }
        history[0]=x;
        memcpy(history[0],colsol,ncols*sizeof(double));
        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
      }
    }
    /* start full try */
    if ((iter % EVERY)==0||doFull) {
      // for next pass
      direction = - direction;
      // randomize. 
      // The cast is to avoid gcc compiler warning
      int kcol = (int)(ncols*CoinDrand48());
      if (kcol==ncols)
        kcol=ncols-1;
      if (direction>0) {
        start[0]=kcol;stop[0]=ncols;
        start[1]=0;stop[1]=kcol;
      } else {
        start[0]=kcol;stop[0]=-1;
        start[1]=ncols-1;stop[1]=kcol;
      }
      int itry=0;
      /*if ((strategy&16)==0) {
        double * x=history[HISTORY-1];
        for (i=HISTORY-1;i>0;i--) {
          history[i]=history[i-1];
        }
        history[0]=x;
        memcpy(history[0],colsol,ncols*sizeof(double));
        memcpy(history[0]+ncols,solExtra,extraBlock*sizeof(double));
        }*/
      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));
          djtot=0.0;
          {
            double * theta=thetaX;
            double * dj=djX;
            double * b=bX;
            double ** a=aX;
            double ** awork=aworkX;
            double * 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;
            memcpy(history[HISTORY-1],colsol,ncols*sizeof(double));
            memcpy(history[HISTORY-1]+ncols,solExtra,extraBlock*sizeof(double));
#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)) {
          // cast is to avoid gcc compiler
          // warning: initialization to 'int' from 'double'
                      int ii=(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  = max (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) {
            printf("%d %g better than %g\n",iter,
                   result.weighted*maxmin-useOffset,bestSol*maxmin-useOffset);
            bestSol=result.weighted;
            memcpy(saveSol,colsol,ncols*sizeof(double));
          }
        }
#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;
            memcpy(rowsol,row2,nrows*sizeof(double));
            memcpy(pi,pi2,nrows*sizeof(double));
            memcpy(colsol,sol2,ncols*sizeof(double));
            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);
            }
              memcpy(colsol,history[0],ncols*sizeof(double));
              memcpy(solExtra,history[0]+ncols,extraBlock*sizeof(double));
              good=1;
            }
          } else if ((strategy&3)==2) {
            if (weightedObj>lastObj+0.1*maxDj) {
              memcpy(colsol,history[0],ncols*sizeof(double));
              memcpy(solExtra,history[0]+ncols,extraBlock*sizeof(double));
              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);
        }
      }
      memcpy(statusWork,statusSave,ncols*sizeof(char));
      nflagged=0;
    }
    nChange=0;
    doFull=0;
    maxDj=0.0;
    // go through forwards or backwards and starting at odd places
    int itry;
    for (itry=0;itry<2;itry++) {
      int icol=start[itry];
      int istop= stop[itry];
      for (;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) {
              double valuep,thetap;
              nChange++;
              if (djval>maxDj) maxDj=djval;
              /*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;
              if ((strategy&4)!=0) {
                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++;
                }
              }
              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;
              objvalue += cost[icol]*value2;
              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;
                nflagged++;
              }
            }
          }
        }
      }
    }
    if (extraBlock) {
      for (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=min(theta,upperExtra[i]-solExtra[i]);
          } else {
            value2=max(theta,-solExtra[i]);
          }
          solExtra[i] += value2;
          rowsol[irow]+=element*value2;
          value=rowsol[irow];
          pi[irow]=-2.0*weight*value;
        }
      }
    }
    if ((iter%10)==2) {
      for (i=DJTEST-1;i>0;i--) {
        djSave[i]=djSave[i-1];
      }
      djSave[0]=maxDj;
      largestDj=max(largestDj,maxDj);
      smallestDj=min(smallestDj,maxDj);
      for (i=DJTEST-1;i>0;i--) {
        maxDj+=djSave[i];
      }
      maxDj = maxDj/(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) {
    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;
      memcpy(saveSol,colsol,ncols*sizeof(double));
    } else {
      printf("restoring previous - now %g best %g\n",
             result.weighted*maxmin-useOffset,bestSol*maxmin-useOffset);
    }
  }
  if (saveSol) {
    if (extraBlock) {
      delete [] useCostExtra;
    }
    memcpy(colsol,saveSol,ncols*sizeof(double));
    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;
  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;
}

---