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

Last change on this file since 1366 was 1366, checked in by forrest, 11 years ago

try and improve quadratic and barrier

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