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

Last change on this file since 1094 was 1094, checked in by forrest, 14 years ago

stop loop in quadratic

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