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

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

out compiler warnings and stability improvements

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