// 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 , 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;
    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)
            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)
      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];
    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  = 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) {
	    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=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 (i=DJTEST-1;i>0;i--) {
	djSave[i]=djSave[i-1];
      }
      djSave[0]=maxDj;
      largestDj=CoinMax(largestDj,maxDj);
      smallestDj=CoinMin(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;
}