source: trunk/Clp/src/ClpSimplexNonlinear.cpp @ 1141

Last change on this file since 1141 was 1141, checked in by forrest, 12 years ago

for threaded random numbers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 122.5 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <math.h>
7#include "CoinHelperFunctions.hpp"
8#include "ClpHelperFunctions.hpp"
9#include "ClpSimplexNonlinear.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpNonLinearCost.hpp"
12#include "ClpLinearObjective.hpp"
13#include "ClpConstraint.hpp"
14#include "ClpQuadraticObjective.hpp"
15#include "CoinPackedMatrix.hpp"
16#include "CoinIndexedVector.hpp"
17#include "ClpPrimalColumnPivot.hpp"
18#include "ClpMessage.hpp"
19#include "ClpEventHandler.hpp"
20#include <cfloat>
21#include <cassert>
22#include <string>
23#include <stdio.h>
24#include <iostream>
25#ifndef NDEBUG
26#define CLP_DEBUG 1
27#endif
28// primal
29int ClpSimplexNonlinear::primal ()
30{
31
32  int ifValuesPass=1;
33  algorithm_ = +3;
34
35  // save data
36  ClpDataSave data = saveData();
37  matrix_->refresh(this); // make sure matrix okay
38
39  // Save objective
40  ClpObjective * saveObjective=NULL;
41  if (objective_->type()>1) {
42    // expand to full if quadratic
43#ifndef NO_RTTI
44    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
45#else
46    ClpQuadraticObjective * quadraticObj = NULL;
47    if (objective_->type()==2)
48      quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
49#endif
50    // for moment only if no scaling
51    // May be faster if switched off - but can't see why
52    if (!quadraticObj->fullMatrix()&&!rowScale_&&objectiveScale_==1.0) {
53      saveObjective = objective_;
54      objective_=new ClpQuadraticObjective(*quadraticObj,1);
55    }
56  }
57  double bestObjectiveWhenFlagged=COIN_DBL_MAX;
58  int pivotMode=15;
59  //pivotMode=20;
60
61  // initialize - maybe values pass and algorithm_ is +1
62  // true does something (? perturbs)
63  if (!startup(true)) {
64   
65    // Set average theta
66    nonLinearCost_->setAverageTheta(1.0e3);
67    int lastCleaned=0; // last time objective or bounds cleaned up
68   
69    // Say no pivot has occurred (for steepest edge and updates)
70    pivotRow_=-2;
71   
72    // This says whether to restore things etc
73    int factorType=0;
74    // Start check for cycles
75    progress_->startCheck();
76    /*
77      Status of problem:
78      0 - optimal
79      1 - infeasible
80      2 - unbounded
81      -1 - iterating
82      -2 - factorization wanted
83      -3 - redo checking without factorization
84      -4 - looks infeasible
85      -5 - looks unbounded
86    */
87    while (problemStatus_<0) {
88      int iRow,iColumn;
89      // clear
90      for (iRow=0;iRow<4;iRow++) {
91        rowArray_[iRow]->clear();
92      }   
93     
94      for (iColumn=0;iColumn<2;iColumn++) {
95        columnArray_[iColumn]->clear();
96      }   
97     
98      // give matrix (and model costs and bounds a chance to be
99      // refreshed (normally null)
100      matrix_->refresh(this);
101      // If getting nowhere - why not give it a kick
102      // If we have done no iterations - special
103      if (lastGoodIteration_==numberIterations_&&factorType)
104        factorType=3;
105     
106      // may factorize, checks if problem finished
107      if (objective_->type()>1&&lastFlaggedIteration_>=0&&
108          numberIterations_>lastFlaggedIteration_+507) {
109        unflag();
110        lastFlaggedIteration_=numberIterations_;
111        if (pivotMode>=10) {
112          pivotMode--;
113#ifdef CLP_DEBUG
114          if (handler_->logLevel()&32) 
115            printf("pivot mode now %d\n",pivotMode);
116#endif
117          if (pivotMode==9) 
118            pivotMode=0;        // switch off fast attempt
119        }
120      }
121      statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,
122                              bestObjectiveWhenFlagged);
123     
124      // Say good factorization
125      factorType=1;
126     
127      // Say no pivot has occurred (for steepest edge and updates)
128      pivotRow_=-2;
129
130      // exit if victory declared
131      if (problemStatus_>=0)
132        break;
133     
134      // test for maximum iterations
135      if (hitMaximumIterations()||(ifValuesPass==2&&firstFree_<0)) {
136        problemStatus_=3;
137        break;
138      }
139
140      if (firstFree_<0) {
141        if (ifValuesPass) {
142          // end of values pass
143          ifValuesPass=0;
144          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
145          if (status>=0) {
146            problemStatus_=5;
147            secondaryStatus_=ClpEventHandler::endOfValuesPass;
148            break;
149          }
150        }
151      }
152      // Check event
153      {
154        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
155        if (status>=0) {
156          problemStatus_=5;
157          secondaryStatus_=ClpEventHandler::endOfFactorization;
158          break;
159        }
160      }
161      // Iterate
162      whileIterating(pivotMode);
163    }
164  }
165  // if infeasible get real values
166  if (problemStatus_==1) {
167    infeasibilityCost_=0.0;
168    createRim(1+4);
169    nonLinearCost_->checkInfeasibilities(0.0);
170    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
171    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
172    // and get good feasible duals
173    computeDuals(NULL);
174  }
175  // correct objective value
176  if (numberColumns_)
177    objectiveValue_ = nonLinearCost_->feasibleCost()+objective_->nonlinearOffset();
178  objectiveValue_ /= (objectiveScale_*rhsScale_);
179  // clean up
180  unflag();
181  finish();
182  restoreData(data);
183  // restore objective if full
184  if (saveObjective) {
185    delete objective_;
186    objective_=saveObjective;
187  }
188  return problemStatus_;
189}
190/*  Refactorizes if necessary
191    Checks if finished.  Updates status.
192    lastCleaned refers to iteration at which some objective/feasibility
193    cleaning too place.
194   
195    type - 0 initial so set up save arrays etc
196    - 1 normal -if good update save
197    - 2 restoring from saved
198*/
199void 
200ClpSimplexNonlinear::statusOfProblemInPrimal(int & lastCleaned, int type,
201                                                 ClpSimplexProgress * progress,
202                                                 bool doFactorization,
203                                                 double & bestObjectiveWhenFlagged)
204{
205  int dummy; // for use in generalExpanded
206  if (type==2) {
207    // trouble - restore solution
208    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
209    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
210           numberRows_*sizeof(double));
211    memcpy(columnActivityWork_,savedSolution_ ,
212           numberColumns_*sizeof(double));
213    // restore extra stuff
214    matrix_->generalExpanded(this,6,dummy);
215    forceFactorization_=1; // a bit drastic but ..
216    pivotRow_=-1; // say no weights update
217    changeMade_++; // say change made
218  }
219  int saveThreshold = factorization_->sparseThreshold();
220  int tentativeStatus = problemStatus_;
221  int numberThrownOut=1; // to loop round on bad factorization in values pass
222  while (numberThrownOut) {
223    if (problemStatus_>-3||problemStatus_==-4) {
224      // factorize
225      // later on we will need to recover from singularities
226      // also we could skip if first time
227      // do weights
228      // This may save pivotRow_ for use
229      if (doFactorization)
230        primalColumnPivot_->saveWeights(this,1);
231     
232      if (type&&doFactorization) {
233        // is factorization okay?
234        int factorStatus = internalFactorize(1);
235        if (factorStatus) {
236          if (type!=1||largestPrimalError_>1.0e3
237              ||largestDualError_>1.0e3) {
238            // was ||largestDualError_>1.0e3||objective_->type()>1) {
239            // switch off dense
240            int saveDense = factorization_->denseThreshold();
241            factorization_->setDenseThreshold(0);
242            // make sure will do safe factorization
243            pivotVariable_[0]=-1;
244            internalFactorize(2);
245            factorization_->setDenseThreshold(saveDense);
246            // Go to safe
247            factorization_->pivotTolerance(0.99);
248            // restore extra stuff
249            matrix_->generalExpanded(this,6,dummy);
250          } else {
251            // no - restore previous basis
252            memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
253            memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
254                   numberRows_*sizeof(double));
255            memcpy(columnActivityWork_,savedSolution_ ,
256                   numberColumns_*sizeof(double));
257            // restore extra stuff
258            matrix_->generalExpanded(this,6,dummy);
259            matrix_->generalExpanded(this,5,dummy);
260            forceFactorization_=1; // a bit drastic but ..
261            type = 2;
262            // Go to safe
263            factorization_->pivotTolerance(0.99);
264            if (internalFactorize(1)!=0)
265               largestPrimalError_=1.0e4; // force other type
266          }
267          // flag incoming
268          if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
269            setFlagged(sequenceIn_);
270            saveStatus_[sequenceIn_]=status_[sequenceIn_];
271          }
272          changeMade_++; // say change made
273        }
274      }
275      if (problemStatus_!=-4)
276        problemStatus_=-3;
277    }
278    // at this stage status is -3 or -5 if looks unbounded
279    // get primal and dual solutions
280    // put back original costs and then check
281    createRim(4);
282    // May need to do more if column generation
283    dummy=4;
284    matrix_->generalExpanded(this,9,dummy);
285    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
286    if (numberThrownOut) {
287      problemStatus_=tentativeStatus;
288      doFactorization=true;
289    }
290  }
291  // Double check reduced costs if no action
292  if (progress->lastIterationNumber(0)==numberIterations_) {
293    if (primalColumnPivot_->looksOptimal()) {
294      numberDualInfeasibilities_ = 0;
295      sumDualInfeasibilities_ = 0.0;
296    }
297  }
298  // Check if looping
299  int loop;
300  if (type!=2) 
301    loop = progress->looping();
302  else
303    loop=-1;
304  if (loop>=0) {
305    if (!problemStatus_) {
306      // declaring victory
307      numberPrimalInfeasibilities_ = 0;
308      sumPrimalInfeasibilities_ = 0.0;
309    } else {
310      problemStatus_ = loop; //exit if in loop
311      problemStatus_ = 10; // instead - try other algorithm
312    }
313    problemStatus_ = 10; // instead - try other algorithm
314    return ;
315  } else if (loop<-1) {
316    // Is it time for drastic measures
317    if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
318        progress->oddState()<10&&progress->oddState()>=0) {
319      progress->newOddState();
320      nonLinearCost_->zapCosts();
321    }
322    // something may have changed
323    gutsOfSolution(NULL,NULL,true);
324  }
325  // If progress then reset costs
326  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
327    createRim(4,false); // costs back
328    delete nonLinearCost_;
329    nonLinearCost_ = new ClpNonLinearCost(this);
330    progress->endOddState();
331    gutsOfSolution(NULL,NULL,true);
332  }
333  // Flag to say whether to go to dual to clean up
334  bool goToDual=false;
335  // really for free variables in
336  //if((progressFlag_&2)!=0)
337  //problemStatus_=-1;;
338  progressFlag_ = 0; //reset progress flag
339
340  handler_->message(CLP_SIMPLEX_STATUS,messages_)
341    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
342  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
343    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
344  handler_->printing(sumDualInfeasibilities_>0.0)
345    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
346  handler_->printing(numberDualInfeasibilitiesWithoutFree_
347                     <numberDualInfeasibilities_)
348                       <<numberDualInfeasibilitiesWithoutFree_;
349  handler_->message()<<CoinMessageEol;
350  if (!primalFeasible()) {
351    nonLinearCost_->checkInfeasibilities(primalTolerance_);
352    gutsOfSolution(NULL,NULL,true);
353    nonLinearCost_->checkInfeasibilities(primalTolerance_);
354  }
355  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
356  if (trueInfeasibility>1.0) {
357    // If infeasibility going up may change weights
358    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
359    if(progress->lastInfeasibility()<testValue) {
360      if (infeasibilityCost_<1.0e14) {
361        infeasibilityCost_ *= 1.5;
362        if (handler_->logLevel()==63)
363          printf("increasing weight to %g\n",infeasibilityCost_);
364        gutsOfSolution(NULL,NULL,true);
365      }
366    }
367  }
368  // we may wish to say it is optimal even if infeasible
369  bool alwaysOptimal = (specialOptions_&1)!=0;
370  // give code benefit of doubt
371  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
372      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
373    // say optimal (with these bounds etc)
374    numberDualInfeasibilities_ = 0;
375    sumDualInfeasibilities_ = 0.0;
376    numberPrimalInfeasibilities_ = 0;
377    sumPrimalInfeasibilities_ = 0.0;
378  }
379  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
380  if (dualFeasible()||problemStatus_==-4) {
381    // see if extra helps
382    if (nonLinearCost_->numberInfeasibilities()&&
383         (nonLinearCost_->sumInfeasibilities()>1.0e-3||sumOfRelaxedPrimalInfeasibilities_)
384        &&!alwaysOptimal) {
385      //may need infeasiblity cost changed
386      // we can see if we can construct a ray
387      // make up a new objective
388      double saveWeight = infeasibilityCost_;
389      // save nonlinear cost as we are going to switch off costs
390      ClpNonLinearCost * nonLinear = nonLinearCost_;
391      // do twice to make sure Primal solution has settled
392      // put non-basics to bounds in case tolerance moved
393      // put back original costs
394      createRim(4);
395      nonLinearCost_->checkInfeasibilities(0.0);
396      gutsOfSolution(NULL,NULL,true);
397
398      infeasibilityCost_=1.0e100;
399      // put back original costs
400      createRim(4);
401      nonLinearCost_->checkInfeasibilities(primalTolerance_);
402      // may have fixed infeasibilities - double check
403      if (nonLinearCost_->numberInfeasibilities()==0) {
404        // carry on
405        problemStatus_ = -1;
406        infeasibilityCost_=saveWeight;
407        nonLinearCost_->checkInfeasibilities(primalTolerance_);
408      } else {
409        nonLinearCost_=NULL;
410        // scale
411        int i;
412        for (i=0;i<numberRows_+numberColumns_;i++) 
413          cost_[i] *= 1.0e-95;
414        gutsOfSolution(NULL,NULL,false);
415        nonLinearCost_=nonLinear;
416        infeasibilityCost_=saveWeight;
417        if ((infeasibilityCost_>=1.0e18||
418             numberDualInfeasibilities_==0)&&perturbation_==101) {
419          goToDual=unPerturb(); // stop any further perturbation
420          if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
421            goToDual=false;
422          nonLinearCost_->checkInfeasibilities(primalTolerance_);
423          numberDualInfeasibilities_=1; // carry on
424          problemStatus_=-1;
425        }
426        if (infeasibilityCost_>=1.0e20||
427            numberDualInfeasibilities_==0) {
428          // we are infeasible - use as ray
429          delete [] ray_;
430          ray_ = new double [numberRows_];
431          memcpy(ray_,dual_,numberRows_*sizeof(double));
432          // and get feasible duals
433          infeasibilityCost_=0.0;
434          createRim(4);
435          nonLinearCost_->checkInfeasibilities(primalTolerance_);
436          gutsOfSolution(NULL,NULL,true);
437          // so will exit
438          infeasibilityCost_=1.0e30;
439          // reset infeasibilities
440          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
441          numberPrimalInfeasibilities_=
442            nonLinearCost_->numberInfeasibilities();
443        }
444        if (infeasibilityCost_<1.0e20) {
445          infeasibilityCost_ *= 5.0;
446          changeMade_++; // say change made
447          unflag();
448          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
449            <<infeasibilityCost_
450            <<CoinMessageEol;
451          // put back original costs and then check
452          createRim(4);
453          nonLinearCost_->checkInfeasibilities(0.0);
454          gutsOfSolution(NULL,NULL,true);
455          problemStatus_=-1; //continue
456          goToDual=false;
457        } else {
458          // say infeasible
459          problemStatus_ = 1;
460        }
461      }
462    } else {
463      // may be optimal
464      if (perturbation_==101) {
465        goToDual=unPerturb(); // stop any further perturbation
466        lastCleaned=-1; // carry on
467      }
468      bool unflagged = unflag()!=0;
469      if ( lastCleaned!=numberIterations_||unflagged) {
470        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
471          <<primalTolerance_
472          <<CoinMessageEol;
473        double current = nonLinearCost_->feasibleReportCost();
474        if (numberTimesOptimal_<4) {
475          if (bestObjectiveWhenFlagged<=current) {
476            numberTimesOptimal_++;
477#ifdef CLP_DEBUG
478            if (handler_->logLevel()&32) 
479              printf("%d times optimal, current %g best %g\n",numberTimesOptimal_,
480                     current,bestObjectiveWhenFlagged);
481#endif
482          } else {
483            bestObjectiveWhenFlagged=current; 
484          }
485          changeMade_++; // say change made
486          if (numberTimesOptimal_==1) {
487            // better to have small tolerance even if slower
488            factorization_->zeroTolerance(1.0e-15);
489          }
490          lastCleaned=numberIterations_;
491          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
492            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
493              <<CoinMessageEol;
494          double oldTolerance = primalTolerance_;
495          primalTolerance_=dblParam_[ClpPrimalTolerance];
496          // put back original costs and then check
497          createRim(4);
498          nonLinearCost_->checkInfeasibilities(oldTolerance);
499          gutsOfSolution(NULL,NULL,true);
500          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
501              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
502            // say optimal (with these bounds etc)
503            numberDualInfeasibilities_ = 0;
504            sumDualInfeasibilities_ = 0.0;
505            numberPrimalInfeasibilities_ = 0;
506            sumPrimalInfeasibilities_ = 0.0;
507          }
508          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
509            problemStatus_=0;
510          else
511            problemStatus_ = -1;
512        } else {
513          problemStatus_=0; // optimal
514          if (lastCleaned<numberIterations_) {
515            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
516              <<CoinMessageEol;
517          }
518        }
519      } else {
520        problemStatus_=0; // optimal
521      }
522    }
523  } else {
524    // see if looks unbounded
525    if (problemStatus_==-5) {
526      if (nonLinearCost_->numberInfeasibilities()) {
527        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
528          // back off weight
529          infeasibilityCost_ = 1.0e13;
530          unPerturb(); // stop any further perturbation
531        }
532        //we need infeasiblity cost changed
533        if (infeasibilityCost_<1.0e20) {
534          infeasibilityCost_ *= 5.0;
535          changeMade_++; // say change made
536          unflag();
537          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
538            <<infeasibilityCost_
539            <<CoinMessageEol;
540          // put back original costs and then check
541          createRim(4);
542          gutsOfSolution(NULL,NULL,true);
543          problemStatus_=-1; //continue
544        } else {
545          // say unbounded
546          problemStatus_ = 2;
547        }
548      } else {
549        // say unbounded
550        problemStatus_ = 2;
551      } 
552    } else {
553      if(type==3&&problemStatus_!=-5)
554        unflag(); // odd
555      // carry on
556      problemStatus_ = -1;
557    }
558  }
559  // save extra stuff
560  matrix_->generalExpanded(this,5,dummy);
561  if (type==0||type==1) {
562    if (type!=1||!saveStatus_) {
563      // create save arrays
564      delete [] saveStatus_;
565      delete [] savedSolution_;
566      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
567      savedSolution_ = new double [numberRows_+numberColumns_];
568    }
569    // save arrays
570    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
571    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
572           numberRows_*sizeof(double));
573    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
574  }
575  if (doFactorization) {
576    // restore weights (if saved) - also recompute infeasibility list
577    if (tentativeStatus>-3) 
578      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
579    else
580      primalColumnPivot_->saveWeights(this,3);
581    if (saveThreshold) {
582      // use default at present
583      factorization_->sparseThreshold(0);
584      factorization_->goSparse();
585    }
586  }
587  if (problemStatus_<0&&!changeMade_) {
588    problemStatus_=4; // unknown
589  }
590  lastGoodIteration_ = numberIterations_;
591  //if (goToDual)
592  //problemStatus_=10; // try dual
593  // Allow matrices to be sorted etc
594  int fake=-999; // signal sort
595  matrix_->correctSequence(this,fake,fake);
596}
597/*
598  Reasons to come out:
599  -1 iterations etc
600  -2 inaccuracy
601  -3 slight inaccuracy (and done iterations)
602  -4 end of values pass and done iterations
603  +0 looks optimal (might be infeasible - but we will investigate)
604  +2 looks unbounded
605  +3 max iterations
606*/
607int
608ClpSimplexNonlinear::whileIterating(int & pivotMode)
609{
610  // Say if values pass
611  //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
612  int ifValuesPass=1;
613  int returnCode=-1;
614  // status stays at -1 while iterating, >=0 finished, -2 to invert
615  // status -3 to go to top without an invert
616  int numberInterior=0;
617  int nextUnflag=10;
618  int nextUnflagIteration=numberIterations_+10;
619  // get two arrays
620  double * array1 = new double[2*(numberRows_+numberColumns_)];
621  double solutionError=-1.0;
622  while (problemStatus_==-1) {
623    int result;
624    rowArray_[1]->clear();
625    //#define CLP_DEBUG
626#if CLP_DEBUG > 1
627    rowArray_[0]->checkClear();
628    rowArray_[1]->checkClear();
629    rowArray_[2]->checkClear();
630    rowArray_[3]->checkClear();
631    columnArray_[0]->checkClear();
632#endif
633    if (numberInterior>=5) {
634      // this can go when we have better minimization
635      if (pivotMode<10)
636        pivotMode=1;
637      unflag();
638#ifdef CLP_DEBUG
639          if (handler_->logLevel()&32) 
640            printf("interior unflag\n");
641#endif
642      numberInterior=0;
643      nextUnflag=10;
644      nextUnflagIteration=numberIterations_+10;
645    } else {
646      if (numberInterior>nextUnflag&&
647          numberIterations_>nextUnflagIteration) {
648        nextUnflagIteration=numberIterations_+10;
649        nextUnflag += 10;
650        unflag();
651#ifdef CLP_DEBUG
652        if (handler_->logLevel()&32) 
653          printf("unflagging as interior\n");
654#endif
655      }
656    } 
657    pivotRow_=-1;
658    result = pivotColumn(rowArray_[3],rowArray_[0],
659                         columnArray_[0],rowArray_[1],pivotMode,solutionError,
660                         array1);
661    if (result) {
662      if (result==3) 
663        break; // null vector not accurate
664#ifdef CLP_DEBUG
665      if (handler_->logLevel()&32) {
666        double currentObj;
667        double thetaObj;
668        double predictedObj;
669        objective_->stepLength(this,solution_,solution_,0.0,
670                               currentObj,predictedObj,thetaObj);
671        printf("obj %g after interior move\n",currentObj);
672      }
673#endif
674      // just move and try again
675      if (pivotMode<10) {
676        pivotMode=result-1;
677        numberInterior++;
678      }
679      continue;
680    } else {
681      if (pivotMode<10) {
682        if (theta_>0.001)
683          pivotMode=0;
684        else if (pivotMode==2)
685          pivotMode=1;
686      }
687      numberInterior=0;
688      nextUnflag=10;
689      nextUnflagIteration=numberIterations_+10;
690    }
691    sequenceOut_=-1;
692    rowArray_[1]->clear();
693    if (sequenceIn_>=0) {
694      // we found a pivot column
695      assert (!flagged(sequenceIn_));
696#ifdef CLP_DEBUG
697      if ((handler_->logLevel()&32)) {
698        char x = isColumn(sequenceIn_) ? 'C' :'R';
699        std::cout<<"pivot column "<<
700          x<<sequenceWithin(sequenceIn_)<<std::endl;
701      }
702#endif
703      // do second half of iteration
704      if (pivotRow_<0&&theta_<1.0e-8) {
705        assert (sequenceIn_>=0);
706        returnCode = pivotResult(ifValuesPass);
707      } else {
708        // specialized code
709        returnCode = pivotNonlinearResult();
710        //printf("odd pivrow %d\n",sequenceOut_);
711        if (sequenceOut_>=0&&theta_<1.0e-5) {
712          if (getStatus(sequenceOut_)!=isFixed) {
713            if (getStatus(sequenceOut_)==atUpperBound)
714              solution_[sequenceOut_]=upper_[sequenceOut_];
715            else if (getStatus(sequenceOut_)==atLowerBound)
716              solution_[sequenceOut_]=lower_[sequenceOut_];
717            setFlagged(sequenceOut_);
718          }
719        }
720      }
721      if (returnCode<-1&&returnCode>-5) {
722        problemStatus_=-2; //
723      } else if (returnCode==-5) {
724        // something flagged - continue;
725      } else if (returnCode==2) {
726        problemStatus_=-5; // looks unbounded
727      } else if (returnCode==4) {
728        problemStatus_=-2; // looks unbounded but has iterated
729      } else if (returnCode!=-1) {
730        assert(returnCode==3);
731        problemStatus_=3;
732      }
733    } else {
734      // no pivot column
735#ifdef CLP_DEBUG
736      if (handler_->logLevel()&32)
737        printf("** no column pivot\n");
738#endif
739      if (pivotMode<10) {
740        // looks optimal
741        primalColumnPivot_->setLooksOptimal(true);
742      } else {
743        pivotMode--;
744#ifdef CLP_DEBUG
745        if (handler_->logLevel()&32) 
746          printf("pivot mode now %d\n",pivotMode);
747#endif
748        if (pivotMode==9) 
749          pivotMode=0;  // switch off fast attempt
750        unflag();
751      }
752      if (nonLinearCost_->numberInfeasibilities())
753        problemStatus_=-4; // might be infeasible
754      returnCode=0;
755      break;
756    }
757  }
758  delete [] array1;
759  return returnCode;
760}
761// Creates direction vector
762void
763ClpSimplexNonlinear::directionVector (CoinIndexedVector * vectorArray,
764                                      CoinIndexedVector * spare1, CoinIndexedVector * spare2,
765                                      int pivotMode2,
766                                      double & normFlagged,double & normUnflagged,
767                                      int & numberNonBasic)
768{
769#if CLP_DEBUG > 1
770  vectorArray->checkClear();
771  spare1->checkClear();
772  spare2->checkClear();
773#endif
774  double *array = vectorArray->denseVector();
775  int * index = vectorArray->getIndices();
776  int number=0;
777  sequenceIn_=-1;
778  normFlagged=0.0;
779  normUnflagged=1.0;
780  double dualTolerance2 = CoinMin(1.0e-8,1.0e-2*dualTolerance_);
781  double dualTolerance3 = CoinMin(1.0e-2,1.0e3*dualTolerance_);
782  if (!numberNonBasic) {
783    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
784    //printf("infeasible\n");
785    if (!pivotMode2||pivotMode2>=10) {
786      normUnflagged=0.0;
787      double bestDj=0.0;
788      double bestSuper=0.0;
789      double sumSuper=0.0;
790      sequenceIn_=-1;
791      int nSuper=0;
792      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
793        array[iSequence]=0.0;
794        if (flagged(iSequence)) {
795          // accumulate  norm
796          switch(getStatus(iSequence)) {
797           
798          case basic:
799          case ClpSimplex::isFixed:
800            break;
801          case atUpperBound:
802            if (dj_[iSequence]>dualTolerance3) 
803              normFlagged += dj_[iSequence]*dj_[iSequence];
804            break;
805          case atLowerBound:
806            if (dj_[iSequence]<-dualTolerance3) 
807              normFlagged += dj_[iSequence]*dj_[iSequence];
808            break;
809          case isFree:
810          case superBasic:
811            if (fabs(dj_[iSequence])>dualTolerance3) 
812              normFlagged += dj_[iSequence]*dj_[iSequence];
813            break;
814          }
815          continue;
816        }
817        switch(getStatus(iSequence)) {
818         
819        case basic:
820        case ClpSimplex::isFixed:
821          break;
822        case atUpperBound:
823          if (dj_[iSequence]>dualTolerance_) {
824            if (dj_[iSequence]>dualTolerance3) 
825              normUnflagged += dj_[iSequence]*dj_[iSequence];
826            if (pivotMode2<10) {
827              array[iSequence]=-dj_[iSequence];
828              index[number++]=iSequence;
829            } else {
830              if (dj_[iSequence]>bestDj) {
831                bestDj=dj_[iSequence];
832                sequenceIn_=iSequence;
833              }
834            }
835          }
836          break;
837        case atLowerBound:
838          if (dj_[iSequence]<-dualTolerance_) {
839            if (dj_[iSequence]<-dualTolerance3) 
840              normUnflagged += dj_[iSequence]*dj_[iSequence];
841            if (pivotMode2<10) {
842              array[iSequence]=-dj_[iSequence];
843              index[number++]=iSequence;
844            } else {
845              if (-dj_[iSequence]>bestDj) {
846                bestDj=-dj_[iSequence];
847                sequenceIn_=iSequence;
848              }
849            }
850          }
851          break;
852        case isFree:
853        case superBasic:
854          //#define ALLSUPER
855#define NOSUPER
856#ifndef ALLSUPER
857          if (fabs(dj_[iSequence])>dualTolerance_) {
858            if (fabs(dj_[iSequence])>dualTolerance3) 
859              normUnflagged += dj_[iSequence]*dj_[iSequence];
860            nSuper++;
861            bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
862            sumSuper += fabs(dj_[iSequence]);
863          }
864          if (fabs(dj_[iSequence])>dualTolerance2) {
865            array[iSequence]=-dj_[iSequence];
866            index[number++]=iSequence;
867          }
868#else
869          array[iSequence]=-dj_[iSequence];
870          index[number++]=iSequence;
871          if (pivotMode2>=10)
872            bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
873#endif
874          break;
875        }
876      }
877#ifdef NOSUPER
878      // redo
879      bestSuper=sumSuper;
880      if(sequenceIn_>=0&&bestDj>bestSuper) {
881        int j;
882        // get rid of superbasics
883        for (j=0;j<number;j++) {
884          int iSequence = index[j];
885          array[iSequence]=0.0;
886        }
887        number=0;
888        array[sequenceIn_]=-dj_[sequenceIn_];
889        index[number++]=sequenceIn_;
890      } else {
891        sequenceIn_=-1;
892      }
893#else
894      if (pivotMode2>=10||!nSuper) {
895        bool takeBest=true;
896        if (pivotMode2==100&&nSuper>1)
897          takeBest=false;
898        if(sequenceIn_>=0&&takeBest) {
899          if (fabs(dj_[sequenceIn_])>bestSuper) {
900            array[sequenceIn_]=-dj_[sequenceIn_];
901            index[number++]=sequenceIn_;
902          } else {
903            sequenceIn_=-1;
904          }
905        } else {
906          sequenceIn_=-1;
907        }
908      }
909#endif
910#ifdef CLP_DEBUG
911      if (handler_->logLevel()&32) {
912        if (sequenceIn_>=0)
913          printf("%d superBasic, chosen %d - dj %g\n",nSuper,sequenceIn_,
914                 dj_[sequenceIn_]);
915        else
916          printf("%d superBasic - none chosen\n",nSuper);
917      }
918#endif
919    } else {
920      double bestDj=0.0;
921      double saveDj=0.0;
922      if (sequenceOut_>=0) {
923        saveDj=dj_[sequenceOut_];
924        dj_[sequenceOut_]=0.0;
925        switch(getStatus(sequenceOut_)) {
926         
927        case basic:
928          sequenceOut_=-1;
929        case ClpSimplex::isFixed:
930          break;
931        case atUpperBound:
932          if (dj_[sequenceOut_]>dualTolerance_) {
933#ifdef CLP_DEBUG
934            if (handler_->logLevel()&32) 
935              printf("after pivot out %d values %g %g %g, dj %g\n",
936                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
937                     upper_[sequenceOut_],dj_[sequenceOut_]);
938#endif
939          }
940          break;
941        case atLowerBound:
942          if (dj_[sequenceOut_]<-dualTolerance_) {
943#ifdef CLP_DEBUG
944            if (handler_->logLevel()&32) 
945              printf("after pivot out %d values %g %g %g, dj %g\n",
946                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
947                     upper_[sequenceOut_],dj_[sequenceOut_]);
948#endif
949          }
950          break;
951        case isFree:
952        case superBasic:
953          if (dj_[sequenceOut_]>dualTolerance_) {
954#ifdef CLP_DEBUG
955            if (handler_->logLevel()&32) 
956              printf("after pivot out %d values %g %g %g, dj %g\n",
957                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
958                     upper_[sequenceOut_],dj_[sequenceOut_]);
959#endif
960          } else if (dj_[sequenceOut_]<-dualTolerance_) {
961#ifdef CLP_DEBUG
962            if (handler_->logLevel()&32) 
963              printf("after pivot out %d values %g %g %g, dj %g\n",
964                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
965                     upper_[sequenceOut_],dj_[sequenceOut_]);
966#endif
967          }
968          break;
969        }
970      }
971      // Go for dj
972      pivotMode2=3;
973      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
974        array[iSequence]=0.0;
975        if (flagged(iSequence))
976          continue;
977        switch(getStatus(iSequence)) {
978         
979        case basic:
980        case ClpSimplex::isFixed:
981          break;
982        case atUpperBound:
983          if (dj_[iSequence]>dualTolerance_) {
984            double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
985            double merit = distance*dj_[iSequence];
986            if (pivotMode2==1)
987              merit *= 1.0e-20; // discourage
988            if (pivotMode2==3)
989              merit = fabs(dj_[iSequence]);
990            if (merit>bestDj) {
991              sequenceIn_=iSequence;
992              bestDj=merit;
993            }
994          }
995          break;
996        case atLowerBound:
997          if (dj_[iSequence]<-dualTolerance_) {
998            double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
999            double merit = -distance*dj_[iSequence];
1000            if (pivotMode2==1)
1001              merit *= 1.0e-20; // discourage
1002            if (pivotMode2==3)
1003              merit = fabs(dj_[iSequence]);
1004            if (merit>bestDj) {
1005              sequenceIn_=iSequence;
1006              bestDj=merit;
1007            }
1008          }
1009          break;
1010        case isFree:
1011        case superBasic:
1012          if (dj_[iSequence]>dualTolerance_) {
1013            double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
1014            distance = CoinMin(solution_[iSequence]-lower_[iSequence],
1015                           upper_[iSequence]-solution_[iSequence]);
1016            double merit = distance*dj_[iSequence];
1017            if (pivotMode2==1)
1018              merit=distance;
1019            if (pivotMode2==3)
1020              merit = fabs(dj_[iSequence]);
1021            if (merit>bestDj) {
1022              sequenceIn_=iSequence;
1023              bestDj=merit;
1024            }
1025          } else if (dj_[iSequence]<-dualTolerance_) {
1026            double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
1027            distance = CoinMin(solution_[iSequence]-lower_[iSequence],
1028                           upper_[iSequence]-solution_[iSequence]);
1029            double merit = -distance*dj_[iSequence];
1030            if (pivotMode2==1)
1031              merit=distance;
1032            if (pivotMode2==3)
1033              merit = fabs(dj_[iSequence]);
1034            if (merit>bestDj) {
1035              sequenceIn_=iSequence;
1036              bestDj=merit;
1037            }
1038          }
1039          break;
1040        }
1041      }
1042      if (sequenceOut_>=0) {
1043        dj_[sequenceOut_]=saveDj;
1044        sequenceOut_=-1;
1045      }
1046      if (sequenceIn_>=0) {
1047        array[sequenceIn_]=-dj_[sequenceIn_];
1048        index[number++]=sequenceIn_;
1049      }
1050    }
1051    numberNonBasic = number;
1052  } else {
1053    // compute norms
1054    normUnflagged=0.0;
1055    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
1056      if (flagged(iSequence)) {
1057        // accumulate  norm
1058        switch(getStatus(iSequence)) {
1059         
1060        case basic:
1061        case ClpSimplex::isFixed:
1062          break;
1063        case atUpperBound:
1064          if (dj_[iSequence]>dualTolerance_) 
1065            normFlagged += dj_[iSequence]*dj_[iSequence];
1066          break;
1067        case atLowerBound:
1068          if (dj_[iSequence]<-dualTolerance_) 
1069            normFlagged += dj_[iSequence]*dj_[iSequence];
1070          break;
1071        case isFree:
1072        case superBasic:
1073          if (fabs(dj_[iSequence])>dualTolerance_) 
1074            normFlagged += dj_[iSequence]*dj_[iSequence];
1075          break;
1076        }
1077      }
1078    }
1079    // re-use list
1080    number=0;
1081    int j;
1082    for (j=0;j<numberNonBasic;j++) {
1083      int iSequence = index[j];
1084      if (flagged(iSequence))
1085        continue;
1086      switch(getStatus(iSequence)) {
1087         
1088      case basic:
1089      case ClpSimplex::isFixed:
1090        continue; //abort();
1091        break;
1092      case atUpperBound:
1093        if (dj_[iSequence]>dualTolerance_) {
1094          number++;
1095          normUnflagged += dj_[iSequence]*dj_[iSequence];
1096        }
1097        break;
1098      case atLowerBound:
1099        if (dj_[iSequence]<-dualTolerance_) {
1100          number++;
1101          normUnflagged += dj_[iSequence]*dj_[iSequence];
1102        }
1103        break;
1104      case isFree:
1105      case superBasic:
1106        if (fabs(dj_[iSequence])>dualTolerance_) {
1107          number++;
1108          normUnflagged += dj_[iSequence]*dj_[iSequence];
1109        }
1110        break;
1111      }
1112      array[iSequence]=-dj_[iSequence];
1113    }
1114    // switch to large
1115    normUnflagged=1.0;
1116    if (!number) {
1117      for (j=0;j<numberNonBasic;j++) {
1118        int iSequence = index[j];
1119        array[iSequence]=0.0;
1120      }
1121      numberNonBasic=0;
1122    }
1123    number=numberNonBasic;
1124  }
1125  if (number) {
1126    int j;
1127    // Now do basic ones - how do I compensate for small basic infeasibilities?
1128    int iRow;
1129    for (iRow=0;iRow<numberRows_;iRow++) {
1130      int iPivot=pivotVariable_[iRow];
1131      double value=0.0;
1132      if (solution_[iPivot]>upper_[iPivot]) {
1133        value = upper_[iPivot]-solution_[iPivot];
1134      } else if (solution_[iPivot]<lower_[iPivot]) {
1135        value = lower_[iPivot]-solution_[iPivot];
1136      }
1137      //if (value)
1138      //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
1139      //     upper_[iPivot]);
1140      //value=0.0;
1141      value *= -1.0e0;
1142      if (value) {
1143        array[iPivot]=value;
1144        index[number++]=iPivot;
1145      }
1146    }
1147    double * array2=spare1->denseVector();
1148    int * index2 = spare1->getIndices();
1149    int number2=0;
1150    times(-1.0,array,array2);
1151    array = array+numberColumns_;
1152    for (iRow=0;iRow<numberRows_;iRow++) {
1153      double value = array2[iRow] + array[iRow];
1154      if (value) {
1155        array2[iRow]=value;
1156        index2[number2++]=iRow;
1157      } else {
1158        array2[iRow]=0.0;
1159      }
1160    }
1161    array -= numberColumns_;
1162    spare1->setNumElements(number2);
1163    // Ftran
1164    factorization_->updateColumn(spare2,spare1);
1165    number2=spare1->getNumElements();
1166    for (j=0;j<number2;j++) {
1167      int iSequence = index2[j];
1168      double value = array2[iSequence];
1169      array2[iSequence]=0.0;
1170      if (value) {
1171        int iPivot=pivotVariable_[iSequence];
1172        double oldValue = array[iPivot];
1173        if (!oldValue) {
1174          array[iPivot] = value;
1175          index[number++]=iPivot;
1176        } else {
1177          // something already there
1178          array[iPivot] = value+oldValue;
1179        }
1180      }
1181    }
1182    spare1->setNumElements(0);
1183  }
1184  vectorArray->setNumElements(number);
1185}
1186#define MINTYPE 1
1187#if MINTYPE==2
1188static double 
1189innerProductIndexed(const double * region1, int size, const double * region2,const int * which)
1190{
1191  int i;
1192  double value=0.0;
1193  for (i=0;i<size;i++) {
1194    int j=which[i];
1195    value += region1[j]*region2[j];
1196  }
1197  return value;
1198}
1199#endif
1200/*
1201   Row array and column array have direction
1202   Returns 0 - can do normal iteration (basis change)
1203   1 - no basis change
1204*/
1205int 
1206ClpSimplexNonlinear::pivotColumn(CoinIndexedVector * longArray,
1207                                 CoinIndexedVector * rowArray,
1208                                 CoinIndexedVector * columnArray,
1209                                 CoinIndexedVector * spare,
1210                                 int & pivotMode,
1211                                 double & solutionError,
1212                                 double * dArray)
1213{
1214  // say not optimal
1215  primalColumnPivot_->setLooksOptimal(false);
1216  double acceptablePivot=1.0e-10;
1217  int lastSequenceIn=-1;
1218  if (pivotMode&&pivotMode<10) {
1219    acceptablePivot=1.0e-6;
1220    if (factorization_->pivots())
1221      acceptablePivot=1.0e-5; // if we have iterated be more strict
1222  }
1223  double acceptableBasic=1.0e-7;
1224 
1225  int number=longArray->getNumElements();
1226  int numberTotal = numberRows_+numberColumns_;
1227  int bestSequence=-1;
1228  int bestBasicSequence=-1;
1229  double eps= 1.0e-1;
1230  eps=1.0e-6;
1231  double basicTheta = 1.0e30;
1232  double objTheta=0.0;
1233  bool finished = false;
1234  sequenceIn_=-1;
1235  int nPasses=0;
1236  int nTotalPasses=0;
1237  int nBigPasses=0;
1238  double djNorm0=0.0;
1239  double djNorm=0.0;
1240  double normFlagged=0.0;
1241  double normUnflagged=0.0;
1242  int localPivotMode=pivotMode;
1243  bool allFinished=false;
1244  bool justOne=false;
1245  int returnCode=1;
1246  double currentObj;
1247  double predictedObj;
1248  double thetaObj;
1249  objective_->stepLength(this,solution_,solution_,0.0,
1250                                               currentObj,predictedObj,thetaObj);
1251  double saveObj=currentObj;
1252#if MINTYPE ==2
1253  // try Shanno's method
1254  //would be memory leak
1255  //double * saveY=new double[numberTotal];
1256  //double * saveS=new double[numberTotal];
1257  //double * saveY2=new double[numberTotal];
1258  //double * saveS2=new double[numberTotal];
1259  double saveY[100];
1260  double saveS[100];
1261  double saveY2[100];
1262  double saveS2[100];
1263  double zz[10000];
1264#endif
1265  double * dArray2 = dArray+numberTotal;
1266  // big big loop
1267  while (!allFinished) {
1268    double * work=longArray->denseVector();
1269    int * which=longArray->getIndices();
1270    allFinished=true;
1271    // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
1272    //#define SMALLTHETA1 1.0e-25
1273    //#define SMALLTHETA2 1.0e-25
1274#define SMALLTHETA1 1.0e-10
1275#define SMALLTHETA2 1.0e-10
1276#define CONJUGATE 2
1277#if CONJUGATE == 0
1278    int conjugate = 0;
1279#elif CONJUGATE == 1
1280    int conjugate = (pivotMode<10) ? MINTYPE : 0;
1281#elif CONJUGATE == 2
1282    int conjugate = MINTYPE;
1283#else
1284    int conjugate = MINTYPE;
1285#endif
1286   
1287    //if (pivotMode==1)
1288    //localPivotMode=11;;
1289#if CLP_DEBUG > 1
1290    longArray->checkClear();
1291#endif
1292    int numberNonBasic=0;
1293    int startLocalMode=-1;
1294    while (!finished) {
1295      double simpleObjective=COIN_DBL_MAX;
1296      returnCode=1;
1297      int iSequence;
1298      objective_->reducedGradient(this,dj_,false);
1299      // get direction vector in longArray
1300      longArray->clear();
1301      // take out comment to try slightly different idea
1302      if (nPasses+nTotalPasses>3000||nBigPasses>100) {
1303        if (factorization_->pivots())
1304          returnCode=3;
1305        break;
1306      }
1307      if (!nPasses) {
1308        numberNonBasic=0;
1309        nBigPasses++;
1310      }
1311      // just do superbasic if in cleanup mode
1312      int local=localPivotMode;
1313      if (!local&&pivotMode>=10&&nBigPasses<10) {
1314        local=100;
1315      } else if (local==-1||nBigPasses>=10) {
1316        local=0;
1317        localPivotMode=0;
1318      } 
1319      if (justOne) {
1320        local=2;
1321        //local=100;
1322        justOne=false;
1323      }
1324      if (!nPasses)
1325        startLocalMode=local;
1326      directionVector(longArray,spare,rowArray,local,
1327                      normFlagged,normUnflagged,numberNonBasic);
1328      {
1329        // check null vector
1330        double * rhs = spare->denseVector();
1331#if CLP_DEBUG > 1
1332        spare->checkClear();
1333#endif
1334        int iRow;
1335        multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,rhs,0.0);
1336        matrix_->times(1.0,solution_,rhs,rowScale_,columnScale_);
1337        double largest=0.0;
1338        int iLargest=-1;
1339        for (iRow=0;iRow<numberRows_;iRow++) {
1340          double value=fabs(rhs[iRow]);
1341          rhs[iRow]=0.0;
1342          if (value>largest) {
1343            largest=value;
1344            iLargest=iRow;
1345          }
1346        }
1347#if CLP_DEBUG > 0
1348        if ((handler_->logLevel()&32)&&largest>1.0e-8) 
1349          printf("largest rhs error %g on row %d\n",largest,iLargest);
1350#endif
1351        if (solutionError<0.0) {
1352          solutionError=largest;
1353        } else if (largest>max(1.0e-8,1.0e2*solutionError)&&
1354                   factorization_->pivots()) {
1355          longArray->clear();
1356          pivotRow_ = -1;
1357          theta_=0.0;
1358          return 3;
1359        }
1360      }
1361      if (sequenceIn_>=0)
1362        lastSequenceIn=sequenceIn_;
1363#if MINTYPE!=2
1364      double djNormSave = djNorm;
1365#endif
1366      djNorm=0.0;
1367      int iIndex;
1368      for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
1369        int iSequence = which[iIndex];
1370        double alpha = work[iSequence];
1371        djNorm += alpha*alpha;
1372      }
1373      // go to conjugate gradient if necessary
1374      if (numberNonBasic&&localPivotMode>=10&&(nPasses>4||sequenceIn_<0)) {
1375        localPivotMode=0;
1376        nTotalPasses+= nPasses;
1377        nPasses=0;
1378      }
1379#if CONJUGATE == 2
1380      conjugate = (localPivotMode<10) ? MINTYPE : 0;
1381#endif
1382      //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
1383      //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
1384      if (!nPasses) {
1385        //printf("numberNon %d\n",numberNonBasic);
1386#if MINTYPE==2
1387        assert (numberNonBasic<100);
1388        memset(zz,0,numberNonBasic*numberNonBasic*sizeof(double));
1389        int put=0;
1390        for (int iVariable=0;iVariable<numberNonBasic;iVariable++) {
1391          zz[put]=1.0;
1392          put+= numberNonBasic+1;
1393        }
1394#endif
1395        djNorm0 = CoinMax(djNorm, 1.0e-20);
1396        memcpy(dArray,work,numberTotal*sizeof(double));
1397        memcpy(dArray2,work,numberTotal*sizeof(double));
1398        if (sequenceIn_>=0&&numberNonBasic==1) {
1399          // see if simple move
1400          double objTheta2 =objective_->stepLength(this,solution_,work,1.0e30,
1401                                               currentObj,predictedObj,thetaObj);
1402          rowArray->clear();
1403         
1404          // update the incoming column
1405          unpackPacked(rowArray);
1406          factorization_->updateColumnFT(spare,rowArray);
1407          theta_=0.0;
1408          double * work2=rowArray->denseVector();
1409          int number=rowArray->getNumElements();
1410          int * which2=rowArray->getIndices();
1411          int iIndex;
1412          bool easyMove=false;
1413          double way;
1414          if (dj_[sequenceIn_]>0.0)
1415            way=-1.0;
1416          else
1417            way=1.0;
1418          double largest=COIN_DBL_MAX;
1419          int kPivot=-1;
1420          for (iIndex=0;iIndex<number;iIndex++) {
1421            int iRow = which2[iIndex];
1422            double alpha = way*work2[iIndex];
1423            int iPivot = pivotVariable_[iRow];
1424            if (alpha<-1.0e-5) {
1425              double distance = upper_[iPivot]-solution_[iPivot];
1426              if (distance<-largest*alpha) {
1427                kPivot=iPivot;
1428                largest=CoinMax(0.0,-distance/alpha);
1429              }
1430              if (distance<-1.0e-12*alpha) {
1431                easyMove=true;
1432                break;
1433              }
1434            } else if (alpha>1.0e-5) {
1435              double distance = solution_[iPivot]-lower_[iPivot];
1436              if (distance<largest*alpha) {
1437                kPivot=iPivot;
1438                largest=CoinMax(0.0,distance/alpha);
1439              }
1440              if (distance<1.0e-12*alpha) {
1441                easyMove=true;
1442                break;
1443              }
1444            }
1445          }
1446          // But largest has to match up with change
1447          assert (work[sequenceIn_]);
1448          largest /= fabs(work[sequenceIn_]);
1449          if (largest<objTheta2) {
1450            easyMove=true;
1451          } else if (!easyMove) {
1452            double objDrop = currentObj-predictedObj;
1453            double th = objective_->stepLength(this,solution_,work,largest,
1454                                   currentObj,predictedObj,simpleObjective);
1455            simpleObjective = CoinMax(simpleObjective,predictedObj);
1456            double easyDrop = currentObj-simpleObjective;
1457            if (easyDrop>1.0e-8&&easyDrop>0.5*objDrop) {
1458              easyMove=true;
1459#ifdef CLP_DEBUG
1460              if (handler_->logLevel()&32) 
1461                printf("easy - obj drop %g, easy drop %g\n",objDrop,easyDrop);
1462#endif
1463              if (easyDrop>objDrop) {
1464                // debug
1465                printf("****** th %g simple %g\n",th,simpleObjective);
1466                objective_->stepLength(this,solution_,work,1.0e30,
1467                                   currentObj,predictedObj,simpleObjective);
1468                objective_->stepLength(this,solution_,work,largest,
1469                                   currentObj,predictedObj,simpleObjective);
1470              }
1471            }
1472          }
1473          rowArray->clear();
1474#ifdef CLP_DEBUG
1475          if (handler_->logLevel()&32)
1476            printf("largest %g piv %d\n",largest,kPivot);
1477#endif
1478          if (easyMove) {
1479            valueIn_=solution_[sequenceIn_];
1480            dualIn_=dj_[sequenceIn_];
1481            lowerIn_=lower_[sequenceIn_];
1482            upperIn_=upper_[sequenceIn_];
1483            if (dualIn_>0.0)
1484              directionIn_ = -1;
1485            else 
1486              directionIn_ = 1;
1487            longArray->clear();
1488            pivotRow_ = -1;
1489            theta_=0.0;
1490            return 0;
1491          }
1492        }
1493      } else {
1494#if MINTYPE==1
1495        if (conjugate) {
1496          double djNorm2 = djNorm;
1497          if (numberNonBasic&&0) {
1498            int iIndex;
1499            djNorm2 = 0.0;
1500            for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
1501              int iSequence = which[iIndex];
1502              double alpha = work[iSequence];
1503              //djNorm2 += alpha*alpha;
1504              double alpha2 = work[iSequence]-dArray2[iSequence];
1505              djNorm2 += alpha*alpha2;
1506            }
1507            //printf("a %.18g b %.18g\n",djNorm,djNorm2);
1508          }
1509          djNorm=djNorm2;
1510          double beta = djNorm2/djNormSave;
1511          // reset beta every so often
1512          //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
1513          //beta=0.0;
1514          //printf("current norm %g, old %g - beta %g\n",
1515          //     djNorm,djNormSave,beta);
1516          for (iSequence=0;iSequence<numberTotal;iSequence++) {
1517            dArray[iSequence] = work[iSequence] + beta *dArray[iSequence];
1518            dArray2[iSequence] = work[iSequence];
1519          }
1520        } else {
1521          for (iSequence=0;iSequence<numberTotal;iSequence++)
1522            dArray[iSequence] = work[iSequence];
1523        }
1524#else
1525        int number2=numberNonBasic;
1526        if (number2) {
1527          // pack down into dArray
1528          int jLast=-1;
1529          for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1530            int j=which[iSequence];
1531            assert(j>jLast);
1532            jLast=j;
1533            double value=work[j];
1534            dArray[iSequence]=-value;
1535          }
1536          // see whether to restart
1537          // check signs - gradient
1538          double g1=innerProduct(dArray,number2,dArray);
1539          double g2=innerProduct(dArray,number2,saveY2);
1540          // Get differences
1541          for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1542            saveY2[iSequence] = dArray[iSequence]-saveY2[iSequence];
1543            saveS2[iSequence] = solution_[iSequence]-saveS2[iSequence];
1544          }
1545          double g3=innerProduct(saveS2,number2,saveY2);
1546          printf("inner %g\n",g3);
1547          //assert(g3>0);
1548          double zzz[10000];
1549          int iVariable;
1550          g2=1.0e50;// temp
1551          if (fabs(g2)>=0.2*fabs(g1)) {
1552            // restart
1553            double delta = innerProduct(saveY2,number2,saveS2)/
1554              innerProduct(saveY2,number2,saveY2);
1555            delta=1.0;//temp
1556            memset(zz,0,number2*sizeof(double));
1557            int put=0;
1558            for (iVariable=0;iVariable<number2;iVariable++) {
1559              zz[put]=delta;
1560              put+= number2+1;
1561            }
1562          } else {
1563          }
1564          memcpy(zzz,zz,number2*number2*sizeof(double));
1565          double ww[100];
1566          // get sk -Hkyk
1567          for (iVariable=0;iVariable<number2;iVariable++) {
1568            double value=0.0;
1569            for (int jVariable=0;jVariable<number2;jVariable++) {
1570              value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
1571            }
1572            ww[iVariable] = saveS2[iVariable]-value;
1573          }
1574          double ys = innerProduct(saveY2,number2,saveS2);
1575          double multiplier1 = 1.0/ys;
1576          double multiplier2= innerProduct(saveY2,number2,ww)/(ys*ys);
1577#if 1
1578          // and second way
1579          // Hy
1580          double h[100];
1581          for (iVariable=0;iVariable<number2;iVariable++) {
1582            double value=0.0;
1583            for (int jVariable=0;jVariable<number2;jVariable++) {
1584              value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
1585            }
1586            h[iVariable] = value;
1587          }
1588          double hh[10000];
1589          double yhy1 = innerProduct(h,number2,saveY2)*multiplier1+1.0;
1590          yhy1 *= multiplier1;
1591          for (iVariable=0;iVariable<number2;iVariable++) {
1592            for (int jVariable=0;jVariable<number2;jVariable++) {
1593              int put = iVariable+number2*jVariable;
1594              double value = zzz[put];
1595              value += yhy1*saveS2[iVariable]*saveS2[jVariable];
1596              hh[put]=value;
1597            }
1598          }
1599          for (iVariable=0;iVariable<number2;iVariable++) {
1600            for (int jVariable=0;jVariable<number2;jVariable++) {
1601              int put = iVariable+number2*jVariable;
1602              double value = hh[put];
1603              value -= multiplier1*(saveS2[iVariable]*h[jVariable]);
1604              value -= multiplier1*(saveS2[jVariable]*h[iVariable]);
1605              hh[put]=value;
1606            }
1607          }
1608#endif
1609          // now update H
1610          for (iVariable=0;iVariable<number2;iVariable++) {
1611            for (int jVariable=0;jVariable<number2;jVariable++) {
1612              int put = iVariable+number2*jVariable;
1613              double value = zzz[put];
1614              value += multiplier1*(ww[iVariable]*saveS2[jVariable]
1615                                    +ww[jVariable]*saveS2[iVariable]);
1616              value -= multiplier2*saveS2[iVariable]*saveS2[jVariable];
1617              zzz[put]=value;
1618            }
1619          }
1620          //memcpy(zzz,hh,size*sizeof(double));
1621          // do search direction
1622          memset(dArray,0,numberTotal*sizeof(double));
1623          for (iVariable=0;iVariable<numberNonBasic;iVariable++) {
1624            double value=0.0;
1625            for (int jVariable=0;jVariable<number2;jVariable++) {
1626              int k=which[jVariable];
1627              value += work[k]*zzz[iVariable+number2*jVariable];
1628            }
1629            int i=which[iVariable];
1630            dArray[i] = value;
1631          }
1632          // Now fill out dArray
1633          {
1634            int j;
1635            // Now do basic ones
1636            int iRow;
1637            CoinIndexedVector * spare1=spare;
1638            CoinIndexedVector * spare2=rowArray;
1639#if CLP_DEBUG > 1
1640            spare1->checkClear();
1641            spare2->checkClear();
1642#endif
1643            double * array2=spare1->denseVector();
1644            int * index2 = spare1->getIndices();
1645            int number2=0;
1646            times(-1.0,dArray,array2);
1647            dArray = dArray+numberColumns_;
1648            for (iRow=0;iRow<numberRows_;iRow++) {
1649              double value = array2[iRow] + dArray[iRow];
1650              if (value) {
1651                array2[iRow]=value;
1652                index2[number2++]=iRow;
1653              } else {
1654                array2[iRow]=0.0;
1655              }
1656            }
1657            dArray -= numberColumns_;
1658            spare1->setNumElements(number2);
1659            // Ftran
1660            factorization_->updateColumn(spare2,spare1);
1661            number2=spare1->getNumElements();
1662            for (j=0;j<number2;j++) {
1663              int iSequence = index2[j];
1664              double value = array2[iSequence];
1665              array2[iSequence]=0.0;
1666              if (value) {
1667                int iPivot=pivotVariable_[iSequence];
1668                double oldValue = dArray[iPivot];
1669                dArray[iPivot] = value+oldValue;
1670              }
1671            }
1672            spare1->setNumElements(0);
1673          }
1674          //assert (innerProductIndexed(dArray,number2,work,which)>0);
1675          //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
1676          printf ("innerP1 %g innerP2 %g\n",innerProduct(dArray,numberTotal,work),
1677                  innerProductIndexed(dArray,numberNonBasic,work,which));
1678          assert (innerProduct(dArray,numberTotal,work)>0);
1679#if 1
1680          {
1681            // check null vector
1682            int iRow;
1683            double qq[10];
1684            memset(qq,0,numberRows_*sizeof(double));
1685            multiplyAdd(dArray+numberColumns_,numberRows_,-1.0,qq,0.0);
1686            matrix_->times(1.0,dArray,qq,rowScale_,columnScale_);
1687            double largest=0.0;
1688            int iLargest=-1;
1689            for (iRow=0;iRow<numberRows_;iRow++) {
1690              double value=fabs(qq[iRow]);
1691              if (value>largest) {
1692                largest=value;
1693                iLargest=iRow;
1694              }
1695            }
1696            printf("largest null error %g on row %d\n",largest,iLargest);
1697            for (iSequence=0;iSequence<numberTotal;iSequence++) {
1698              if (getStatus(iSequence)==basic)
1699                assert (fabs(dj_[iSequence])<1.0e-3);
1700            }
1701          }
1702#endif
1703          memcpy(saveY,saveY2,numberNonBasic*sizeof(double));
1704          memcpy(saveS,saveS2,numberNonBasic*sizeof(double));
1705        }
1706#endif
1707      }
1708#if MINTYPE==2
1709      for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1710        int j=which[iSequence];
1711        saveY2[iSequence]=-work[j];
1712        saveS2[iSequence]=solution_[j];
1713      }
1714#endif
1715      if (djNorm<eps*djNorm0||(nPasses>100&&djNorm<CoinMin(1.0e-1*djNorm0, 1.0e-12))) {
1716#ifdef CLP_DEBUG
1717        if (handler_->logLevel()&32) 
1718          printf("dj norm reduced from %g to %g\n",djNorm0,djNorm);
1719#endif
1720        if (pivotMode<10||!numberNonBasic) {
1721          finished=true;
1722        } else {
1723          localPivotMode=pivotMode;
1724          nTotalPasses+= nPasses;
1725          nPasses=0;
1726          continue;
1727        }
1728      }
1729      //if (nPasses==100||nPasses==50)
1730      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
1731      //         normFlagged);
1732      if (nPasses>100&&djNorm<1.0e-2*normFlagged&&!startLocalMode) {
1733#ifdef CLP_DEBUG
1734        if (handler_->logLevel()&32) 
1735          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n",djNorm0,djNorm,
1736                 normFlagged);
1737#endif
1738        unflag();
1739        localPivotMode=0;
1740        nTotalPasses+= nPasses;
1741        nPasses=0;
1742        continue;
1743      }
1744      if (djNorm>0.99*djNorm0&&nPasses>1500) {
1745        finished=true;
1746#ifdef CLP_DEBUG
1747        if (handler_->logLevel()&32) 
1748          printf("dj norm NOT reduced from %g to %g\n",djNorm0,djNorm);
1749#endif
1750        djNorm=1.2345e-20;
1751      }
1752      number=longArray->getNumElements();
1753      if (!numberNonBasic) {
1754        // looks optimal
1755        // check if any dj look plausible
1756        int nSuper=0;
1757        int nFlagged=0;
1758        int nNormal=0;
1759        for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
1760          if (flagged(iSequence)) {
1761            switch(getStatus(iSequence)) {
1762             
1763            case basic:
1764            case ClpSimplex::isFixed:
1765              break;
1766            case atUpperBound:
1767              if (dj_[iSequence]>dualTolerance_) 
1768                nFlagged++;
1769              break;
1770            case atLowerBound:
1771              if (dj_[iSequence]<-dualTolerance_) 
1772                nFlagged++;
1773              break;
1774            case isFree:
1775            case superBasic:
1776              if (fabs(dj_[iSequence])>dualTolerance_) 
1777                nFlagged++;
1778              break;
1779            }
1780            continue;
1781          }
1782          switch(getStatus(iSequence)) {
1783           
1784          case basic:
1785          case ClpSimplex::isFixed:
1786            break;
1787          case atUpperBound:
1788            if (dj_[iSequence]>dualTolerance_) 
1789              nNormal++;
1790            break;
1791          case atLowerBound:
1792            if (dj_[iSequence]<-dualTolerance_) 
1793              nNormal++;
1794            break;
1795          case isFree:
1796          case superBasic:
1797            if (fabs(dj_[iSequence])>dualTolerance_) 
1798              nSuper++;
1799            break;
1800          }
1801        }
1802#ifdef CLP_DEBUG
1803        if (handler_->logLevel()&32) 
1804          printf("%d super, %d normal, %d flagged\n",
1805                 nSuper,nNormal,nFlagged);
1806#endif
1807       
1808        int nFlagged2=1;
1809        if (lastSequenceIn<0&&!nNormal&&!nSuper) {
1810          nFlagged2=unflag();
1811          if (pivotMode>=10) {
1812            pivotMode--;
1813#ifdef CLP_DEBUG
1814            if (handler_->logLevel()&32) 
1815              printf("pivot mode now %d\n",pivotMode);
1816#endif
1817            if (pivotMode==9) 
1818              pivotMode=0;      // switch off fast attempt
1819          }
1820        } else {
1821          lastSequenceIn=-1;
1822        }
1823        if (!nFlagged2||(normFlagged+normUnflagged<1.0e-8)) {
1824          primalColumnPivot_->setLooksOptimal(true);
1825          return 0;
1826        } else {
1827          localPivotMode=-1;
1828          nTotalPasses+= nPasses;
1829          nPasses=0;
1830          finished=false;
1831          continue;
1832        }
1833      }
1834      bestSequence=-1;
1835      bestBasicSequence=-1;
1836     
1837      // temp
1838      nPasses++;
1839      if (nPasses>2000)
1840        finished=true;
1841      double theta = 1.0e30;
1842      basicTheta = 1.0e30;
1843      theta_=1.0e30;
1844      double basicTolerance = 1.0e-4*primalTolerance_;
1845      for (iSequence=0;iSequence<numberTotal;iSequence++) {
1846        //if (flagged(iSequence)
1847        //  continue;
1848        double alpha = dArray[iSequence];
1849        Status thisStatus = getStatus(iSequence);
1850        double oldValue = solution_[iSequence];
1851        if (thisStatus!=basic) {
1852          if (fabs(alpha)>=acceptablePivot) {
1853            if (alpha<0.0) {
1854              // variable going towards lower bound
1855              double bound = lower_[iSequence];
1856              oldValue -= bound;
1857              if (oldValue+theta*alpha<0.0) {
1858                bestSequence=iSequence;
1859                theta = CoinMax(0.0,oldValue/(-alpha));
1860              }
1861            } else {
1862              // variable going towards upper bound
1863              double bound = upper_[iSequence];
1864              oldValue = bound-oldValue;
1865              if (oldValue-theta*alpha<0.0) {
1866                bestSequence=iSequence;
1867                theta = CoinMax(0.0,oldValue/alpha);
1868              }
1869            }
1870          }
1871        } else {
1872          if (fabs(alpha)>=acceptableBasic) {
1873            if (alpha<0.0) {
1874              // variable going towards lower bound
1875              double bound = lower_[iSequence];
1876              oldValue -= bound;
1877              if (oldValue+basicTheta*alpha<-basicTolerance) {
1878                bestBasicSequence=iSequence;
1879                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/(-alpha));
1880              }
1881            } else {
1882              // variable going towards upper bound
1883              double bound = upper_[iSequence];
1884              oldValue = bound-oldValue;
1885              if (oldValue-basicTheta*alpha<-basicTolerance) {
1886                bestBasicSequence=iSequence;
1887                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/alpha);
1888              }
1889            }
1890          }
1891        }
1892      }
1893      theta_ = CoinMin(theta,basicTheta);
1894      // Now find minimum of function
1895      double objTheta2 =objective_->stepLength(this,solution_,dArray,CoinMin(theta,basicTheta),
1896                                               currentObj,predictedObj,thetaObj);
1897#ifdef CLP_DEBUG
1898      if (handler_->logLevel()&32) 
1899        printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
1900#endif
1901#if MINTYPE==1
1902      if (conjugate) {
1903        double offset;
1904        const double * gradient = objective_->gradient(this,
1905                                                       dArray, offset,
1906                                                       true,0);
1907        double product=0.0;
1908        for (iSequence = 0;iSequence<numberColumns_;iSequence++) {
1909          double alpha = dArray[iSequence];
1910          double value = alpha*gradient[iSequence];
1911          product += value;
1912        }
1913        //#define INCLUDESLACK
1914#ifdef INCLUDESLACK
1915        for (;iSequence<numberColumns_+numberRows_;iSequence++) {
1916          double alpha = dArray[iSequence];
1917          double value = alpha*cost_[iSequence];
1918          product += value;
1919        }
1920#endif
1921        if (product>0.0)
1922          objTheta = djNorm/product;
1923        else
1924          objTheta=COIN_DBL_MAX;
1925#ifndef NDEBUG
1926        if (product<-1.0e-8&&handler_->logLevel()>1) 
1927          printf("bad product %g\n",product);
1928#endif
1929        product = CoinMax(product,0.0);
1930      } else {
1931        objTheta =objTheta2;
1932      }
1933#else
1934      objTheta =objTheta2;
1935#endif
1936      // if very small difference then take pivot (depends on djNorm?)
1937      // distinguish between basic and non-basic
1938      bool chooseObjTheta=objTheta<theta_;
1939      if (chooseObjTheta) {
1940        if (thetaObj<currentObj-1.0e-12&&objTheta+1.0e-10>theta_)
1941          chooseObjTheta=false;
1942        //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
1943        //chooseObjTheta=false;
1944      }
1945      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
1946      if (chooseObjTheta) {
1947        theta_ = objTheta;
1948      } else {
1949        objTheta = CoinMax(objTheta,1.00000001*theta_+1.0e-12);
1950        //if (theta+1.0e-13>basicTheta) {
1951        //theta = CoinMax(theta,1.00000001*basicTheta);
1952        //theta_ = basicTheta;
1953        //}
1954      }
1955      // always out if one variable in and zero move
1956      if (theta_==basicTheta||(sequenceIn_>=0&&!theta_))
1957        finished=true; // come out
1958#ifdef CLP_DEBUG
1959      if (handler_->logLevel()&32) {
1960        printf("%d pass,",nPasses);
1961        if (sequenceIn_>=0)
1962          printf (" Sin %d,",sequenceIn_);
1963        if (basicTheta==theta_)
1964          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1965        else
1966          printf(" basicTheta %g",basicTheta);
1967        if (theta==theta_)
1968          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1969        else
1970          printf(" non-basicTheta %g",theta);
1971        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1972        printf(" djNorm %g\n",djNorm);
1973      }
1974#endif
1975      if (handler_->logLevel()>3&&objTheta!=theta_) {
1976        printf("%d pass obj %g,",nPasses,currentObj);
1977        if (sequenceIn_>=0)
1978          printf (" Sin %d,",sequenceIn_);
1979        if (basicTheta==theta_)
1980          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1981        else
1982          printf(" basicTheta %g",basicTheta);
1983        if (theta==theta_)
1984          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1985        else
1986          printf(" non-basicTheta %g",theta);
1987        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1988        printf(" djNorm %g\n",djNorm);
1989      }
1990      if (objTheta!=theta_) {
1991        //printf("hit boundary after %d passes\n",nPasses);
1992        nTotalPasses+= nPasses;
1993        nPasses=0; // start again
1994      }
1995      if (localPivotMode<10||lastSequenceIn==bestSequence) {
1996        if (theta_==theta&&theta_<basicTheta&&theta_<1.0e-5)
1997          setFlagged(bestSequence); // out of active set
1998      }
1999      // Update solution
2000      for (iSequence=0;iSequence<numberTotal;iSequence++) {
2001        //for (iIndex=0;iIndex<number;iIndex++) {
2002       
2003        //int iSequence = which[iIndex];
2004        double alpha = dArray[iSequence];
2005        if (alpha) {
2006          double value = solution_[iSequence]+theta_*alpha;
2007          solution_[iSequence]=value;
2008          switch(getStatus(iSequence)) {
2009           
2010          case basic:
2011          case isFixed:
2012          case isFree:
2013          case atUpperBound:
2014          case atLowerBound:
2015            nonLinearCost_->setOne(iSequence,value);
2016            break;
2017          case superBasic:
2018            // To get correct action
2019            setStatus(iSequence,isFixed);
2020            nonLinearCost_->setOne(iSequence,value);
2021            //assert (getStatus(iSequence)!=isFixed);
2022            break;
2023          }
2024        }
2025      }
2026      if (objTheta<SMALLTHETA2&&objTheta==theta_) {
2027        if (sequenceIn_>=0&&basicTheta<1.0e-9) {
2028          // try just one
2029          localPivotMode=0;
2030          sequenceIn_=-1;
2031          nTotalPasses+= nPasses;
2032          nPasses=0;
2033          //finished=true;
2034          //objTheta=0.0;
2035          //theta_=0.0;
2036        }
2037        //if (objTheta<1.0e-12)
2038        //finished=true;
2039        //printf("zero move\n");
2040      }
2041#ifdef CLP_DEBUG
2042      if (handler_->logLevel()&32) {
2043        objective_->stepLength(this,solution_,work,0.0,
2044                               currentObj,predictedObj,thetaObj);
2045        printf("current obj %g after update - simple was %g\n",currentObj,simpleObjective);
2046      }
2047#endif
2048      if (sequenceIn_>=0&&!finished&&objTheta>1.0e-4) {
2049        // we made some progress - back to normal
2050        if (localPivotMode<10) {
2051          localPivotMode=0;
2052          sequenceIn_=-1;
2053          nTotalPasses+= nPasses;
2054          nPasses=0;
2055        }
2056#ifdef CLP_DEBUG
2057        if (handler_->logLevel()&32) 
2058          printf("some progress?\n");
2059#endif
2060      }
2061#if CLP_DEBUG > 1
2062      longArray->checkClean();
2063#endif
2064    }
2065#ifdef CLP_DEBUG
2066    if (handler_->logLevel()&32) 
2067      printf("out of loop after %d passes\n",nPasses);
2068#endif
2069    bool ordinaryDj=false;
2070    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
2071    //printf("could try ordinary iteration %g\n",theta_);
2072    if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-15) {
2073      //printf("trying ordinary iteration\n");
2074      ordinaryDj=true;
2075    }
2076    if (!basicTheta&&!ordinaryDj) {
2077      //returnCode=2;
2078      //objTheta=-1.0; // so we fall through
2079    }
2080    assert (theta_<1.0e30); // for now
2081    // See if we need to pivot
2082    if (theta_==basicTheta||ordinaryDj) {
2083      if (!ordinaryDj) {
2084        // find an incoming column which will force pivot
2085        int iRow;
2086        pivotRow_=-1;
2087        for (iRow=0;iRow<numberRows_;iRow++) {
2088          if (pivotVariable_[iRow]==bestBasicSequence) {
2089            pivotRow_=iRow;
2090            break;
2091          }
2092        }
2093        assert (pivotRow_>=0);
2094        // Get good size for pivot
2095        double acceptablePivot=1.0e-7;
2096        if (factorization_->pivots()>10)
2097          acceptablePivot=1.0e-5; // if we have iterated be more strict
2098        else if (factorization_->pivots()>5)
2099          acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
2100        // should be dArray but seems better this way!
2101        double direction = work[bestBasicSequence]>0.0 ? -1.0 : 1.0;
2102        // create as packed
2103        rowArray->createPacked(1,&pivotRow_,&direction);
2104        factorization_->updateColumnTranspose(spare,rowArray);
2105        // put row of tableau in rowArray and columnArray
2106        matrix_->transposeTimes(this,-1.0,
2107                                rowArray,spare,columnArray);
2108        // choose one futhest away from bound which has reasonable pivot
2109        // If increasing we want negative alpha
2110       
2111        double * work2;
2112        int iSection;
2113       
2114        sequenceIn_=-1;
2115        double bestValue=-1.0;
2116        double bestDirection=0.0;
2117        // First pass we take correct direction and large pivots
2118        // then correct direction
2119        // then any
2120        double check[]={1.0e-8,-1.0e-12,-1.0e30};
2121        double mult[]={100.0,1.0,1.0};
2122        for (int iPass=0;iPass<3;iPass++) {
2123          //if (!bestValue&&iPass==2)
2124          //bestValue=-1.0;
2125          double acceptable=acceptablePivot*mult[iPass];
2126          double checkValue=check[iPass];
2127          for (iSection=0;iSection<2;iSection++) {
2128           
2129            int addSequence;
2130           
2131            if (!iSection) {
2132              work2 = rowArray->denseVector();
2133              number = rowArray->getNumElements();
2134              which = rowArray->getIndices();
2135              addSequence = numberColumns_;
2136            } else {
2137              work2 = columnArray->denseVector();
2138              number = columnArray->getNumElements();
2139              which = columnArray->getIndices();
2140              addSequence = 0;
2141            }
2142            int i;
2143           
2144            for (i=0;i<number;i++) {
2145              int iSequence = which[i]+addSequence;
2146              if (flagged(iSequence))
2147                continue;
2148              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
2149              //                  upper_[iSequence]-solution_[iSequence]);
2150              double alpha=work2[i];
2151              // should be dArray but seems better this way!
2152              double change=work[iSequence];
2153              Status thisStatus = getStatus(iSequence);
2154              double direction=0;;
2155              switch(thisStatus) {
2156               
2157              case basic:
2158              case ClpSimplex::isFixed:
2159                break;
2160              case isFree:
2161              case superBasic:
2162                if (alpha<-acceptable&&change>checkValue)
2163                  direction=1.0;
2164                else if (alpha>acceptable&&change<-checkValue)
2165                  direction=-1.0;
2166                break;
2167              case atUpperBound:
2168                if (alpha>acceptable&&change<-checkValue)
2169                  direction=-1.0;
2170                else if (iPass==2&&alpha<-acceptable&&change<-checkValue)
2171                  direction=1.0;
2172                break;
2173              case atLowerBound:
2174                if (alpha<-acceptable&&change>checkValue)
2175                  direction=1.0;
2176                else if (iPass==2&&alpha>acceptable&&change>checkValue)
2177                  direction=-1.0;
2178                break;
2179              }
2180              if (direction) {
2181                if (sequenceIn_!=lastSequenceIn||localPivotMode<10) { 
2182                  if (CoinMin(solution_[iSequence]-lower_[iSequence],
2183                          upper_[iSequence]-solution_[iSequence])>bestValue) {
2184                    bestValue=CoinMin(solution_[iSequence]-lower_[iSequence],
2185                                  upper_[iSequence]-solution_[iSequence]);
2186                    sequenceIn_=iSequence;
2187                    bestDirection=direction;
2188                  }
2189                } else {
2190                  // choose
2191                  bestValue=COIN_DBL_MAX;
2192                  sequenceIn_=iSequence;
2193                  bestDirection=direction;
2194                }
2195              }
2196            }
2197          }
2198          if (sequenceIn_>=0&&bestValue>0.0)
2199            break;
2200        }
2201        if (sequenceIn_>=0) {
2202          valueIn_=solution_[sequenceIn_];
2203          dualIn_=dj_[sequenceIn_];
2204          if (bestDirection<0.0) {
2205            // we want positive dj
2206            if (dualIn_<=0.0) {
2207              //printf("bad dj - xx %g\n",dualIn_);
2208              dualIn_=1.0;
2209            }
2210          } else {
2211            // we want negative dj
2212            if (dualIn_>=0.0) {
2213              //printf("bad dj - xx %g\n",dualIn_);
2214              dualIn_=-1.0;
2215            }
2216          }
2217          lowerIn_=lower_[sequenceIn_];
2218          upperIn_=upper_[sequenceIn_];
2219          if (dualIn_>0.0)
2220            directionIn_ = -1;
2221          else 
2222            directionIn_ = 1;
2223        } else {
2224          //ordinaryDj=true;
2225#ifdef CLP_DEBUG
2226          if (handler_->logLevel()&32) {
2227            printf("no easy pivot - norm %g mode %d\n",djNorm,localPivotMode);
2228            if (rowArray->getNumElements()+columnArray->getNumElements()<12) {
2229              for (iSection=0;iSection<2;iSection++) {
2230               
2231                int addSequence;
2232               
2233                if (!iSection) {
2234                  work2 = rowArray->denseVector();
2235                  number = rowArray->getNumElements();
2236                  which = rowArray->getIndices();
2237                  addSequence = numberColumns_;
2238                } else {
2239                  work2 = columnArray->denseVector();
2240                  number = columnArray->getNumElements();
2241                  which = columnArray->getIndices();
2242                  addSequence = 0;
2243                }
2244                int i;
2245                char section[]={'R','C'};
2246                for (i=0;i<number;i++) {
2247                  int iSequence = which[i]+addSequence;
2248                  if (flagged(iSequence)) {
2249                    printf("%c%d flagged\n",section[iSection],which[i]);
2250                    continue;
2251                  } else {
2252                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
2253                           section[iSection],which[i],status_[iSequence],
2254                           lower_[iSequence],solution_[iSequence],upper_[iSequence],
2255                           work2[i],work[iSequence]);
2256                  }
2257                }
2258              }
2259            }
2260          }
2261#endif
2262          assert (sequenceIn_<0);
2263          justOne=true;
2264          allFinished=false; // Round again
2265          finished=false;
2266          nTotalPasses+= nPasses;
2267          nPasses=0;
2268          if (djNorm<0.9*djNorm0&&djNorm<1.0e-3&&!localPivotMode) {
2269#ifdef CLP_DEBUG
2270            if (handler_->logLevel()&32) 
2271              printf("no pivot - mode %d norms %g %g - unflagging\n",
2272                     localPivotMode,djNorm0,djNorm);
2273#endif
2274            unflag(); //unflagging
2275            returnCode=1;
2276          } else {
2277            returnCode=2; // do single incoming
2278            //returnCode=1;
2279          }
2280        }
2281      }
2282      rowArray->clear();
2283      columnArray->clear();
2284      longArray->clear();
2285      if (ordinaryDj) {
2286        valueIn_=solution_[sequenceIn_];
2287        dualIn_=dj_[sequenceIn_];
2288        lowerIn_=lower_[sequenceIn_];
2289        upperIn_=upper_[sequenceIn_];
2290        if (dualIn_>0.0)
2291          directionIn_ = -1;
2292        else 
2293          directionIn_ = 1;
2294      }
2295      if (returnCode==1)
2296        returnCode=0;
2297    } else {
2298      // round again
2299      longArray->clear();
2300      if (djNorm<1.0e-3&&!localPivotMode) {
2301        if (djNorm==1.2345e-20&&djNorm0>1.0e-4) {
2302#ifdef CLP_DEBUG
2303          if (handler_->logLevel()&32) 
2304            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n",djNorm0,nPasses,
2305                   localPivotMode,returnCode);
2306#endif
2307          //if (!localPivotMode)
2308          //returnCode=2; // force singleton
2309        } else {
2310#ifdef CLP_DEBUG
2311          if (handler_->logLevel()&32) 
2312            printf("unflagging as djNorm %g %g, %d passes\n",djNorm,djNorm0,nPasses);
2313#endif
2314          if (pivotMode>=10) {
2315            pivotMode--;
2316#ifdef CLP_DEBUG
2317            if (handler_->logLevel()&32) 
2318              printf("pivot mode now %d\n",pivotMode);
2319#endif
2320            if (pivotMode==9) 
2321              pivotMode=0;      // switch off fast attempt
2322          }
2323          bool unflagged = unflag()!=0;
2324          if (!unflagged&&djNorm<1.0e-10) {
2325            // ? declare victory
2326            sequenceIn_=-1;
2327            returnCode=0;
2328          }
2329        }
2330      }
2331    }
2332  }
2333  if (djNorm0<1.0e-12*normFlagged) {
2334#ifdef CLP_DEBUG
2335    if (handler_->logLevel()&32) 
2336      printf("unflagging as djNorm %g %g and flagged norm %g\n",djNorm,djNorm0,normFlagged);
2337#endif
2338    unflag();
2339  }
2340  if (saveObj-currentObj<1.0e-5&&nTotalPasses>2000) {
2341    normUnflagged=0.0;
2342    double dualTolerance3 = CoinMin(1.0e-2,1.0e3*dualTolerance_);
2343    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
2344      switch(getStatus(iSequence)) {
2345       
2346      case basic:
2347      case ClpSimplex::isFixed:
2348        break;
2349      case atUpperBound:
2350        if (dj_[iSequence]>dualTolerance3) 
2351          normFlagged += dj_[iSequence]*dj_[iSequence];
2352        break;
2353      case atLowerBound:
2354        if (dj_[iSequence]<-dualTolerance3) 
2355          normFlagged += dj_[iSequence]*dj_[iSequence];
2356        break;
2357      case isFree:
2358      case superBasic:
2359        if (fabs(dj_[iSequence])>dualTolerance3) 
2360          normFlagged += dj_[iSequence]*dj_[iSequence];
2361        break;
2362      }
2363    }
2364    if (handler_->logLevel()>2)
2365      printf("possible optimal  %d %d %g %g\n",
2366             nBigPasses,nTotalPasses,saveObj-currentObj,normFlagged);
2367    if (normFlagged<1.0e-5) {
2368      unflag();
2369      primalColumnPivot_->setLooksOptimal(true);
2370      returnCode=0;
2371    }
2372  }
2373  return returnCode;
2374}
2375/* Do last half of an iteration.
2376   Return codes
2377   Reasons to come out normal mode
2378   -1 normal
2379   -2 factorize now - good iteration
2380   -3 slight inaccuracy - refactorize - iteration done
2381   -4 inaccuracy - refactorize - no iteration
2382   -5 something flagged - go round again
2383   +2 looks unbounded
2384   +3 max iterations (iteration done)
2385   
2386*/
2387int 
2388ClpSimplexNonlinear::pivotNonlinearResult()
2389{
2390
2391  int returnCode=-1;
2392
2393  rowArray_[1]->clear();
2394 
2395  // we found a pivot column
2396  // update the incoming column
2397  unpackPacked(rowArray_[1]);
2398  factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2399  theta_=0.0;
2400  double * work=rowArray_[1]->denseVector();
2401  int number=rowArray_[1]->getNumElements();
2402  int * which=rowArray_[1]->getIndices();
2403  bool keepValue=false;
2404  double saveValue=0.0;
2405  if (pivotRow_>=0) {
2406    sequenceOut_=pivotVariable_[pivotRow_];;
2407    valueOut_ = solution(sequenceOut_);
2408    keepValue=true;
2409    saveValue=valueOut_;
2410    lowerOut_=lower_[sequenceOut_];
2411    upperOut_=upper_[sequenceOut_];
2412    int iIndex;
2413    for (iIndex=0;iIndex<number;iIndex++) {
2414     
2415      int iRow = which[iIndex];
2416      if (iRow==pivotRow_) {
2417        alpha_ = work[iIndex];
2418        break;
2419      }
2420    }
2421  } else {
2422    int iIndex;
2423    double smallest=COIN_DBL_MAX;
2424    for (iIndex=0;iIndex<number;iIndex++) {
2425      int iRow = which[iIndex];
2426      double alpha = work[iIndex];
2427      if (fabs(alpha)>1.0e-6) {
2428        int iPivot = pivotVariable_[iRow];
2429        double distance = CoinMin(upper_[iPivot]-solution_[iPivot],
2430                              solution_[iPivot]-lower_[iPivot]);
2431        if (distance<smallest) {
2432          pivotRow_=iRow;
2433          alpha_=alpha;
2434          smallest=distance;
2435        }
2436      }
2437    }
2438    if (smallest>primalTolerance_) {
2439      smallest=COIN_DBL_MAX;
2440      for (iIndex=0;iIndex<number;iIndex++) {
2441        int iRow = which[iIndex];
2442        double alpha = work[iIndex];
2443        if (fabs(alpha)>1.0e-6) {
2444          double distance = randomNumberGenerator_.randomDouble();
2445          if (distance<smallest) {
2446            pivotRow_=iRow;
2447            alpha_=alpha;
2448            smallest=distance;
2449          }
2450        }
2451      }
2452    }
2453    assert (pivotRow_>=0);
2454    sequenceOut_=pivotVariable_[pivotRow_];;
2455    valueOut_ = solution(sequenceOut_);
2456    lowerOut_=lower_[sequenceOut_];
2457    upperOut_=upper_[sequenceOut_];
2458  }
2459  double newValue = valueOut_ - theta_*alpha_;
2460  bool isSuperBasic=false;
2461  if (valueOut_>=upperOut_-primalTolerance_) {
2462    directionOut_=-1;      // to upper bound
2463    upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2464    upperOut_ = newValue;
2465  } else if (valueOut_<=lowerOut_+primalTolerance_) {
2466    directionOut_=1;      // to lower bound
2467    lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2468  } else {
2469    lowerOut_=valueOut_;
2470    upperOut_=valueOut_;
2471    isSuperBasic=true;
2472    //printf("XX superbasic out\n");
2473  }
2474  dualOut_ = dj_[sequenceOut_];
2475  double checkValue=1.0e-2;
2476  if (largestDualError_>1.0e-5)
2477    checkValue=1.0e-1;
2478  // if stable replace in basis
2479 
2480  int updateStatus = factorization_->replaceColumn(this,
2481                                                   rowArray_[2],
2482                                                   rowArray_[1],
2483                                                   pivotRow_,
2484                                                   alpha_);
2485 
2486  // if no pivots, bad update but reasonable alpha - take and invert
2487  if (updateStatus==2&&
2488      lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2489    updateStatus=4;
2490  if (updateStatus==1||updateStatus==4) {
2491    // slight error
2492    if (factorization_->pivots()>5||updateStatus==4) {
2493      returnCode=-3;
2494    }
2495  } else if (updateStatus==2) {
2496    // major error
2497    // better to have small tolerance even if slower
2498    factorization_->zeroTolerance(1.0e-15);
2499    int maxFactor = factorization_->maximumPivots();
2500    if (maxFactor>10) {
2501      if (forceFactorization_<0)
2502        forceFactorization_= maxFactor;
2503      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2504    } 
2505    // later we may need to unwind more e.g. fake bounds
2506    if(lastGoodIteration_ != numberIterations_) {
2507      clearAll();
2508      pivotRow_=-1;
2509      returnCode=-4;
2510    } else {
2511      // need to reject something
2512      char x = isColumn(sequenceIn_) ? 'C' :'R';
2513      handler_->message(CLP_SIMPLEX_FLAG,messages_)
2514        <<x<<sequenceWithin(sequenceIn_)
2515        <<CoinMessageEol;
2516      setFlagged(sequenceIn_);
2517      progress_->clearBadTimes();
2518      lastBadIteration_ = numberIterations_; // say be more cautious
2519      clearAll();
2520      pivotRow_=-1;
2521      sequenceOut_=-1;
2522      returnCode = -5;
2523     
2524    }
2525    return returnCode;
2526  } else if (updateStatus==3) {
2527    // out of memory
2528    // increase space if not many iterations
2529    if (factorization_->pivots()<
2530        0.5*factorization_->maximumPivots()&&
2531        factorization_->pivots()<200)
2532      factorization_->areaFactor(
2533                                 factorization_->areaFactor() * 1.1);
2534    returnCode =-2; // factorize now
2535  } else if (updateStatus==5) {
2536    problemStatus_=-2; // factorize now
2537  }
2538 
2539  // update primal solution
2540 
2541  double objectiveChange=0.0;
2542  // after this rowArray_[1] is not empty - used to update djs
2543  // If pivot row >= numberRows then may be gub
2544  updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1);
2545 
2546  double oldValue = valueIn_;
2547  if (directionIn_==-1) {
2548    // as if from upper bound
2549    if (sequenceIn_!=sequenceOut_) {
2550      // variable becoming basic
2551      valueIn_ -= fabs(theta_);
2552    } else {
2553      valueIn_=lowerIn_;
2554    }
2555  } else {
2556    // as if from lower bound
2557    if (sequenceIn_!=sequenceOut_) {
2558      // variable becoming basic
2559      valueIn_ += fabs(theta_);
2560    } else {
2561      valueIn_=upperIn_;
2562    }
2563  }
2564  objectiveChange += dualIn_*(valueIn_-oldValue);
2565  // outgoing
2566  if (sequenceIn_!=sequenceOut_) {
2567    if (directionOut_>0) {
2568      valueOut_ = lowerOut_;
2569    } else {
2570      valueOut_ = upperOut_;
2571    }
2572    if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2573      valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2574    else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2575      valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2576    // may not be exactly at bound and bounds may have changed
2577    // Make sure outgoing looks feasible
2578    if (!isSuperBasic)
2579      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2580    solution_[sequenceOut_]=valueOut_;
2581  }
2582  // change cost and bounds on incoming if primal
2583  nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2584  int whatNext=housekeeping(objectiveChange);
2585  if (keepValue)
2586    solution_[sequenceOut_]=saveValue;
2587  if (isSuperBasic)
2588    setStatus(sequenceOut_,superBasic);
2589  //#define CLP_DEBUG
2590#if CLP_DEBUG > 1
2591  {
2592    int ninf= matrix_->checkFeasible(this);
2593    if (ninf)
2594      printf("infeas %d\n",ninf);
2595  }
2596#endif
2597  if (whatNext==1) {
2598    returnCode =-2; // refactorize
2599  } else if (whatNext==2) {
2600    // maximum iterations or equivalent
2601    returnCode=3;
2602  } else if(numberIterations_ == lastGoodIteration_
2603            + 2 * factorization_->maximumPivots()) {
2604    // done a lot of flips - be safe
2605    returnCode =-2; // refactorize
2606  }
2607  // Check event
2608  {
2609    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2610    if (status>=0) {
2611      problemStatus_=5;
2612      secondaryStatus_=ClpEventHandler::endOfIteration;
2613      returnCode=4;
2614    }
2615  }
2616  return returnCode;
2617}
2618// A sequential LP method
2619int 
2620ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
2621{
2622  // Are we minimizing or maximizing
2623  double whichWay=optimizationDirection();
2624  if (whichWay<0.0)
2625    whichWay=-1.0;
2626  else if (whichWay>0.0)
2627    whichWay=1.0;
2628
2629
2630  int numberColumns = this->numberColumns();
2631  int numberRows = this->numberRows();
2632  double * columnLower = this->columnLower();
2633  double * columnUpper = this->columnUpper();
2634  double * solution = this->primalColumnSolution();
2635 
2636  if (objective_->type()<2) {
2637    // no nonlinear part
2638    return ClpSimplex::primal(0);
2639  }
2640  // Get list of non linear columns
2641  char * markNonlinear = new char[numberColumns];
2642  memset(markNonlinear,0,numberColumns);
2643  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
2644 
2645  if (!numberNonLinearColumns) {
2646    delete [] markNonlinear;
2647    // no nonlinear part
2648    return ClpSimplex::primal(0);
2649  }
2650  int iColumn;
2651  int * listNonLinearColumn = new int[numberNonLinearColumns];
2652  numberNonLinearColumns=0;
2653  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2654    if(markNonlinear[iColumn])
2655      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
2656  }
2657  // Replace objective
2658  ClpObjective * trueObjective = objective_;
2659  objective_=new ClpLinearObjective(NULL,numberColumns);
2660  double * objective = this->objective();
2661 
2662  // get feasible
2663  if (this->status()<0||numberPrimalInfeasibilities())
2664    ClpSimplex::primal(1);
2665  // still infeasible
2666  if (numberPrimalInfeasibilities()) {
2667    delete [] listNonLinearColumn;
2668    return 0;
2669  }
2670  algorithm_=1;
2671  int jNon;
2672  int * last[3];
2673 
2674  double * trust = new double[numberNonLinearColumns];
2675  double * trueLower = new double[numberNonLinearColumns];
2676  double * trueUpper = new double[numberNonLinearColumns];
2677  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2678    iColumn=listNonLinearColumn[jNon];
2679    trust[jNon]=0.5;
2680    trueLower[jNon]=columnLower[iColumn];
2681    trueUpper[jNon]=columnUpper[iColumn];
2682    if (solution[iColumn]<trueLower[jNon])
2683      solution[iColumn]=trueLower[jNon];
2684    else if (solution[iColumn]>trueUpper[jNon])
2685      solution[iColumn]=trueUpper[jNon];
2686  }
2687  int saveLogLevel=logLevel();
2688  int iPass;
2689  double lastObjective=1.0e31;
2690  double * saveSolution = new double [numberColumns];
2691  double * saveRowSolution = new double [numberRows];
2692  memset(saveRowSolution,0,numberRows*sizeof(double));
2693  double * savePi = new double [numberRows];
2694  double * safeSolution = new double [numberColumns];
2695  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
2696#define MULTIPLE 0
2697#if MULTIPLE>2
2698  // Duplication but doesn't really matter
2699  double * saveSolutionM[MULTIPLE};
2700  for (jNon=0;jNon<MULTIPLE;jNon++) {
2701    saveSolutionM[jNon]=new double[numberColumns];
2702    memcpy(saveSolutionM,solution,numberColumns*sizeof(double));
2703  }
2704#endif
2705  double targetDrop=1.0e31;
2706  double objectiveOffset;
2707  getDblParam(ClpObjOffset,objectiveOffset);
2708  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
2709  for (iPass=0;iPass<3;iPass++) {
2710    last[iPass]=new int[numberNonLinearColumns];
2711    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
2712      last[iPass][jNon]=0;
2713  }
2714  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2715  int goodMove=-2;
2716  char * statusCheck = new char[numberColumns];
2717  double * changeRegion = new double [numberColumns];
2718  double offset=0.0;
2719  double objValue=0.0;
2720  int exitPass = 2*numberPasses+10;
2721  for (iPass=0;iPass<numberPasses;iPass++) {
2722    exitPass--;
2723    // redo objective
2724    offset=0.0;
2725    objValue=-objectiveOffset;
2726    // make sure x updated
2727    trueObjective->newXValues();
2728    double theta=-1.0;
2729    double maxTheta=COIN_DBL_MAX;
2730    //maxTheta=1.0;
2731    if (iPass) {
2732      int jNon=0;
2733      for (iColumn=0;iColumn<numberColumns;iColumn++) { 
2734        changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
2735        double alpha = changeRegion[iColumn];
2736        double oldValue = saveSolution[iColumn];
2737        if (markNonlinear[iColumn]==0) {
2738          // linear
2739          if (alpha<-1.0e-15) {
2740            // variable going towards lower bound
2741            double bound = columnLower[iColumn];
2742            oldValue -= bound;
2743            if (oldValue+maxTheta*alpha<0.0) {
2744              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2745            }
2746          } else if (alpha>1.0e-15) {
2747            // variable going towards upper bound
2748            double bound = columnUpper[iColumn];
2749            oldValue = bound-oldValue;
2750            if (oldValue-maxTheta*alpha<0.0) {
2751              maxTheta = CoinMax(0.0,oldValue/alpha);
2752            }
2753          }
2754        } else {
2755          // nonlinear
2756          if (alpha<-1.0e-15) {
2757            // variable going towards lower bound
2758            double bound = trueLower[jNon];
2759            oldValue -= bound;
2760            if (oldValue+maxTheta*alpha<0.0) {
2761              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2762            }
2763          } else if (alpha>1.0e-15) {
2764            // variable going towards upper bound
2765            double bound = trueUpper[jNon];
2766            oldValue = bound-oldValue;
2767            if (oldValue-maxTheta*alpha<0.0) {
2768              maxTheta = CoinMax(0.0,oldValue/alpha);
2769            }
2770          }
2771          jNon++;
2772        }
2773      }
2774      // make sure both accurate
2775      memset(rowActivity_,0,numberRows_*sizeof(double));
2776      times(1.0,solution,rowActivity_);
2777      memset(saveRowSolution,0,numberRows_*sizeof(double));
2778      times(1.0,saveSolution,saveRowSolution);
2779      for (int iRow=0;iRow<numberRows;iRow++) { 
2780        double alpha =rowActivity_[iRow]-saveRowSolution[iRow];
2781        double oldValue = saveRowSolution[iRow];
2782        if (alpha<-1.0e-15) {
2783          // variable going towards lower bound
2784          double bound = rowLower_[iRow];
2785          oldValue -= bound;
2786          if (oldValue+maxTheta*alpha<0.0) {
2787            maxTheta = CoinMax(0.0,oldValue/(-alpha));
2788          }
2789        } else if (alpha>1.0e-15) {
2790          // variable going towards upper bound
2791          double bound = rowUpper_[iRow];
2792          oldValue = bound-oldValue;
2793          if (oldValue-maxTheta*alpha<0.0) {
2794            maxTheta = CoinMax(0.0,oldValue/alpha);
2795          }
2796        }
2797      }
2798    } else {
2799      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2800        changeRegion[iColumn]=0.0;
2801        saveSolution[iColumn]=solution[iColumn];
2802      }
2803      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
2804    }
2805    // get current value anyway
2806    double predictedObj,thetaObj;
2807    double maxTheta2 = 2.0; // to work out a b c
2808    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
2809                                             objValue,predictedObj,thetaObj);
2810    int lastMoveStatus = goodMove;
2811    if (goodMove>=0) {
2812      theta = CoinMin(theta2,maxTheta);
2813#ifdef CLP_DEBUG
2814      if (handler_->logLevel()&32) 
2815        printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
2816               theta,objValue,thetaObj,predictedObj);
2817#endif
2818      if (theta>0.0&&theta<=1.0) {
2819        // update solution
2820        double lambda = 1.0-theta;
2821        for (iColumn=0;iColumn<numberColumns;iColumn++) 
2822          solution[iColumn] = lambda * saveSolution[iColumn] 
2823            + theta * solution[iColumn];
2824        memset(rowActivity_,0,numberRows_*sizeof(double));
2825        times(1.0,solution,rowActivity_);
2826        if (lambda>0.999) {
2827          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
2828          memcpy(status_,saveStatus,numberRows+numberColumns);
2829        }
2830        // Do local minimization
2831#define LOCAL
2832#ifdef LOCAL
2833        bool absolutelyOptimal=false;
2834        int saveScaling = scalingFlag();
2835        scaling(0);
2836        int savePerturbation=perturbation_;
2837        perturbation_=100;
2838        if (saveLogLevel==1)
2839          setLogLevel(0);
2840        int status=startup(1);
2841        if (!status) {
2842          int numberTotal = numberRows_+numberColumns_;
2843          // resize arrays
2844          for (int i=0;i<4;i++) {
2845            rowArray_[i]->reserve(CoinMax(numberRows_+numberColumns_,rowArray_[i]->capacity()));
2846          }
2847          CoinIndexedVector * longArray = rowArray_[3];
2848          CoinIndexedVector * rowArray = rowArray_[0];
2849          //CoinIndexedVector * columnArray = columnArray_[0];
2850          CoinIndexedVector * spare = rowArray_[1];
2851          double * work=longArray->denseVector();
2852          int * which=longArray->getIndices();
2853          int nPass=100;
2854          //bool conjugate=false;
2855          // Put back true bounds
2856          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2857            int iColumn=listNonLinearColumn[jNon];
2858            double value;
2859            value=trueLower[jNon];
2860            trueLower[jNon]=lower_[iColumn];
2861            lower_[iColumn]=value;
2862            value=trueUpper[jNon];
2863            trueUpper[jNon]=upper_[iColumn];
2864            upper_[iColumn]=value;
2865            switch(getStatus(iColumn)) {
2866             
2867            case basic:
2868            case isFree:
2869            case superBasic:
2870              break;
2871            case isFixed:
2872            case atUpperBound:
2873            case atLowerBound:
2874              if (solution_[iColumn]>lower_[iColumn]+primalTolerance_) {
2875                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2876                  setStatus(iColumn,superBasic);
2877                } else {
2878                  setStatus(iColumn,atUpperBound);
2879                }
2880              } else {
2881                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2882                  setStatus(iColumn,atLowerBound);
2883                } else {
2884                  setStatus(iColumn,isFixed);
2885                }
2886              }
2887              break;
2888            }
2889          }
2890          for (int jPass=0;jPass<nPass;jPass++) {
2891            trueObjective->reducedGradient(this,dj_,true);
2892            // get direction vector in longArray
2893            longArray->clear();
2894            double normFlagged=0.0;
2895            double normUnflagged=0.0;
2896            int numberNonBasic=0;
2897            directionVector(longArray,spare,rowArray,0,
2898                            normFlagged,normUnflagged,numberNonBasic);
2899            if (normFlagged+normUnflagged<1.0e-8) {
2900              absolutelyOptimal=true;
2901              break;  //looks optimal
2902            }
2903            double djNorm=0.0;
2904            int iIndex;
2905            for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
2906              int iSequence = which[iIndex];
2907              double alpha = work[iSequence];
2908              djNorm += alpha*alpha;
2909            }
2910            //if (!jPass)
2911            //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
2912            //int number=longArray->getNumElements();
2913            if (!numberNonBasic) {
2914              // looks optimal
2915              absolutelyOptimal=true;
2916              break;
2917            }
2918            int bestSequence=-1;
2919            double theta = 1.0e30;
2920            int iSequence;
2921            for (iSequence=0;iSequence<numberTotal;iSequence++) {
2922              double alpha = work[iSequence];
2923              double oldValue = solution_[iSequence];
2924              if (alpha<-1.0e-15) {
2925                // variable going towards lower bound
2926                double bound = lower_[iSequence];
2927                oldValue -= bound;
2928                if (oldValue+theta*alpha<0.0) {
2929                  bestSequence=iSequence;
2930                  theta = CoinMax(0.0,oldValue/(-alpha));
2931                }
2932              } else if (alpha>1.0e-15) {
2933                // variable going towards upper bound
2934                double bound = upper_[iSequence];
2935                oldValue = bound-oldValue;
2936                if (oldValue-theta*alpha<0.0) {
2937                  bestSequence=iSequence;
2938                  theta = CoinMax(0.0,oldValue/alpha);
2939                }
2940              }
2941            }
2942            // Now find minimum of function
2943            double currentObj;
2944            double predictedObj;
2945            double thetaObj;
2946            // need to use true objective
2947            double * saveCost = cost_;
2948            cost_=NULL;
2949            double objTheta =trueObjective->stepLength(this,solution_,work,theta,
2950                                                       currentObj,predictedObj,thetaObj);
2951            cost_=saveCost;
2952#ifdef CLP_DEBUG
2953            if (handler_->logLevel()&32) 
2954              printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2955#endif
2956            //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2957            //printf("objTheta %g theta %g\n",objTheta,theta);
2958            if (theta>objTheta) { 
2959              theta=objTheta;
2960              thetaObj=predictedObj;
2961            }
2962            // update one used outside
2963            objValue=currentObj;
2964            if (theta>1.0e-9&&
2965                (currentObj-thetaObj<-CoinMax(1.0e-8,1.0e-15*fabs(currentObj))||jPass<5)) {
2966              // Update solution
2967              for (iSequence=0;iSequence<numberTotal;iSequence++) {
2968                double alpha = work[iSequence];
2969                if (alpha) {
2970                  double value = solution_[iSequence]+theta*alpha;
2971                  solution_[iSequence]=value;
2972                  switch(getStatus(iSequence)) {
2973                   
2974                  case basic:
2975                  case isFixed:
2976                  case isFree:
2977                    break;
2978                  case atUpperBound:
2979                  case atLowerBound:
2980                  case superBasic:
2981                    if (value>lower_[iSequence]+primalTolerance_) {
2982                      if(value<upper_[iSequence]-primalTolerance_) {
2983                        setStatus(iSequence,superBasic);
2984                      } else {
2985                        setStatus(iSequence,atUpperBound);
2986                      }
2987                    } else {
2988                      if(value<upper_[iSequence]-primalTolerance_) {
2989                        setStatus(iSequence,atLowerBound);
2990                      } else {
2991                        setStatus(iSequence,isFixed);
2992                      }
2993                    }
2994                    break;
2995                  }
2996                }
2997              }
2998            } else {
2999              break;
3000            }
3001          }
3002          // Put back fake bounds
3003          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3004            int iColumn=listNonLinearColumn[jNon];
3005            double value;
3006            value=trueLower[jNon];
3007            trueLower[jNon]=lower_[iColumn];
3008            lower_[iColumn]=value;
3009            value=trueUpper[jNon];
3010            trueUpper[jNon]=upper_[iColumn];
3011            upper_[iColumn]=value;
3012          }
3013        }
3014        problemStatus_=0;
3015        finish();
3016        scaling(saveScaling);
3017        perturbation_=savePerturbation;
3018        setLogLevel(saveLogLevel);
3019#endif
3020        // redo rowActivity
3021        memset(rowActivity_,0,numberRows_*sizeof(double));
3022        times(1.0,solution,rowActivity_);
3023        if (theta>0.99999&&theta2<1.9&&!absolutelyOptimal) {
3024          // If big changes then tighten
3025          /*  thetaObj is objvalue + a*2*2 +b*2
3026              predictedObj is objvalue + a*theta2*theta2 +b*theta2
3027          */
3028          double rhs1 = thetaObj-objValue;
3029          double rhs2 = predictedObj-objValue;
3030          double subtractB = theta2*0.5;
3031          double a = (rhs2-subtractB*rhs1)/(theta2*theta2-4.0*subtractB);
3032          double b = 0.5*(rhs1 - 4.0*a);
3033          if (fabs(a+b)>1.0e-2) {
3034            // tighten all
3035            goodMove=-1;
3036          }
3037        }
3038      }
3039    }
3040    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3041           numberColumns*sizeof(double));
3042    //printf("offset comp %g orig %g\n",offset,objectiveOffset);
3043    // could do faster
3044    trueObjective->stepLength(this,solution,changeRegion,0.0,
3045                                             objValue,predictedObj,thetaObj);
3046    printf("offset comp %g orig %g - obj from stepLength %g\n",offset,objectiveOffset,objValue);
3047    setDblParam(ClpObjOffset,objectiveOffset+offset);
3048    int * temp=last[2];
3049    last[2]=last[1];
3050    last[1]=last[0];
3051    last[0]=temp;
3052    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3053      iColumn=listNonLinearColumn[jNon];
3054      double change = solution[iColumn]-saveSolution[iColumn];
3055      if (change<-1.0e-5) {
3056        if (fabs(change+trust[jNon])<1.0e-5) 
3057          temp[jNon]=-1;
3058        else
3059          temp[jNon]=-2;
3060      } else if(change>1.0e-5) {
3061        if (fabs(change-trust[jNon])<1.0e-5) 
3062          temp[jNon]=1;
3063        else
3064          temp[jNon]=2;
3065      } else {
3066        temp[jNon]=0;
3067      }
3068    } 
3069    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3070    double maxDelta=0.0;
3071    if (goodMove>=0) {
3072      if (objValue-lastObjective<=1.0e-15*fabs(lastObjective)) 
3073        goodMove=1;
3074      else
3075        goodMove=0;
3076    } else {
3077      maxDelta=1.0e10;
3078    }
3079    double maxGap=0.0;
3080    int numberSmaller=0;
3081    int numberSmaller2=0;
3082    int numberLarger=0;
3083    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3084      iColumn=listNonLinearColumn[jNon];
3085      maxDelta = CoinMax(maxDelta,
3086                     fabs(solution[iColumn]-saveSolution[iColumn]));
3087      if (goodMove>0) {
3088        if (last[0][jNon]*last[1][jNon]<0) {
3089          // halve
3090          trust[jNon] *= 0.5;
3091          numberSmaller2++;
3092        } else {
3093          if (last[0][jNon]==last[1][jNon]&&
3094              last[0][jNon]==last[2][jNon])
3095            trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6); 
3096          numberLarger++;
3097        }
3098      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
3099        trust[jNon] *= 0.2;
3100        numberSmaller++;
3101      }
3102      maxGap = CoinMax(maxGap,trust[jNon]);
3103    }
3104#ifdef CLP_DEBUG
3105          if (handler_->logLevel()&32) 
3106            std::cout<<"largest gap is "<<maxGap<<" "
3107                     <<numberSmaller+numberSmaller2<<" reduced ("
3108                     <<numberSmaller<<" badMove ), "
3109                     <<numberLarger<<" increased"<<std::endl;
3110#endif
3111    if (iPass>10000) {
3112      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3113        trust[jNon] *=0.0001;
3114    }
3115    if(lastMoveStatus==-1&&goodMove==-1)
3116      goodMove=1; // to force solve
3117    if (goodMove>0) {
3118      double drop = lastObjective-objValue;
3119      handler_->message(CLP_SLP_ITER,messages_)
3120        <<iPass<<objValue-objectiveOffset
3121        <<drop<<maxDelta
3122        <<CoinMessageEol;
3123      if (iPass>20&&drop<1.0e-12*fabs(objValue))
3124        drop=0.999e-4; // so will exit
3125      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999) {
3126        if (handler_->logLevel()>1) 
3127          std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
3128        break;
3129      }
3130    } else if (!numberSmaller&&iPass>1) {
3131      if (handler_->logLevel()>1) 
3132          std::cout<<"Exiting as all gaps small"<<std::endl;
3133        break;
3134    }
3135    if (!iPass)
3136      goodMove=1;
3137    targetDrop=0.0;
3138    double * r = this->dualColumnSolution();
3139    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3140      iColumn=listNonLinearColumn[jNon];
3141      columnLower[iColumn]=CoinMax(solution[iColumn]
3142                               -trust[jNon],
3143                               trueLower[jNon]);
3144      columnUpper[iColumn]=CoinMin(solution[iColumn]
3145                               +trust[jNon],
3146                               trueUpper[jNon]);
3147    }
3148    if (iPass) {
3149      // get reduced costs
3150      this->matrix()->transposeTimes(savePi,
3151                                     this->dualColumnSolution());
3152      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3153        iColumn=listNonLinearColumn[jNon];
3154        double dj = objective[iColumn]-r[iColumn];
3155        r[iColumn]=dj;
3156        if (dj<-dualTolerance_) 
3157          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
3158        else if (dj>dualTolerance_)
3159          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
3160      }
3161    } else {
3162      memset(r,0,numberColumns*sizeof(double));
3163    }
3164#if 0
3165    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3166      iColumn=listNonLinearColumn[jNon];
3167      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
3168        columnLower[iColumn]=CoinMax(solution[iColumn],
3169                                 trueLower[jNon]);
3170        columnUpper[iColumn]=CoinMin(solution[iColumn]
3171                                 +trust[jNon],
3172                                 trueUpper[jNon]);
3173      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
3174        columnLower[iColumn]=CoinMax(solution[iColumn]
3175                                 -trust[jNon],
3176                                 trueLower[jNon]);
3177        columnUpper[iColumn]=CoinMin(solution[iColumn],
3178                                 trueUpper[jNon]);
3179      } else {
3180        columnLower[iColumn]=CoinMax(solution[iColumn]
3181                                 -trust[jNon],
3182                                 trueLower[jNon]);
3183        columnUpper[iColumn]=CoinMin(solution[iColumn]
3184                                 +trust[jNon],
3185                                 trueUpper[jNon]);
3186      }
3187    }
3188#endif
3189    if (goodMove>0) {
3190      memcpy(saveSolution,solution,numberColumns*sizeof(double));
3191      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
3192      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
3193      memcpy(saveStatus,status_,numberRows+numberColumns);
3194#if MULTIPLE>2
3195      double * tempSol=saveSolutionM[0];
3196      for (jNon=0;jNon<MULTIPLE-1;jNon++) {
3197        saveSolutionM[jNon]=saveSolutionM[jNon+1];
3198      }
3199      saveSolutionM[MULTIPLE-1]=tempSol;
3200      memcpy(tempSol,solution,numberColumns*sizeof(double));
3201     
3202#endif
3203     
3204#ifdef CLP_DEBUG
3205      if (handler_->logLevel()&32) 
3206        std::cout<<"Pass - "<<iPass
3207                 <<", target drop is "<<targetDrop
3208                 <<std::endl;
3209#endif
3210      lastObjective = objValue;
3211      if (targetDrop<CoinMax(1.0e-8,CoinMin(1.0e-6,1.0e-6*fabs(objValue)))&&goodMove&&iPass>3) {
3212        if (handler_->logLevel()>1) 
3213        printf("Exiting on target drop %g\n",targetDrop);
3214        break;
3215      }
3216#ifdef CLP_DEBUG
3217      {
3218        double * r = this->dualColumnSolution();
3219        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3220          iColumn=listNonLinearColumn[jNon];
3221          if (handler_->logLevel()&32) 
3222            printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
3223                   jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
3224                   r[iColumn],statusCheck[iColumn],columnLower[iColumn],
3225                   columnUpper[iColumn]);
3226        }
3227      }
3228#endif
3229      //setLogLevel(63);
3230      this->scaling(false);
3231      if (saveLogLevel==1)
3232        setLogLevel(0);
3233      ClpSimplex::primal(1);
3234      algorithm_=1;
3235      setLogLevel(saveLogLevel);
3236#ifdef CLP_DEBUG
3237      if (this->status()) {
3238        writeMps("xx.mps");
3239      }
3240#endif
3241      if (this->status()==1) {
3242        // not feasible ! - backtrack and exit
3243        // use safe solution
3244        memcpy(solution,safeSolution,numberColumns*sizeof(double));
3245        memcpy(saveSolution,solution,numberColumns*sizeof(double));
3246        memset(rowActivity_,0,numberRows_*sizeof(double));
3247        times(1.0,solution,rowActivity_);
3248        memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
3249        memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
3250        memcpy(status_,saveStatus,numberRows+numberColumns);
3251        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3252          iColumn=listNonLinearColumn[jNon];
3253          columnLower[iColumn]=CoinMax(solution[iColumn]
3254                                       -trust[jNon],
3255                                       trueLower[jNon]);
3256          columnUpper[iColumn]=CoinMin(solution[iColumn]
3257                                       +trust[jNon],
3258                                       trueUpper[jNon]);
3259        }
3260        break;
3261      } else {
3262        // save in case problems
3263        memcpy(safeSolution,solution,numberColumns*sizeof(double));
3264      }
3265      if (problemStatus_==3)
3266        break; // probably user interrupt
3267      goodMove=1;
3268    } else {
3269      // bad pass - restore solution
3270#ifdef CLP_DEBUG
3271      if (handler_->logLevel()&32) 
3272        printf("Backtracking\n");
3273#endif
3274      memcpy(solution,saveSolution,numberColumns*sizeof(double));
3275      memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
3276      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
3277      memcpy(status_,saveStatus,numberRows+numberColumns);
3278      iPass--;
3279      assert (exitPass>0);
3280      goodMove=-1;
3281    }
3282  }
3283#if MULTIPLE>2
3284  for (jNon=0;jNon<MULTIPLE;jNon++) {
3285    delete [] saveSolutionM[jNon];
3286  }
3287#endif
3288  // restore solution
3289  memcpy(solution,saveSolution,numberColumns*sizeof(double));
3290  memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
3291  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3292    iColumn=listNonLinearColumn[jNon];
3293    columnLower[iColumn]=CoinMax(solution[iColumn],
3294                             trueLower[jNon]);
3295    columnUpper[iColumn]=CoinMin(solution[iColumn],
3296                             trueUpper[jNon]);
3297  }
3298  delete [] markNonlinear;
3299  delete [] statusCheck;
3300  delete [] savePi;
3301  delete [] saveStatus;
3302  // redo objective
3303  memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3304         numberColumns*sizeof(double));
3305  ClpSimplex::primal(1);
3306  delete objective_;
3307  objective_=trueObjective;
3308  // redo values
3309  setDblParam(ClpObjOffset,objectiveOffset);
3310  objectiveValue_ += whichWay*offset;
3311  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3312    iColumn=listNonLinearColumn[jNon];
3313    columnLower[iColumn]= trueLower[jNon];
3314    columnUpper[iColumn]= trueUpper[jNon];
3315  }
3316  delete [] saveSolution;
3317  delete [] safeSolution;
3318  delete [] saveRowSolution;
3319  for (iPass=0;iPass<3;iPass++) 
3320    delete [] last[iPass];
3321  delete [] trust;
3322  delete [] trueUpper;
3323  delete [] trueLower;
3324  delete [] listNonLinearColumn;
3325  delete [] changeRegion;
3326  // temp
3327  //setLogLevel(63);
3328  return 0;
3329}
3330/* Primal algorithm for nonlinear constraints
3331   Using a semi-trust region approach as for pooling problem
3332   This is in because I have it lying around
3333   
3334*/
3335int 
3336ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
3337                               int numberPasses, double deltaTolerance)
3338{
3339  if (!numberConstraints) {
3340    // no nonlinear constraints - may be nonlinear objective
3341    return primalSLP(numberPasses,deltaTolerance);
3342  }
3343  // Are we minimizing or maximizing
3344  double whichWay=optimizationDirection();
3345  if (whichWay<0.0)
3346    whichWay=-1.0;
3347  else if (whichWay>0.0)
3348    whichWay=1.0;
3349  // check all matrix for odd rows is empty
3350  int iConstraint;
3351  int numberBad=0;
3352  CoinPackedMatrix * columnCopy = matrix();
3353  // Get a row copy in standard format
3354  CoinPackedMatrix copy;
3355  copy.reverseOrderedCopyOf(*columnCopy);
3356  // get matrix data pointers
3357  //const int * column = copy.getIndices();
3358  //const CoinBigIndex * rowStart = copy.getVectorStarts();
3359  const int * rowLength = copy.getVectorLengths(); 
3360  //const double * elementByRow = copy.getElements();
3361  int numberArtificials=0;
3362  // We could use nonlinearcost to do segments - maybe later
3363#define SEGMENTS 3
3364  // Penalties may be adjusted by duals
3365  // Both these should be modified depending on problem
3366  // Possibly start with big bounds
3367  //double penalties[]={1.0e-3,1.0e7,1.0e9};
3368  double penalties[]={1.0e7,1.0e8,1.0e9};
3369  double bounds[] = {1.0e-2,1.0e2,COIN_DBL_MAX};
3370  // see how many extra we need
3371  CoinBigIndex numberExtra=0;
3372  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3373    ClpConstraint * constraint = constraints[iConstraint];
3374    int iRow = constraint->rowNumber();
3375    assert (iRow>=0);
3376    int n = constraint->numberCoefficients() -rowLength[iRow];
3377    numberExtra += n;
3378    if (iRow>=numberRows_)
3379      numberBad++;
3380    else if (rowLength[iRow]&&n)
3381      numberBad++;
3382    if (rowLower_[iRow]>-1.0e20)
3383      numberArtificials++;
3384    if (rowUpper_[iRow]<1.0e20)
3385      numberArtificials++;
3386  }
3387  if (numberBad)
3388    return numberBad;
3389  ClpObjective * trueObjective = NULL;
3390  if (objective_->type()>=2) {
3391    // Replace objective
3392    trueObjective = objective_;
3393    objective_=new ClpLinearObjective(NULL,numberColumns_);
3394  }
3395  ClpSimplex newModel(*this);
3396  // we can put back true objective
3397  if (trueObjective) {
3398    // Replace objective
3399    delete objective_;
3400    objective_=trueObjective;
3401  }
3402  int numberColumns2 = numberColumns_;
3403  int numberSmallGap = numberArtificials;
3404  if (numberArtificials) {
3405    numberArtificials *= SEGMENTS;
3406    numberColumns2 += numberArtificials;
3407    int * addStarts = new int [numberArtificials+1];
3408    int * addRow = new int[numberArtificials];
3409    double * addElement = new double[numberArtificials];
3410    double * addUpper = new double[numberArtificials];
3411    addStarts[0]=0;
3412    double * addCost = new double [numberArtificials];
3413    numberArtificials=0;
3414    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3415      ClpConstraint * constraint = constraints[iConstraint];
3416      int iRow = constraint->rowNumber();
3417      if (rowLower_[iRow]>-1.0e20) {
3418        for (int k=0;k<SEGMENTS;k++) {
3419          addRow[numberArtificials]=iRow;
3420          addElement[numberArtificials]=1.0;
3421          addCost[numberArtificials]=penalties[k];
3422          addUpper[numberArtificials]=bounds[k];
3423          numberArtificials++;
3424          addStarts[numberArtificials]=numberArtificials;
3425        }
3426      }
3427      if (rowUpper_[iRow]<1.0e20) {
3428        for (int k=0;k<SEGMENTS;k++) {
3429          addRow[numberArtificials]=iRow;
3430          addElement[numberArtificials]=-1.0;
3431          addCost[numberArtificials]=penalties[k];
3432          addUpper[numberArtificials]=bounds[k];
3433          numberArtificials++;
3434          addStarts[numberArtificials]=numberArtificials;
3435        }
3436      }
3437    }
3438    newModel.addColumns(numberArtificials,NULL,addUpper,addCost,
3439                       addStarts,addRow,addElement);
3440    delete [] addStarts;
3441    delete [] addRow;
3442    delete [] addElement;
3443    delete [] addUpper;
3444    delete [] addCost;
3445    //    newModel.primal(1);
3446  }
3447  // find nonlinear columns
3448  int * listNonLinearColumn = new int [numberColumns_+numberSmallGap];
3449  char * mark = new char[numberColumns_];
3450  memset(mark,0,numberColumns_);
3451  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3452    ClpConstraint * constraint = constraints[iConstraint];
3453    constraint->markNonlinear(mark);
3454  }
3455  if (trueObjective)
3456    trueObjective->markNonlinear(mark);
3457  int iColumn;
3458  int numberNonLinearColumns=0;
3459  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3460    if (mark[iColumn])
3461      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
3462  }
3463  // and small gap artificials
3464  for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
3465    listNonLinearColumn[numberNonLinearColumns++]=iColumn;
3466  }
3467  // build row copy of matrix with space for nonzeros
3468  // Get column copy
3469  columnCopy = newModel.matrix();
3470  copy.reverseOrderedCopyOf(*columnCopy);
3471  // get matrix data pointers
3472  const int * column = copy.getIndices();
3473  const CoinBigIndex * rowStart = copy.getVectorStarts();
3474  rowLength = copy.getVectorLengths(); 
3475  const double * elementByRow = copy.getElements();
3476  numberExtra +=copy.getNumElements();
3477  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
3478  int * newColumn = new int[numberExtra];
3479  double * newElement = new double[numberExtra];
3480  newStarts[0]=0;
3481  int * backRow = new int [numberRows_];
3482  int iRow;
3483  for (iRow=0;iRow<numberRows_;iRow++) 
3484    backRow[iRow]=-1;
3485  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3486    ClpConstraint * constraint = constraints[iConstraint];
3487    iRow = constraint->rowNumber();
3488    backRow[iRow]=iConstraint;
3489  }
3490  numberExtra=0;
3491  for (iRow=0;iRow<numberRows_;iRow++) {
3492    if (backRow[iRow]<0) {
3493      // copy normal
3494      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3495           j++) {
3496        newColumn[numberExtra]=column[j];
3497        newElement[numberExtra++] = elementByRow[j];
3498      }
3499    } else {
3500      ClpConstraint * constraint = constraints[backRow[iRow]];
3501      assert(iRow == constraint->rowNumber());
3502      int numberArtificials=0;
3503      if (rowLower_[iRow]>-1.0e20)
3504        numberArtificials += SEGMENTS;
3505      if (rowUpper_[iRow]<1.0e20)
3506        numberArtificials += SEGMENTS;
3507      if (numberArtificials==rowLength[iRow]) {
3508        // all possible
3509        memset(mark,0,numberColumns_);
3510        constraint->markNonzero(mark);
3511        for (int k=0;k<numberColumns_;k++) {
3512          if (mark[k]) {
3513            newColumn[numberExtra]=k;
3514            newElement[numberExtra++] = 1.0;
3515          }
3516        }
3517        // copy artificials
3518        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3519             j++) {
3520          newColumn[numberExtra]=column[j];
3521          newElement[numberExtra++] = elementByRow[j];
3522        }
3523      } else {
3524        // there already
3525        // copy
3526        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3527             j++) {
3528          newColumn[numberExtra]=column[j];
3529          assert (elementByRow[j]);
3530          newElement[numberExtra++] = elementByRow[j];
3531        }
3532      }
3533    }
3534    newStarts[iRow+1]=numberExtra;
3535  }
3536  delete [] backRow;
3537  CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_,
3538                              numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0);
3539  delete [] newStarts;
3540  delete [] newColumn;
3541  delete [] newElement;
3542  delete [] mark;
3543  // get feasible
3544  if (whichWay<0.0) {
3545    newModel.setOptimizationDirection(1.0);
3546    double * objective = newModel.objective();
3547    for (int iColumn=0;iColumn<numberColumns_;iColumn++) 
3548      objective[iColumn] = - objective[iColumn];
3549  }
3550  newModel.primal(1);
3551  // still infeasible
3552  if (newModel.problemStatus()==1) {
3553    delete [] listNonLinearColumn;
3554    return 0;
3555  } else if (newModel.problemStatus()==2) {
3556    // unbounded - add bounds
3557    double * columnLower = newModel.columnLower();
3558    double * columnUpper = newModel.columnUpper();
3559    for (int i=0;i<numberColumns_;i++) {
3560      columnLower[i]= CoinMax(-1.0e8,columnLower[i]);
3561      columnUpper[i]= CoinMin(1.0e8,columnUpper[i]);
3562    }
3563    newModel.primal(1);
3564  }
3565  int numberRows = newModel.numberRows();
3566  double * columnLower = newModel.columnLower();
3567  double * columnUpper = newModel.columnUpper();
3568  double * objective = newModel.objective();
3569  double * rowLower = newModel.rowLower();
3570  double * rowUpper = newModel.rowUpper(); 
3571  double * solution = newModel.primalColumnSolution();
3572  int jNon;
3573  int * last[3];
3574 
3575  double * trust = new double[numberNonLinearColumns];
3576  double * trueLower = new double[numberNonLinearColumns];
3577  double * trueUpper = new double[numberNonLinearColumns];
3578  double objectiveOffset;
3579  double objectiveOffset2;
3580  getDblParam(ClpObjOffset,objectiveOffset);
3581  objectiveOffset *= whichWay;
3582  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3583    iColumn=listNonLinearColumn[jNon];
3584    double upper = columnUpper[iColumn];
3585    double lower = columnLower[iColumn];
3586    if (solution[iColumn]<lower)
3587      solution[iColumn]=lower;
3588    else if (solution[iColumn]>upper)
3589      solution[iColumn]=upper;
3590#if 0
3591    double large = CoinMax(1000.0,10.0*fabs(solution[iColumn]));
3592    if (upper>1.0e10)
3593      upper = solution[iColumn]+large;
3594    if (lower<-1.0e10)
3595      lower = solution[iColumn]-large;
3596#else
3597    upper = solution[iColumn]+0.5;
3598    lower = solution[iColumn]-0.5;
3599#endif
3600    //columnUpper[iColumn]=upper;
3601    trust[jNon]=0.05*(1.0+upper-lower);
3602    trueLower[jNon]=columnLower[iColumn];
3603    //trueUpper[jNon]=upper;
3604    trueUpper[jNon]=columnUpper[iColumn];
3605  }
3606  bool tryFix=false;
3607  int iPass;
3608  double lastObjective=1.0e31;
3609  double lastGoodObjective=1.0e31;
3610  double * bestSolution = NULL;
3611  double * saveSolution = new double [numberColumns2+numberRows];
3612  char * saveStatus = new char [numberColumns2+numberRows];
3613  double targetDrop=1.0e31;
3614  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
3615  for (iPass=0;iPass<3;iPass++) {
3616    last[iPass]=new int[numberNonLinearColumns];
3617    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3618      last[iPass][jNon]=0;
3619  }
3620  int numberZeroPasses=0;
3621  bool zeroTargetDrop=false;
3622  double * gradient = new double [numberColumns_];
3623  bool goneFeasible=false;
3624  // keep sum of artificials
3625#define KEEP_SUM 5
3626  double sumArt[KEEP_SUM];
3627  for (jNon=0;jNon<KEEP_SUM;jNon++)
3628    sumArt[jNon]=COIN_DBL_MAX;
3629#define SMALL_FIX 0.0
3630  for (iPass=0;iPass<numberPasses;iPass++) {
3631    objectiveOffset2 = objectiveOffset;
3632    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3633      iColumn=listNonLinearColumn[jNon];
3634      if (solution[iColumn]<trueLower[jNon])
3635        solution[iColumn]=trueLower[jNon];
3636      else if (solution[iColumn]>trueUpper[jNon])
3637        solution[iColumn]=trueUpper[jNon];
3638      columnLower[iColumn]=CoinMax(solution[iColumn]
3639                                   -trust[jNon],
3640                                   trueLower[jNon]);
3641      if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
3642        columnLower[iColumn]=SMALL_FIX;
3643      columnUpper[iColumn]=CoinMin(solution[iColumn]
3644                                   +trust[jNon],
3645                                   trueUpper[jNon]);
3646      if (!trueLower[jNon]) 
3647        columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
3648                                       columnLower[iColumn]+SMALL_FIX);
3649      if (!trueLower[jNon]&&tryFix&&
3650          columnLower[iColumn]==SMALL_FIX&&
3651          columnUpper[iColumn]<3.0*SMALL_FIX) {
3652        columnLower[iColumn]=0.0;
3653        solution[iColumn]=0.0;
3654        columnUpper[iColumn]=0.0;
3655        printf("fixing %d\n",iColumn);
3656      }
3657    }
3658    // redo matrix
3659    double offset;
3660    CoinPackedMatrix newMatrix(saveMatrix);
3661    // get matrix data pointers
3662    column = newMatrix.getIndices();
3663    rowStart = newMatrix.getVectorStarts();
3664    rowLength = newMatrix.getVectorLengths(); 
3665    // make sure x updated
3666    if (numberConstraints)
3667      constraints[0]->newXValues();
3668    else
3669      trueObjective->newXValues();
3670    double * changeableElement = newMatrix.getMutableElements();
3671    if (trueObjective) {
3672      memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3673             numberColumns_*sizeof(double));
3674    } else {
3675      memcpy(objective,objective_->gradient(this,solution,offset,true,2),
3676             numberColumns_*sizeof(double));
3677    }
3678    if (whichWay<0.0) {
3679      for (int iColumn=0;iColumn<numberColumns_;iColumn++) 
3680        objective[iColumn] = - objective[iColumn];
3681    }
3682    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3683      ClpConstraint * constraint = constraints[iConstraint];
3684      int iRow = constraint->rowNumber();
3685      double functionValue;
3686#ifndef NDEBUG
3687      int numberErrors = 
3688#endif
3689      constraint->gradient(&newModel,solution,gradient,functionValue,offset);
3690      assert (!numberErrors);
3691      // double dualValue = newModel.dualRowSolution()[iRow];
3692      int numberCoefficients = constraint->numberCoefficients();
3693      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
3694        int iColumn = column[j];
3695        changeableElement[j] = gradient[iColumn];
3696        //objective[iColumn] -= dualValue*gradient[iColumn];
3697        gradient[iColumn]=0.0;
3698      }
3699      for (int k=0;k<numberColumns_;k++)
3700        assert (!gradient[k]);
3701      if (rowLower_[iRow]>-1.0e20)
3702        rowLower[iRow] = rowLower_[iRow] - offset;
3703      if (rowUpper_[iRow]<1.0e20)
3704        rowUpper[iRow] = rowUpper_[iRow] - offset;
3705    }
3706    // Replace matrix
3707    // Get a column copy in standard format
3708    CoinPackedMatrix * columnCopy = new CoinPackedMatrix();;
3709    columnCopy->reverseOrderedCopyOf(newMatrix);
3710    newModel.replaceMatrix(columnCopy,true);
3711    // solve
3712    newModel.primal(1);
3713    if (newModel.status()==1) { 
3714      // Infeasible!
3715      newModel.allSlackBasis();
3716      newModel.primal();
3717      newModel.writeMps("infeas.mps");
3718      assert(!newModel.status());
3719    }
3720    double sumInfeas=0.0;
3721    int numberInfeas=0;
3722    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
3723      if (solution[iColumn]>1.0e-8) {
3724        numberInfeas++;
3725        sumInfeas += solution[iColumn];
3726      }
3727    }
3728    printf("%d artificial infeasibilities - summing to %g\n",
3729           numberInfeas,sumInfeas);
3730    for (jNon=0;jNon<KEEP_SUM-1;jNon++)
3731      sumArt[jNon]=sumArt[jNon+1];
3732    sumArt[KEEP_SUM-1]=sumInfeas;
3733    if (sumInfeas>0.01&&sumInfeas*1.1>sumArt[0]&&penalties[1]<1.0e7) {
3734      // not doing very well - increase - be more sophisticated later
3735      lastObjective=COIN_DBL_MAX;
3736      for (jNon=0;jNon<KEEP_SUM;jNon++)
3737        sumArt[jNon]=COIN_DBL_MAX;
3738      //for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
3739      //objective[iColumn+1] *= 1.5;
3740      //}
3741      penalties[1] *= 1.5;
3742      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3743        if (trust[jNon]>0.1)
3744          trust[jNon] *= 2.0;
3745        else
3746          trust[jNon] = 0.1;
3747    }
3748    if (sumInfeas<0.001&&!goneFeasible) {
3749      goneFeasible=true;
3750      penalties[0]=1.0e-3;
3751      penalties[1]=1.0e6;
3752      printf("Got feasible\n");
3753    }
3754    double infValue=0.0;
3755    double objValue=0.0;
3756    // make sure x updated
3757    if (numberConstraints)
3758      constraints[0]->newXValues();
3759    else
3760      trueObjective->newXValues();
3761    if (trueObjective) {
3762      objValue=trueObjective->objectiveValue(this,solution);
3763      printf("objective offset %g\n",offset);
3764      objectiveOffset2 =objectiveOffset+offset; // ? sign
3765      newModel.setObjectiveOffset(objectiveOffset2);
3766    } else {
3767      objValue=objective_->objectiveValue(this,solution);
3768    }
3769    objValue *= whichWay;
3770    double infPenalty=0.0;
3771    // This penalty is for target drop
3772    double infPenalty2=0.0;
3773    //const int * row = columnCopy->getIndices();
3774    //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3775    //const int * columnLength = columnCopy->getVectorLengths();
3776    //const double * element = columnCopy->getElements();
3777    double * cost = newModel.objective();
3778    column = newMatrix.getIndices();
3779    rowStart = newMatrix.getVectorStarts();
3780    rowLength = newMatrix.getVectorLengths(); 
3781    elementByRow = newMatrix.getElements();
3782    int jColumn = numberColumns_;
3783    double objectiveAdjustment=0.0;
3784    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3785      ClpConstraint * constraint = constraints[iConstraint];
3786      int iRow = constraint->rowNumber();
3787      double functionValue=constraint->functionValue(this,solution);
3788      double dualValue = newModel.dualRowSolution()[iRow];
3789      if (numberConstraints<-50)
3790        printf("For row %d current value is %g (row activity %g) , dual is %g\n",iRow,functionValue,
3791               newModel.primalRowSolution()[iRow],
3792               dualValue);
3793      double movement = newModel.primalRowSolution()[iRow]+constraint->offset();
3794      movement = fabs((movement-functionValue)*dualValue);
3795      infPenalty2 += movement;
3796      double sumOfActivities=0.0;
3797      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3798        int iColumn = column[j];
3799        sumOfActivities += fabs(solution[iColumn]*elementByRow[j]); 
3800      }
3801      if (rowLower_[iRow]>-1.0e20) {
3802        if (functionValue<rowLower_[iRow]-1.0e-5) {
3803          double infeasibility = rowLower_[iRow]-functionValue;
3804          double thisPenalty=0.0;
3805          infValue += infeasibility;
3806          double boundMultiplier = 1.0;
3807          if (sumOfActivities<0.001)
3808            boundMultiplier = 0.1;
3809          else if (sumOfActivities>100.0)
3810            boundMultiplier=10.0;
3811          int k;
3812          assert (dualValue>=-1.0e-5);
3813          dualValue = CoinMax(dualValue,0.0);
3814          for ( k=0;k<SEGMENTS;k++) {
3815            if (infeasibility<=0)
3816              break;
3817            double thisPart = CoinMin(infeasibility,bounds[k]);
3818            thisPenalty += thisPart*cost[jColumn+k];
3819            infeasibility -= thisPart;
3820          }
3821          infeasibility = functionValue-rowUpper_[iRow];
3822          double newPenalty=0.0;
3823          for ( k=0;k<SEGMENTS;k++) {
3824            double thisPart = CoinMin(infeasibility,bounds[k]);
3825            cost[jColumn+k] = CoinMax(penalties[k],dualValue+1.0e-3);
3826            newPenalty += thisPart*cost[jColumn+k];
3827            infeasibility -= thisPart;
3828          }
3829          infPenalty += thisPenalty;
3830          objectiveAdjustment += CoinMax(0.0,newPenalty-thisPenalty);
3831        }
3832        jColumn += SEGMENTS;
3833      }
3834      if (rowUpper_[iRow]<1.0e20) {
3835        if (functionValue>rowUpper_[iRow]+1.0e-5) {
3836          double infeasibility = functionValue-rowUpper_[iRow];
3837          double thisPenalty=0.0;
3838          infValue += infeasibility;
3839          double boundMultiplier = 1.0;
3840          if (sumOfActivities<0.001)
3841            boundMultiplier = 0.1;
3842          else if (sumOfActivities>100.0)
3843            boundMultiplier=10.0;
3844          int k;
3845          dualValue = -dualValue;
3846          assert (dualValue>=-1.0e-5);
3847          dualValue = CoinMax(dualValue,0.0);
3848          for ( k=0;k<SEGMENTS;k++) {
3849            if (infeasibility<=0)
3850              break;
3851            double thisPart = CoinMin(infeasibility,bounds[k]);
3852            thisPenalty += thisPart*cost[jColumn+k];
3853            infeasibility -= thisPart;
3854          }
3855          infeasibility = functionValue-rowUpper_[iRow];
3856          double newPenalty=0.0;
3857          for ( k=0;k<SEGMENTS;k++) {
3858            double thisPart = CoinMin(infeasibility,bounds[k]);
3859            cost[jColumn+k] = CoinMax(penalties[k],dualValue+1.0e-3);
3860            newPenalty += thisPart*cost[jColumn+k];
3861            infeasibility -= thisPart;
3862          }
3863          infPenalty += thisPenalty;
3864          objectiveAdjustment += CoinMax(0.0,newPenalty-thisPenalty);
3865        }
3866        jColumn += SEGMENTS;
3867      }
3868    }
3869    // adjust last objective value
3870    lastObjective += objectiveAdjustment;
3871    if (infValue)
3872      printf("Sum infeasibilities %g - penalty %g ",infValue,infPenalty);
3873    if (objectiveOffset2)
3874      printf("offset2 %g ",objectiveOffset2);
3875    objValue -= objectiveOffset2;
3876    printf("True objective %g or maybe %g (with penalty %g) -pen2 %g %g\n",objValue,
3877           objValue+objectiveOffset2,objValue+objectiveOffset2+infPenalty,infPenalty2,penalties[1]);
3878    double useObjValue = objValue+objectiveOffset2+infPenalty;
3879    objValue += infPenalty+infPenalty2;
3880    objValue = useObjValue;
3881    if (iPass) {
3882      double drop = lastObjective-objValue;
3883      std::cout<<"True drop was "<<drop<<std::endl;
3884      if (drop<-0.05*fabs(objValue)-1.0e-4) {
3885        // pretty bad - go back and halve
3886        memcpy(solution,saveSolution,numberColumns2*sizeof(double));
3887        memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2,
3888               numberRows*sizeof(double));
3889        memcpy(newModel.statusArray(),saveStatus,
3890               numberColumns2+numberRows);
3891        for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3892          if (trust[jNon]>0.1)
3893            trust[jNon] *= 0.5;
3894          else
3895            trust[jNon] *= 0.9;
3896       
3897        printf("** Halving trust\n");
3898        objValue=lastObjective;
3899        continue;
3900      } else if ((iPass%25)==-1) {
3901        for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3902          trust[jNon] *= 2.0;
3903        printf("** Doubling trust\n");
3904      }
3905      int * temp=last[2];
3906      last[2]=last[1];
3907      last[1]=last[0];
3908      last[0]=temp;
3909      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3910        iColumn=listNonLinearColumn[jNon];
3911        double change = solution[iColumn]-saveSolution[iColumn];
3912        if (change<-1.0e-5) {
3913          if (fabs(change+trust[jNon])<1.0e-5) 
3914            temp[jNon]=-1;
3915          else
3916            temp[jNon]=-2;
3917        } else if(change>1.0e-5) {
3918          if (fabs(change-trust[jNon])<1.0e-5) 
3919            temp[jNon]=1;
3920          else
3921            temp[jNon]=2;
3922        } else {
3923          temp[jNon]=0;
3924        }
3925      } 
3926      double maxDelta=0.0;
3927      double smallestTrust=1.0e31;
3928      double smallestNonLinearGap=1.0e31;
3929      double smallestGap=1.0e31;
3930      bool increasing=false;
3931      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3932        double gap = columnUpper[iColumn]-columnLower[iColumn];
3933        assert (gap>=0.0);
3934        if (gap)
3935          smallestGap = CoinMin(smallestGap,gap);
3936      }
3937      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3938        iColumn=listNonLinearColumn[jNon];
3939        double gap = columnUpper[iColumn]-columnLower[iColumn];
3940        assert (gap>=0.0);
3941        if (gap) {
3942          smallestNonLinearGap = CoinMin(smallestNonLinearGap,gap);
3943          if (gap<1.0e-7&&iPass==1) {
3944            printf("Small gap %d %d %g %g %g\n",
3945                   jNon,iColumn,columnLower[iColumn],columnUpper[iColumn],
3946                   gap);
3947            //trueUpper[jNon]=trueLower[jNon];
3948            //columnUpper[iColumn]=columnLower[iColumn];
3949          }
3950        }
3951        maxDelta = CoinMax(maxDelta,
3952                       fabs(solution[iColumn]-saveSolution[iColumn]));
3953        if (last[0][jNon]*last[1][jNon]<0) {
3954          // halve
3955          if (trust[jNon]>1.0)
3956            trust[jNon] *= 0.5;
3957          else
3958            trust[jNon] *= 0.7;
3959        } else {
3960          // ? only increase if +=1 ?
3961          if (last[0][jNon]==last[1][jNon]&&
3962              last[0][jNon]==last[2][jNon]&&
3963              last[0][jNon]) {
3964            trust[jNon] *= 1.8; 
3965            increasing=true;
3966          }
3967        }
3968        smallestTrust = CoinMin(smallestTrust,trust[jNon]);
3969      }
3970      std::cout<<"largest delta is "<<maxDelta
3971               <<", smallest trust is "<<smallestTrust
3972               <<", smallest gap is "<<smallestGap
3973               <<", smallest nonlinear gap is "<<smallestNonLinearGap
3974               <<std::endl;
3975      if (iPass>200) {
3976        //double useObjValue = objValue+objectiveOffset2+infPenalty;
3977        if (useObjValue+1.0e-4> lastGoodObjective&&iPass>250) {
3978          std::cout<<"Exiting as objective not changing much"<<std::endl;
3979          break;
3980        } else if (useObjValue<lastGoodObjective) {
3981          lastGoodObjective = useObjValue;
3982          if (!bestSolution)
3983            bestSolution = new double [numberColumns2];
3984          memcpy(bestSolution,solution,numberColumns2*sizeof(double));
3985        }
3986      }
3987      if (maxDelta<deltaTolerance&&!increasing&&iPass>100) {
3988        numberZeroPasses++;
3989        if (numberZeroPasses==3) {
3990          if (tryFix) {
3991            std::cout<<"Exiting"<<std::endl;
3992            break;
3993          } else {
3994            tryFix=true;
3995            for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3996              trust[jNon] = CoinMax(trust[jNon],1.0e-1);
3997            numberZeroPasses=0;
3998          }
3999        }
4000      } else {
4001        numberZeroPasses=0;
4002      }
4003    }
4004    memcpy(saveSolution,solution,numberColumns2*sizeof(double));
4005    memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(),
4006           numberRows*sizeof(double));
4007    memcpy(saveStatus,newModel.statusArray(),
4008           numberColumns2+numberRows);
4009   
4010    targetDrop=infPenalty+infPenalty2;
4011    if (iPass) {
4012      // get reduced costs
4013      const double * pi = newModel.dualRowSolution();
4014      newModel.matrix()->transposeTimes(pi,
4015                                        newModel.dualColumnSolution());
4016      double * r = newModel.dualColumnSolution();
4017      for (iColumn=0;iColumn<numberColumns_;iColumn++) 
4018        r[iColumn] = objective[iColumn]-r[iColumn];
4019      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
4020        iColumn=listNonLinearColumn[jNon];
4021        double dj = r[iColumn];
4022        if (dj<-1.0e-6) {
4023          double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
4024          //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
4025          //double drop2 = -dj*(upper-solution[iColumn]);
4026#if 0
4027          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2&&iPass>100))
4028            printf("Big drop %d %g %g %g %g T %g\n",
4029                   iColumn,columnLower[iColumn],solution[iColumn],
4030                   columnUpper[iColumn],dj,trueUpper[jNon]);
4031#endif
4032          targetDrop += drop;
4033          if (dj<-1.0e-1&&trust[jNon]<1.0e-3
4034              &&trueUpper[jNon]-solution[iColumn]>1.0e-2) {
4035            trust[jNon] *= 1.5;
4036            //printf("Increasing trust on %d to %g\n",
4037            //     iColumn,trust[jNon]);
4038          }
4039        } else if (dj>1.0e-6) {
4040          double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
4041          //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
4042          //double drop2 = -dj*(lower-solution[iColumn]);
4043#if 0
4044          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2))
4045            printf("Big drop %d %g %g %g %g T %g\n",
4046                   iColumn,columnLower[iColumn],solution[iColumn],
4047                   columnUpper[iColumn],dj,trueLower[jNon]);
4048#endif
4049          targetDrop += drop;
4050          if (dj>1.0e-1&&trust[jNon]<1.0e-3
4051              &&solution[iColumn]-trueLower[jNon]>1.0e-2) {
4052            trust[jNon] *= 1.5;
4053            printf("Increasing trust on %d to %g\n",
4054                   iColumn,trust[jNon]);
4055          }
4056        }
4057      }
4058    }
4059    std::cout<<"Pass - "<<iPass
4060             <<", target drop is "<<targetDrop
4061             <<std::endl;
4062    if (iPass>1&&targetDrop<1.0e-5&&zeroTargetDrop)
4063      break;
4064    if (iPass>1&&targetDrop<1.0e-5)
4065      zeroTargetDrop = true;
4066    else
4067      zeroTargetDrop = false;
4068    //if (iPass==5)
4069    //newModel.setLogLevel(63);
4070    lastObjective = objValue;
4071    // take out when ClpPackedMatrix changed
4072    //newModel.scaling(false);
4073#if 0
4074    CoinMpsIO writer;
4075    writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
4076                      newModel.getColLower(), newModel.getColUpper(),
4077                      newModel.getObjCoefficients(),
4078                      (const char*) 0 /*integrality*/,
4079                      newModel.getRowLower(), newModel.getRowUpper(),
4080                      NULL,NULL);
4081    writer.writeMps("xx.mps");
4082#endif
4083  }
4084  delete [] saveSolution;
4085  delete [] saveStatus;
4086  for (iPass=0;iPass<3;iPass++) 
4087    delete [] last[iPass];
4088  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
4089    iColumn=listNonLinearColumn[jNon];
4090    columnLower[iColumn]=trueLower[jNon];
4091    columnUpper[iColumn]=trueUpper[jNon];
4092  }
4093  delete [] trust;
4094  delete [] trueUpper;
4095  delete [] trueLower;
4096  objectiveValue_=newModel.objectiveValue();
4097  if (bestSolution) {
4098    memcpy(solution,bestSolution,numberColumns2*sizeof(double));
4099    delete [] bestSolution;
4100    printf("restoring objective of %g\n",lastGoodObjective);
4101    objectiveValue_=lastGoodObjective;
4102  }
4103  // Simplest way to get true row activity ?
4104  double * rowActivity = newModel.primalRowSolution();
4105  for (iRow=0;iRow<numberRows;iRow++) {
4106    double difference;
4107    if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow]))
4108      difference = rowLower_[iRow]-rowLower[iRow];
4109    else
4110      difference = rowUpper_[iRow]-rowUpper[iRow];
4111    rowLower[iRow]=rowLower_[iRow];
4112    rowUpper[iRow]=rowUpper_[iRow];
4113    if (difference) {
4114      if (numberRows<50)
4115        printf("For row %d activity changes from %g to %g\n",
4116               iRow,rowActivity[iRow],rowActivity[iRow]+difference);
4117      rowActivity[iRow]+= difference;
4118    }
4119  }
4120  delete [] listNonLinearColumn;
4121  delete [] gradient;
4122  printf("solution still in newModel - do objective etc!\n");
4123  numberIterations_ =newModel.numberIterations();
4124  problemStatus_ =newModel.problemStatus();
4125  secondaryStatus_ =newModel.secondaryStatus();
4126  memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double));
4127  // should do status region
4128  CoinZeroN(rowActivity_,numberRows_);
4129  matrix_->times(1.0,columnActivity_,rowActivity_);
4130  return 0;
4131}
Note: See TracBrowser for help on using the repository browser.