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

Last change on this file since 1034 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 115.9 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)
1302        break;
1303      if (!nPasses) {
1304        numberNonBasic=0;
1305        nBigPasses++;
1306      }
1307      // just do superbasic if in cleanup mode
1308      int local=localPivotMode;
1309      if (!local&&pivotMode>=10&&nBigPasses<10) {
1310        local=100;
1311      } else if (local==-1||nBigPasses>=10) {
1312        local=0;
1313        localPivotMode=0;
1314      } 
1315      if (justOne) {
1316        local=2;
1317        //local=100;
1318        justOne=false;
1319      }
1320      if (!nPasses)
1321        startLocalMode=local;
1322      directionVector(longArray,spare,rowArray,local,
1323                      normFlagged,normUnflagged,numberNonBasic);
1324      {
1325        // check null vector
1326        double * rhs = spare->denseVector();
1327#if CLP_DEBUG > 1
1328        spare->checkClear();
1329#endif
1330        int iRow;
1331        multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,rhs,0.0);
1332        matrix_->times(1.0,solution_,rhs,rowScale_,columnScale_);
1333        double largest=0.0;
1334        int iLargest=-1;
1335        for (iRow=0;iRow<numberRows_;iRow++) {
1336          double value=fabs(rhs[iRow]);
1337          rhs[iRow]=0.0;
1338          if (value>largest) {
1339            largest=value;
1340            iLargest=iRow;
1341          }
1342        }
1343#if CLP_DEBUG > 0
1344        if ((handler_->logLevel()&32)&&largest>1.0e-8) 
1345          printf("largest rhs error %g on row %d\n",largest,iLargest);
1346#endif
1347        if (solutionError<0.0) {
1348          solutionError=largest;
1349        } else if (largest>max(1.0e-8,1.0e2*solutionError)&&
1350                   factorization_->pivots()) {
1351          longArray->clear();
1352          pivotRow_ = -1;
1353          theta_=0.0;
1354          return 3;
1355        }
1356      }
1357      if (sequenceIn_>=0)
1358        lastSequenceIn=sequenceIn_;
1359#if MINTYPE!=2
1360      double djNormSave = djNorm;
1361#endif
1362      djNorm=0.0;
1363      int iIndex;
1364      for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
1365        int iSequence = which[iIndex];
1366        double alpha = work[iSequence];
1367        djNorm += alpha*alpha;
1368      }
1369      // go to conjugate gradient if necessary
1370      if (numberNonBasic&&localPivotMode>=10&&(nPasses>4||sequenceIn_<0)) {
1371        localPivotMode=0;
1372        nTotalPasses+= nPasses;
1373        nPasses=0;
1374      }
1375#if CONJUGATE == 2
1376      conjugate = (localPivotMode<10) ? MINTYPE : 0;
1377#endif
1378      //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
1379      //     nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
1380      if (!nPasses) {
1381        //printf("numberNon %d\n",numberNonBasic);
1382#if MINTYPE==2
1383        assert (numberNonBasic<100);
1384        memset(zz,0,numberNonBasic*numberNonBasic*sizeof(double));
1385        int put=0;
1386        for (int iVariable=0;iVariable<numberNonBasic;iVariable++) {
1387          zz[put]=1.0;
1388          put+= numberNonBasic+1;
1389        }
1390#endif
1391        djNorm0 = CoinMax(djNorm, 1.0e-20);
1392        memcpy(dArray,work,numberTotal*sizeof(double));
1393        memcpy(dArray2,work,numberTotal*sizeof(double));
1394        if (sequenceIn_>=0&&numberNonBasic==1) {
1395          // see if simple move
1396          double objTheta2 =objective_->stepLength(this,solution_,work,1.0e30,
1397                                               currentObj,predictedObj,thetaObj);
1398          rowArray->clear();
1399         
1400          // update the incoming column
1401          unpackPacked(rowArray);
1402          factorization_->updateColumnFT(spare,rowArray);
1403          theta_=0.0;
1404          double * work2=rowArray->denseVector();
1405          int number=rowArray->getNumElements();
1406          int * which2=rowArray->getIndices();
1407          int iIndex;
1408          bool easyMove=false;
1409          double way;
1410          if (dj_[sequenceIn_]>0.0)
1411            way=-1.0;
1412          else
1413            way=1.0;
1414          double largest=COIN_DBL_MAX;
1415          int kPivot=-1;
1416          for (iIndex=0;iIndex<number;iIndex++) {
1417            int iRow = which2[iIndex];
1418            double alpha = way*work2[iIndex];
1419            int iPivot = pivotVariable_[iRow];
1420            if (alpha<-1.0e-5) {
1421              double distance = upper_[iPivot]-solution_[iPivot];
1422              if (distance<-largest*alpha) {
1423                kPivot=iPivot;
1424                largest=CoinMax(0.0,-distance/alpha);
1425              }
1426              if (distance<-1.0e-12*alpha) {
1427                easyMove=true;
1428                break;
1429              }
1430            } else if (alpha>1.0e-5) {
1431              double distance = solution_[iPivot]-lower_[iPivot];
1432              if (distance<largest*alpha) {
1433                kPivot=iPivot;
1434                largest=CoinMax(0.0,distance/alpha);
1435              }
1436              if (distance<1.0e-12*alpha) {
1437                easyMove=true;
1438                break;
1439              }
1440            }
1441          }
1442          // But largest has to match up with change
1443          assert (work[sequenceIn_]);
1444          largest /= fabs(work[sequenceIn_]);
1445          if (largest<objTheta2) {
1446            easyMove=true;
1447          } else if (!easyMove) {
1448            double objDrop = currentObj-predictedObj;
1449            double th = objective_->stepLength(this,solution_,work,largest,
1450                                   currentObj,predictedObj,simpleObjective);
1451            simpleObjective = CoinMax(simpleObjective,predictedObj);
1452            double easyDrop = currentObj-simpleObjective;
1453            if (easyDrop>1.0e-8&&easyDrop>0.5*objDrop) {
1454              easyMove=true;
1455#ifdef CLP_DEBUG
1456              if (handler_->logLevel()&32) 
1457                printf("easy - obj drop %g, easy drop %g\n",objDrop,easyDrop);
1458#endif
1459              if (easyDrop>objDrop) {
1460                // debug
1461                printf("****** th %g simple %g\n",th,simpleObjective);
1462                objective_->stepLength(this,solution_,work,1.0e30,
1463                                   currentObj,predictedObj,simpleObjective);
1464                objective_->stepLength(this,solution_,work,largest,
1465                                   currentObj,predictedObj,simpleObjective);
1466              }
1467            }
1468          }
1469          rowArray->clear();
1470#ifdef CLP_DEBUG
1471          if (handler_->logLevel()&32)
1472            printf("largest %g piv %d\n",largest,kPivot);
1473#endif
1474          if (easyMove) {
1475            valueIn_=solution_[sequenceIn_];
1476            dualIn_=dj_[sequenceIn_];
1477            lowerIn_=lower_[sequenceIn_];
1478            upperIn_=upper_[sequenceIn_];
1479            if (dualIn_>0.0)
1480              directionIn_ = -1;
1481            else 
1482              directionIn_ = 1;
1483            longArray->clear();
1484            pivotRow_ = -1;
1485            theta_=0.0;
1486            return 0;
1487          }
1488        }
1489      } else {
1490#if MINTYPE==1
1491        if (conjugate) {
1492          double djNorm2 = djNorm;
1493          if (numberNonBasic&&0) {
1494            int iIndex;
1495            djNorm2 = 0.0;
1496            for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
1497              int iSequence = which[iIndex];
1498              double alpha = work[iSequence];
1499              //djNorm2 += alpha*alpha;
1500              double alpha2 = work[iSequence]-dArray2[iSequence];
1501              djNorm2 += alpha*alpha2;
1502            }
1503            //printf("a %.18g b %.18g\n",djNorm,djNorm2);
1504          }
1505          djNorm=djNorm2;
1506          double beta = djNorm2/djNormSave;
1507          // reset beta every so often
1508          //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1)
1509          //beta=0.0;
1510          //printf("current norm %g, old %g - beta %g\n",
1511          //     djNorm,djNormSave,beta);
1512          for (iSequence=0;iSequence<numberTotal;iSequence++) {
1513            dArray[iSequence] = work[iSequence] + beta *dArray[iSequence];
1514            dArray2[iSequence] = work[iSequence];
1515          }
1516        } else {
1517          for (iSequence=0;iSequence<numberTotal;iSequence++)
1518            dArray[iSequence] = work[iSequence];
1519        }
1520#else
1521        int number2=numberNonBasic;
1522        if (number2) {
1523          // pack down into dArray
1524          int jLast=-1;
1525          for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1526            int j=which[iSequence];
1527            assert(j>jLast);
1528            jLast=j;
1529            double value=work[j];
1530            dArray[iSequence]=-value;
1531          }
1532          // see whether to restart
1533          // check signs - gradient
1534          double g1=innerProduct(dArray,number2,dArray);
1535          double g2=innerProduct(dArray,number2,saveY2);
1536          // Get differences
1537          for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1538            saveY2[iSequence] = dArray[iSequence]-saveY2[iSequence];
1539            saveS2[iSequence] = solution_[iSequence]-saveS2[iSequence];
1540          }
1541          double g3=innerProduct(saveS2,number2,saveY2);
1542          printf("inner %g\n",g3);
1543          //assert(g3>0);
1544          double zzz[10000];
1545          int iVariable;
1546          g2=1.0e50;// temp
1547          if (fabs(g2)>=0.2*fabs(g1)) {
1548            // restart
1549            double delta = innerProduct(saveY2,number2,saveS2)/
1550              innerProduct(saveY2,number2,saveY2);
1551            delta=1.0;//temp
1552            memset(zz,0,number2*sizeof(double));
1553            int put=0;
1554            for (iVariable=0;iVariable<number2;iVariable++) {
1555              zz[put]=delta;
1556              put+= number2+1;
1557            }
1558          } else {
1559          }
1560          memcpy(zzz,zz,number2*number2*sizeof(double));
1561          double ww[100];
1562          // get sk -Hkyk
1563          for (iVariable=0;iVariable<number2;iVariable++) {
1564            double value=0.0;
1565            for (int jVariable=0;jVariable<number2;jVariable++) {
1566              value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
1567            }
1568            ww[iVariable] = saveS2[iVariable]-value;
1569          }
1570          double ys = innerProduct(saveY2,number2,saveS2);
1571          double multiplier1 = 1.0/ys;
1572          double multiplier2= innerProduct(saveY2,number2,ww)/(ys*ys);
1573#if 1
1574          // and second way
1575          // Hy
1576          double h[100];
1577          for (iVariable=0;iVariable<number2;iVariable++) {
1578            double value=0.0;
1579            for (int jVariable=0;jVariable<number2;jVariable++) {
1580              value += saveY2[jVariable]*zzz[iVariable+number2*jVariable];
1581            }
1582            h[iVariable] = value;
1583          }
1584          double hh[10000];
1585          double yhy1 = innerProduct(h,number2,saveY2)*multiplier1+1.0;
1586          yhy1 *= multiplier1;
1587          for (iVariable=0;iVariable<number2;iVariable++) {
1588            for (int jVariable=0;jVariable<number2;jVariable++) {
1589              int put = iVariable+number2*jVariable;
1590              double value = zzz[put];
1591              value += yhy1*saveS2[iVariable]*saveS2[jVariable];
1592              hh[put]=value;
1593            }
1594          }
1595          for (iVariable=0;iVariable<number2;iVariable++) {
1596            for (int jVariable=0;jVariable<number2;jVariable++) {
1597              int put = iVariable+number2*jVariable;
1598              double value = hh[put];
1599              value -= multiplier1*(saveS2[iVariable]*h[jVariable]);
1600              value -= multiplier1*(saveS2[jVariable]*h[iVariable]);
1601              hh[put]=value;
1602            }
1603          }
1604#endif
1605          // now update H
1606          for (iVariable=0;iVariable<number2;iVariable++) {
1607            for (int jVariable=0;jVariable<number2;jVariable++) {
1608              int put = iVariable+number2*jVariable;
1609              double value = zzz[put];
1610              value += multiplier1*(ww[iVariable]*saveS2[jVariable]
1611                                    +ww[jVariable]*saveS2[iVariable]);
1612              value -= multiplier2*saveS2[iVariable]*saveS2[jVariable];
1613              zzz[put]=value;
1614            }
1615          }
1616          //memcpy(zzz,hh,size*sizeof(double));
1617          // do search direction
1618          memset(dArray,0,numberTotal*sizeof(double));
1619          for (iVariable=0;iVariable<numberNonBasic;iVariable++) {
1620            double value=0.0;
1621            for (int jVariable=0;jVariable<number2;jVariable++) {
1622              int k=which[jVariable];
1623              value += work[k]*zzz[iVariable+number2*jVariable];
1624            }
1625            int i=which[iVariable];
1626            dArray[i] = value;
1627          }
1628          // Now fill out dArray
1629          {
1630            int j;
1631            // Now do basic ones
1632            int iRow;
1633            CoinIndexedVector * spare1=spare;
1634            CoinIndexedVector * spare2=rowArray;
1635#if CLP_DEBUG > 1
1636            spare1->checkClear();
1637            spare2->checkClear();
1638#endif
1639            double * array2=spare1->denseVector();
1640            int * index2 = spare1->getIndices();
1641            int number2=0;
1642            times(-1.0,dArray,array2);
1643            dArray = dArray+numberColumns_;
1644            for (iRow=0;iRow<numberRows_;iRow++) {
1645              double value = array2[iRow] + dArray[iRow];
1646              if (value) {
1647                array2[iRow]=value;
1648                index2[number2++]=iRow;
1649              } else {
1650                array2[iRow]=0.0;
1651              }
1652            }
1653            dArray -= numberColumns_;
1654            spare1->setNumElements(number2);
1655            // Ftran
1656            factorization_->updateColumn(spare2,spare1);
1657            number2=spare1->getNumElements();
1658            for (j=0;j<number2;j++) {
1659              int iSequence = index2[j];
1660              double value = array2[iSequence];
1661              array2[iSequence]=0.0;
1662              if (value) {
1663                int iPivot=pivotVariable_[iSequence];
1664                double oldValue = dArray[iPivot];
1665                dArray[iPivot] = value+oldValue;
1666              }
1667            }
1668            spare1->setNumElements(0);
1669          }
1670          //assert (innerProductIndexed(dArray,number2,work,which)>0);
1671          //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
1672          printf ("innerP1 %g innerP2 %g\n",innerProduct(dArray,numberTotal,work),
1673                  innerProductIndexed(dArray,numberNonBasic,work,which));
1674          assert (innerProduct(dArray,numberTotal,work)>0);
1675#if 1
1676          {
1677            // check null vector
1678            int iRow;
1679            double qq[10];
1680            memset(qq,0,numberRows_*sizeof(double));
1681            multiplyAdd(dArray+numberColumns_,numberRows_,-1.0,qq,0.0);
1682            matrix_->times(1.0,dArray,qq,rowScale_,columnScale_);
1683            double largest=0.0;
1684            int iLargest=-1;
1685            for (iRow=0;iRow<numberRows_;iRow++) {
1686              double value=fabs(qq[iRow]);
1687              if (value>largest) {
1688                largest=value;
1689                iLargest=iRow;
1690              }
1691            }
1692            printf("largest null error %g on row %d\n",largest,iLargest);
1693            for (iSequence=0;iSequence<numberTotal;iSequence++) {
1694              if (getStatus(iSequence)==basic)
1695                assert (fabs(dj_[iSequence])<1.0e-3);
1696            }
1697          }
1698#endif
1699          memcpy(saveY,saveY2,numberNonBasic*sizeof(double));
1700          memcpy(saveS,saveS2,numberNonBasic*sizeof(double));
1701        }
1702#endif
1703      }
1704#if MINTYPE==2
1705      for (iSequence=0;iSequence<numberNonBasic;iSequence++) {
1706        int j=which[iSequence];
1707        saveY2[iSequence]=-work[j];
1708        saveS2[iSequence]=solution_[j];
1709      }
1710#endif
1711      if (djNorm<eps*djNorm0||(nPasses>100&&djNorm<CoinMin(1.0e-1*djNorm0, 1.0e-12))) {
1712#ifdef CLP_DEBUG
1713        if (handler_->logLevel()&32) 
1714          printf("dj norm reduced from %g to %g\n",djNorm0,djNorm);
1715#endif
1716        if (pivotMode<10||!numberNonBasic) {
1717          finished=true;
1718        } else {
1719          localPivotMode=pivotMode;
1720          nTotalPasses+= nPasses;
1721          nPasses=0;
1722          continue;
1723        }
1724      }
1725      //if (nPasses==100||nPasses==50)
1726      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
1727      //         normFlagged);
1728      if (nPasses>100&&djNorm<1.0e-2*normFlagged&&!startLocalMode) {
1729#ifdef CLP_DEBUG
1730        if (handler_->logLevel()&32) 
1731          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n",djNorm0,djNorm,
1732                 normFlagged);
1733#endif
1734        unflag();
1735        localPivotMode=0;
1736        nTotalPasses+= nPasses;
1737        nPasses=0;
1738        continue;
1739      }
1740      if (djNorm>0.99*djNorm0&&nPasses>1500) {
1741        finished=true;
1742#ifdef CLP_DEBUG
1743        if (handler_->logLevel()&32) 
1744          printf("dj norm NOT reduced from %g to %g\n",djNorm0,djNorm);
1745#endif
1746        djNorm=1.2345e-20;
1747      }
1748      number=longArray->getNumElements();
1749      if (!numberNonBasic) {
1750        // looks optimal
1751        // check if any dj look plausible
1752        int nSuper=0;
1753        int nFlagged=0;
1754        int nNormal=0;
1755        for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
1756          if (flagged(iSequence)) {
1757            switch(getStatus(iSequence)) {
1758             
1759            case basic:
1760            case ClpSimplex::isFixed:
1761              break;
1762            case atUpperBound:
1763              if (dj_[iSequence]>dualTolerance_) 
1764                nFlagged++;
1765              break;
1766            case atLowerBound:
1767              if (dj_[iSequence]<-dualTolerance_) 
1768                nFlagged++;
1769              break;
1770            case isFree:
1771            case superBasic:
1772              if (fabs(dj_[iSequence])>dualTolerance_) 
1773                nFlagged++;
1774              break;
1775            }
1776            continue;
1777          }
1778          switch(getStatus(iSequence)) {
1779           
1780          case basic:
1781          case ClpSimplex::isFixed:
1782            break;
1783          case atUpperBound:
1784            if (dj_[iSequence]>dualTolerance_) 
1785              nNormal++;
1786            break;
1787          case atLowerBound:
1788            if (dj_[iSequence]<-dualTolerance_) 
1789              nNormal++;
1790            break;
1791          case isFree:
1792          case superBasic:
1793            if (fabs(dj_[iSequence])>dualTolerance_) 
1794              nSuper++;
1795            break;
1796          }
1797        }
1798#ifdef CLP_DEBUG
1799        if (handler_->logLevel()&32) 
1800          printf("%d super, %d normal, %d flagged\n",
1801                 nSuper,nNormal,nFlagged);
1802#endif
1803       
1804        int nFlagged2=1;
1805        if (lastSequenceIn<0&&!nNormal&&!nSuper) {
1806          nFlagged2=unflag();
1807          if (pivotMode>=10) {
1808            pivotMode--;
1809#ifdef CLP_DEBUG
1810            if (handler_->logLevel()&32) 
1811              printf("pivot mode now %d\n",pivotMode);
1812#endif
1813            if (pivotMode==9) 
1814              pivotMode=0;      // switch off fast attempt
1815          }
1816        } else {
1817          lastSequenceIn=-1;
1818        }
1819        if (!nFlagged2||(normFlagged+normUnflagged<1.0e-8)) {
1820          primalColumnPivot_->setLooksOptimal(true);
1821          return 0;
1822        } else {
1823          localPivotMode=-1;
1824          nTotalPasses+= nPasses;
1825          nPasses=0;
1826          finished=false;
1827          continue;
1828        }
1829      }
1830      bestSequence=-1;
1831      bestBasicSequence=-1;
1832     
1833      // temp
1834      nPasses++;
1835      if (nPasses>2000)
1836        finished=true;
1837      double theta = 1.0e30;
1838      basicTheta = 1.0e30;
1839      theta_=1.0e30;
1840      double basicTolerance = 1.0e-4*primalTolerance_;
1841      for (iSequence=0;iSequence<numberTotal;iSequence++) {
1842        //if (flagged(iSequence)
1843        //  continue;
1844        double alpha = dArray[iSequence];
1845        Status thisStatus = getStatus(iSequence);
1846        double oldValue = solution_[iSequence];
1847        if (thisStatus!=basic) {
1848          if (fabs(alpha)>=acceptablePivot) {
1849            if (alpha<0.0) {
1850              // variable going towards lower bound
1851              double bound = lower_[iSequence];
1852              oldValue -= bound;
1853              if (oldValue+theta*alpha<0.0) {
1854                bestSequence=iSequence;
1855                theta = CoinMax(0.0,oldValue/(-alpha));
1856              }
1857            } else {
1858              // variable going towards upper bound
1859              double bound = upper_[iSequence];
1860              oldValue = bound-oldValue;
1861              if (oldValue-theta*alpha<0.0) {
1862                bestSequence=iSequence;
1863                theta = CoinMax(0.0,oldValue/alpha);
1864              }
1865            }
1866          }
1867        } else {
1868          if (fabs(alpha)>=acceptableBasic) {
1869            if (alpha<0.0) {
1870              // variable going towards lower bound
1871              double bound = lower_[iSequence];
1872              oldValue -= bound;
1873              if (oldValue+basicTheta*alpha<-basicTolerance) {
1874                bestBasicSequence=iSequence;
1875                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/(-alpha));
1876              }
1877            } else {
1878              // variable going towards upper bound
1879              double bound = upper_[iSequence];
1880              oldValue = bound-oldValue;
1881              if (oldValue-basicTheta*alpha<-basicTolerance) {
1882                bestBasicSequence=iSequence;
1883                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/alpha);
1884              }
1885            }
1886          }
1887        }
1888      }
1889      theta_ = CoinMin(theta,basicTheta);
1890      // Now find minimum of function
1891      double objTheta2 =objective_->stepLength(this,solution_,dArray,CoinMin(theta,basicTheta),
1892                                               currentObj,predictedObj,thetaObj);
1893#ifdef CLP_DEBUG
1894      if (handler_->logLevel()&32) 
1895        printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
1896#endif
1897#if MINTYPE==1
1898      if (conjugate) {
1899        double offset;
1900        const double * gradient = objective_->gradient(this,
1901                                                       dArray, offset,
1902                                                       true,0);
1903        double product=0.0;
1904        for (iSequence = 0;iSequence<numberColumns_;iSequence++) {
1905          double alpha = dArray[iSequence];
1906          double value = alpha*gradient[iSequence];
1907          product += value;
1908        }
1909        //#define INCLUDESLACK
1910#ifdef INCLUDESLACK
1911        for (;iSequence<numberColumns_+numberRows_;iSequence++) {
1912          double alpha = dArray[iSequence];
1913          double value = alpha*cost_[iSequence];
1914          product += value;
1915        }
1916#endif
1917        if (product>0.0)
1918          objTheta = djNorm/product;
1919        else
1920          objTheta=COIN_DBL_MAX;
1921#ifndef NDEBUG
1922        if (product<-1.0e-8&&handler_->logLevel()>1) 
1923          printf("bad product %g\n",product);
1924#endif
1925        product = CoinMax(product,0.0);
1926      } else {
1927        objTheta =objTheta2;
1928      }
1929#else
1930      objTheta =objTheta2;
1931#endif
1932      // if very small difference then take pivot (depends on djNorm?)
1933      // distinguish between basic and non-basic
1934      bool chooseObjTheta=objTheta<theta_;
1935      if (chooseObjTheta) {
1936        if (thetaObj<currentObj-1.0e-12&&objTheta+1.0e-10>theta_)
1937          chooseObjTheta=false;
1938        //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
1939        //chooseObjTheta=false;
1940      }
1941      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
1942      if (chooseObjTheta) {
1943        theta_ = objTheta;
1944      } else {
1945        objTheta = CoinMax(objTheta,1.00000001*theta_+1.0e-12);
1946        //if (theta+1.0e-13>basicTheta) {
1947        //theta = CoinMax(theta,1.00000001*basicTheta);
1948        //theta_ = basicTheta;
1949        //}
1950      }
1951      // always out if one variable in and zero move
1952      if (theta_==basicTheta||(sequenceIn_>=0&&!theta_))
1953        finished=true; // come out
1954#ifdef CLP_DEBUG
1955      if (handler_->logLevel()&32) {
1956        printf("%d pass,",nPasses);
1957        if (sequenceIn_>=0)
1958          printf (" Sin %d,",sequenceIn_);
1959        if (basicTheta==theta_)
1960          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1961        else
1962          printf(" basicTheta %g",basicTheta);
1963        if (theta==theta_)
1964          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1965        else
1966          printf(" non-basicTheta %g",theta);
1967        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1968        printf(" djNorm %g\n",djNorm);
1969      }
1970#endif
1971      if (handler_->logLevel()>3&&objTheta!=theta_) {
1972        printf("%d pass obj %g,",nPasses,currentObj);
1973        if (sequenceIn_>=0)
1974          printf (" Sin %d,",sequenceIn_);
1975        if (basicTheta==theta_)
1976          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1977        else
1978          printf(" basicTheta %g",basicTheta);
1979        if (theta==theta_)
1980          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1981        else
1982          printf(" non-basicTheta %g",theta);
1983        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1984        printf(" djNorm %g\n",djNorm);
1985      }
1986      if (objTheta!=theta_) {
1987        //printf("hit boundary after %d passes\n",nPasses);
1988        nTotalPasses+= nPasses;
1989        nPasses=0; // start again
1990      }
1991      if (localPivotMode<10||lastSequenceIn==bestSequence) {
1992        if (theta_==theta&&theta_<basicTheta&&theta_<1.0e-5)
1993          setFlagged(bestSequence); // out of active set
1994      }
1995      // Update solution
1996      for (iSequence=0;iSequence<numberTotal;iSequence++) {
1997        //for (iIndex=0;iIndex<number;iIndex++) {
1998       
1999        //int iSequence = which[iIndex];
2000        double alpha = dArray[iSequence];
2001        if (alpha) {
2002          double value = solution_[iSequence]+theta_*alpha;
2003          solution_[iSequence]=value;
2004          switch(getStatus(iSequence)) {
2005           
2006          case basic:
2007          case isFixed:
2008          case isFree:
2009          case atUpperBound:
2010          case atLowerBound:
2011            nonLinearCost_->setOne(iSequence,value);
2012            break;
2013          case superBasic:
2014            // To get correct action
2015            setStatus(iSequence,isFixed);
2016            nonLinearCost_->setOne(iSequence,value);
2017            //assert (getStatus(iSequence)!=isFixed);
2018            break;
2019          }
2020        }
2021      }
2022      if (objTheta<SMALLTHETA2&&objTheta==theta_) {
2023        if (sequenceIn_>=0&&basicTheta<1.0e-9) {
2024          // try just one
2025          localPivotMode=0;
2026          sequenceIn_=-1;
2027          nTotalPasses+= nPasses;
2028          nPasses=0;
2029          //finished=true;
2030          //objTheta=0.0;
2031          //theta_=0.0;
2032        }
2033        //if (objTheta<1.0e-12)
2034        //finished=true;
2035        //printf("zero move\n");
2036      }
2037#ifdef CLP_DEBUG
2038      if (handler_->logLevel()&32) {
2039        objective_->stepLength(this,solution_,work,0.0,
2040                               currentObj,predictedObj,thetaObj);
2041        printf("current obj %g after update - simple was %g\n",currentObj,simpleObjective);
2042      }
2043#endif
2044      if (sequenceIn_>=0&&!finished&&objTheta>1.0e-4) {
2045        // we made some progress - back to normal
2046        if (localPivotMode<10) {
2047          localPivotMode=0;
2048          sequenceIn_=-1;
2049          nTotalPasses+= nPasses;
2050          nPasses=0;
2051        }
2052#ifdef CLP_DEBUG
2053        if (handler_->logLevel()&32) 
2054          printf("some progress?\n");
2055#endif
2056      }
2057#if CLP_DEBUG > 1
2058      longArray->checkClean();
2059#endif
2060    }
2061#ifdef CLP_DEBUG
2062    if (handler_->logLevel()&32) 
2063      printf("out of loop after %d passes\n",nPasses);
2064#endif
2065    bool ordinaryDj=false;
2066    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
2067    //printf("could try ordinary iteration %g\n",theta_);
2068    if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-15) {
2069      //printf("trying ordinary iteration\n");
2070      ordinaryDj=true;
2071    }
2072    if (!basicTheta&&!ordinaryDj) {
2073      //returnCode=2;
2074      //objTheta=-1.0; // so we fall through
2075    }
2076    assert (theta_<1.0e30); // for now
2077    // See if we need to pivot
2078    if (theta_==basicTheta||ordinaryDj) {
2079      if (!ordinaryDj) {
2080        // find an incoming column which will force pivot
2081        int iRow;
2082        pivotRow_=-1;
2083        for (iRow=0;iRow<numberRows_;iRow++) {
2084          if (pivotVariable_[iRow]==bestBasicSequence) {
2085            pivotRow_=iRow;
2086            break;
2087          }
2088        }
2089        assert (pivotRow_>=0);
2090        // Get good size for pivot
2091        double acceptablePivot=1.0e-7;
2092        if (factorization_->pivots()>10)
2093          acceptablePivot=1.0e-5; // if we have iterated be more strict
2094        else if (factorization_->pivots()>5)
2095          acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
2096        // should be dArray but seems better this way!
2097        double direction = work[bestBasicSequence]>0.0 ? -1.0 : 1.0;
2098        // create as packed
2099        rowArray->createPacked(1,&pivotRow_,&direction);
2100        factorization_->updateColumnTranspose(spare,rowArray);
2101        // put row of tableau in rowArray and columnArray
2102        matrix_->transposeTimes(this,-1.0,
2103                                rowArray,spare,columnArray);
2104        // choose one futhest away from bound which has reasonable pivot
2105        // If increasing we want negative alpha
2106       
2107        double * work2;
2108        int iSection;
2109       
2110        sequenceIn_=-1;
2111        double bestValue=-1.0;
2112        double bestDirection=0.0;
2113        // First pass we take correct direction and large pivots
2114        // then correct direction
2115        // then any
2116        double check[]={1.0e-8,-1.0e-12,-1.0e30};
2117        double mult[]={100.0,1.0,1.0};
2118        for (int iPass=0;iPass<3;iPass++) {
2119          //if (!bestValue&&iPass==2)
2120          //bestValue=-1.0;
2121          double acceptable=acceptablePivot*mult[iPass];
2122          double checkValue=check[iPass];
2123          for (iSection=0;iSection<2;iSection++) {
2124           
2125            int addSequence;
2126           
2127            if (!iSection) {
2128              work2 = rowArray->denseVector();
2129              number = rowArray->getNumElements();
2130              which = rowArray->getIndices();
2131              addSequence = numberColumns_;
2132            } else {
2133              work2 = columnArray->denseVector();
2134              number = columnArray->getNumElements();
2135              which = columnArray->getIndices();
2136              addSequence = 0;
2137            }
2138            int i;
2139           
2140            for (i=0;i<number;i++) {
2141              int iSequence = which[i]+addSequence;
2142              if (flagged(iSequence))
2143                continue;
2144              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
2145              //                  upper_[iSequence]-solution_[iSequence]);
2146              double alpha=work2[i];
2147              // should be dArray but seems better this way!
2148              double change=work[iSequence];
2149              Status thisStatus = getStatus(iSequence);
2150              double direction=0;;
2151              switch(thisStatus) {
2152               
2153              case basic:
2154              case ClpSimplex::isFixed:
2155                break;
2156              case isFree:
2157              case superBasic:
2158                if (alpha<-acceptable&&change>checkValue)
2159                  direction=1.0;
2160                else if (alpha>acceptable&&change<-checkValue)
2161                  direction=-1.0;
2162                break;
2163              case atUpperBound:
2164                if (alpha>acceptable&&change<-checkValue)
2165                  direction=-1.0;
2166                else if (iPass==2&&alpha<-acceptable&&change<-checkValue)
2167                  direction=1.0;
2168                break;
2169              case atLowerBound:
2170                if (alpha<-acceptable&&change>checkValue)
2171                  direction=1.0;
2172                else if (iPass==2&&alpha>acceptable&&change>checkValue)
2173                  direction=-1.0;
2174                break;
2175              }
2176              if (direction) {
2177                if (sequenceIn_!=lastSequenceIn||localPivotMode<10) { 
2178                  if (CoinMin(solution_[iSequence]-lower_[iSequence],
2179                          upper_[iSequence]-solution_[iSequence])>bestValue) {
2180                    bestValue=CoinMin(solution_[iSequence]-lower_[iSequence],
2181                                  upper_[iSequence]-solution_[iSequence]);
2182                    sequenceIn_=iSequence;
2183                    bestDirection=direction;
2184                  }
2185                } else {
2186                  // choose
2187                  bestValue=COIN_DBL_MAX;
2188                  sequenceIn_=iSequence;
2189                  bestDirection=direction;
2190                }
2191              }
2192            }
2193          }
2194          if (sequenceIn_>=0&&bestValue>0.0)
2195            break;
2196        }
2197        if (sequenceIn_>=0) {
2198          valueIn_=solution_[sequenceIn_];
2199          dualIn_=dj_[sequenceIn_];
2200          if (bestDirection<0.0) {
2201            // we want positive dj
2202            if (dualIn_<=0.0) {
2203              //printf("bad dj - xx %g\n",dualIn_);
2204              dualIn_=1.0;
2205            }
2206          } else {
2207            // we want negative dj
2208            if (dualIn_>=0.0) {
2209              //printf("bad dj - xx %g\n",dualIn_);
2210              dualIn_=-1.0;
2211            }
2212          }
2213          lowerIn_=lower_[sequenceIn_];
2214          upperIn_=upper_[sequenceIn_];
2215          if (dualIn_>0.0)
2216            directionIn_ = -1;
2217          else 
2218            directionIn_ = 1;
2219        } else {
2220          //ordinaryDj=true;
2221#ifdef CLP_DEBUG
2222          if (handler_->logLevel()&32) {
2223            printf("no easy pivot - norm %g mode %d\n",djNorm,localPivotMode);
2224            if (rowArray->getNumElements()+columnArray->getNumElements()<12) {
2225              for (iSection=0;iSection<2;iSection++) {
2226               
2227                int addSequence;
2228               
2229                if (!iSection) {
2230                  work2 = rowArray->denseVector();
2231                  number = rowArray->getNumElements();
2232                  which = rowArray->getIndices();
2233                  addSequence = numberColumns_;
2234                } else {
2235                  work2 = columnArray->denseVector();
2236                  number = columnArray->getNumElements();
2237                  which = columnArray->getIndices();
2238                  addSequence = 0;
2239                }
2240                int i;
2241                char section[]={'R','C'};
2242                for (i=0;i<number;i++) {
2243                  int iSequence = which[i]+addSequence;
2244                  if (flagged(iSequence)) {
2245                    printf("%c%d flagged\n",section[iSection],which[i]);
2246                    continue;
2247                  } else {
2248                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
2249                           section[iSection],which[i],status_[iSequence],
2250                           lower_[iSequence],solution_[iSequence],upper_[iSequence],
2251                           work2[i],work[iSequence]);
2252                  }
2253                }
2254              }
2255            }
2256          }
2257#endif
2258          assert (sequenceIn_<0);
2259          justOne=true;
2260          allFinished=false; // Round again
2261          finished=false;
2262          nTotalPasses+= nPasses;
2263          nPasses=0;
2264          if (djNorm<0.9*djNorm0&&djNorm<1.0e-3&&!localPivotMode) {
2265#ifdef CLP_DEBUG
2266            if (handler_->logLevel()&32) 
2267              printf("no pivot - mode %d norms %g %g - unflagging\n",
2268                     localPivotMode,djNorm0,djNorm);
2269#endif
2270            unflag(); //unflagging
2271            returnCode=1;
2272          } else {
2273            returnCode=2; // do single incoming
2274            //returnCode=1;
2275          }
2276        }
2277      }
2278      rowArray->clear();
2279      columnArray->clear();
2280      longArray->clear();
2281      if (ordinaryDj) {
2282        valueIn_=solution_[sequenceIn_];
2283        dualIn_=dj_[sequenceIn_];
2284        lowerIn_=lower_[sequenceIn_];
2285        upperIn_=upper_[sequenceIn_];
2286        if (dualIn_>0.0)
2287          directionIn_ = -1;
2288        else 
2289          directionIn_ = 1;
2290      }
2291      if (returnCode==1)
2292        returnCode=0;
2293    } else {
2294      // round again
2295      longArray->clear();
2296      if (djNorm<1.0e-3&&!localPivotMode) {
2297        if (djNorm==1.2345e-20&&djNorm0>1.0e-4) {
2298#ifdef CLP_DEBUG
2299          if (handler_->logLevel()&32) 
2300            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n",djNorm0,nPasses,
2301                   localPivotMode,returnCode);
2302#endif
2303          //if (!localPivotMode)
2304          //returnCode=2; // force singleton
2305        } else {
2306#ifdef CLP_DEBUG
2307          if (handler_->logLevel()&32) 
2308            printf("unflagging as djNorm %g %g, %d passes\n",djNorm,djNorm0,nPasses);
2309#endif
2310          if (pivotMode>=10) {
2311            pivotMode--;
2312#ifdef CLP_DEBUG
2313            if (handler_->logLevel()&32) 
2314              printf("pivot mode now %d\n",pivotMode);
2315#endif
2316            if (pivotMode==9) 
2317              pivotMode=0;      // switch off fast attempt
2318          }
2319          bool unflagged = unflag()!=0;
2320          if (!unflagged&&djNorm<1.0e-10) {
2321            // ? declare victory
2322            sequenceIn_=-1;
2323            returnCode=0;
2324          }
2325        }
2326      }
2327    }
2328  }
2329  if (djNorm0<1.0e-12*normFlagged) {
2330#ifdef CLP_DEBUG
2331    if (handler_->logLevel()&32) 
2332      printf("unflagging as djNorm %g %g and flagged norm %g\n",djNorm,djNorm0,normFlagged);
2333#endif
2334    unflag();
2335  }
2336  if (saveObj-currentObj<1.0e-5&&nTotalPasses>2000) {
2337    normUnflagged=0.0;
2338    double dualTolerance3 = CoinMin(1.0e-2,1.0e3*dualTolerance_);
2339    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
2340      switch(getStatus(iSequence)) {
2341       
2342      case basic:
2343      case ClpSimplex::isFixed:
2344        break;
2345      case atUpperBound:
2346        if (dj_[iSequence]>dualTolerance3) 
2347          normFlagged += dj_[iSequence]*dj_[iSequence];
2348        break;
2349      case atLowerBound:
2350        if (dj_[iSequence]<-dualTolerance3) 
2351          normFlagged += dj_[iSequence]*dj_[iSequence];
2352        break;
2353      case isFree:
2354      case superBasic:
2355        if (fabs(dj_[iSequence])>dualTolerance3) 
2356          normFlagged += dj_[iSequence]*dj_[iSequence];
2357        break;
2358      }
2359    }
2360    if (handler_->logLevel()>2)
2361      printf("possible optimal  %d %d %g %g\n",
2362             nBigPasses,nTotalPasses,saveObj-currentObj,normFlagged);
2363    if (normFlagged<1.0e-5) {
2364      unflag();
2365      primalColumnPivot_->setLooksOptimal(true);
2366      returnCode=0;
2367    }
2368  }
2369  return returnCode;
2370}
2371/* Do last half of an iteration.
2372   Return codes
2373   Reasons to come out normal mode
2374   -1 normal
2375   -2 factorize now - good iteration
2376   -3 slight inaccuracy - refactorize - iteration done
2377   -4 inaccuracy - refactorize - no iteration
2378   -5 something flagged - go round again
2379   +2 looks unbounded
2380   +3 max iterations (iteration done)
2381   
2382*/
2383int 
2384ClpSimplexNonlinear::pivotNonlinearResult()
2385{
2386
2387  int returnCode=-1;
2388
2389  rowArray_[1]->clear();
2390 
2391  // we found a pivot column
2392  // update the incoming column
2393  unpackPacked(rowArray_[1]);
2394  factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2395  theta_=0.0;
2396  double * work=rowArray_[1]->denseVector();
2397  int number=rowArray_[1]->getNumElements();
2398  int * which=rowArray_[1]->getIndices();
2399  bool keepValue=false;
2400  double saveValue=0.0;
2401  if (pivotRow_>=0) {
2402    sequenceOut_=pivotVariable_[pivotRow_];;
2403    valueOut_ = solution(sequenceOut_);
2404    keepValue=true;
2405    saveValue=valueOut_;
2406    lowerOut_=lower_[sequenceOut_];
2407    upperOut_=upper_[sequenceOut_];
2408    int iIndex;
2409    for (iIndex=0;iIndex<number;iIndex++) {
2410     
2411      int iRow = which[iIndex];
2412      if (iRow==pivotRow_) {
2413        alpha_ = work[iIndex];
2414        break;
2415      }
2416    }
2417  } else {
2418    int iIndex;
2419    double smallest=COIN_DBL_MAX;
2420    for (iIndex=0;iIndex<number;iIndex++) {
2421      int iRow = which[iIndex];
2422      double alpha = work[iIndex];
2423      if (fabs(alpha)>1.0e-6) {
2424        int iPivot = pivotVariable_[iRow];
2425        double distance = CoinMin(upper_[iPivot]-solution_[iPivot],
2426                              solution_[iPivot]-lower_[iPivot]);
2427        if (distance<smallest) {
2428          pivotRow_=iRow;
2429          alpha_=alpha;
2430          smallest=distance;
2431        }
2432      }
2433    }
2434    if (smallest>primalTolerance_) {
2435      smallest=COIN_DBL_MAX;
2436      for (iIndex=0;iIndex<number;iIndex++) {
2437        int iRow = which[iIndex];
2438        double alpha = work[iIndex];
2439        if (fabs(alpha)>1.0e-6) {
2440          double distance = CoinDrand48();
2441          if (distance<smallest) {
2442            pivotRow_=iRow;
2443            alpha_=alpha;
2444            smallest=distance;
2445          }
2446        }
2447      }
2448    }
2449    assert (pivotRow_>=0);
2450    sequenceOut_=pivotVariable_[pivotRow_];;
2451    valueOut_ = solution(sequenceOut_);
2452    lowerOut_=lower_[sequenceOut_];
2453    upperOut_=upper_[sequenceOut_];
2454  }
2455  double newValue = valueOut_ - theta_*alpha_;
2456  bool isSuperBasic=false;
2457  if (valueOut_>=upperOut_-primalTolerance_) {
2458    directionOut_=-1;      // to upper bound
2459    upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2460    upperOut_ = newValue;
2461  } else if (valueOut_<=lowerOut_+primalTolerance_) {
2462    directionOut_=1;      // to lower bound
2463    lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2464  } else {
2465    lowerOut_=valueOut_;
2466    upperOut_=valueOut_;
2467    isSuperBasic=true;
2468    //printf("XX superbasic out\n");
2469  }
2470  dualOut_ = dj_[sequenceOut_];
2471  double checkValue=1.0e-2;
2472  if (largestDualError_>1.0e-5)
2473    checkValue=1.0e-1;
2474  // if stable replace in basis
2475 
2476  int updateStatus = factorization_->replaceColumn(this,
2477                                                   rowArray_[2],
2478                                                   rowArray_[1],
2479                                                   pivotRow_,
2480                                                   alpha_);
2481 
2482  // if no pivots, bad update but reasonable alpha - take and invert
2483  if (updateStatus==2&&
2484      lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2485    updateStatus=4;
2486  if (updateStatus==1||updateStatus==4) {
2487    // slight error
2488    if (factorization_->pivots()>5||updateStatus==4) {
2489      returnCode=-3;
2490    }
2491  } else if (updateStatus==2) {
2492    // major error
2493    // better to have small tolerance even if slower
2494    factorization_->zeroTolerance(1.0e-15);
2495    int maxFactor = factorization_->maximumPivots();
2496    if (maxFactor>10) {
2497      if (forceFactorization_<0)
2498        forceFactorization_= maxFactor;
2499      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2500    } 
2501    // later we may need to unwind more e.g. fake bounds
2502    if(lastGoodIteration_ != numberIterations_) {
2503      clearAll();
2504      pivotRow_=-1;
2505      returnCode=-4;
2506    } else {
2507      // need to reject something
2508      char x = isColumn(sequenceIn_) ? 'C' :'R';
2509      handler_->message(CLP_SIMPLEX_FLAG,messages_)
2510        <<x<<sequenceWithin(sequenceIn_)
2511        <<CoinMessageEol;
2512      setFlagged(sequenceIn_);
2513      progress_->clearBadTimes();
2514      lastBadIteration_ = numberIterations_; // say be more cautious
2515      clearAll();
2516      pivotRow_=-1;
2517      sequenceOut_=-1;
2518      returnCode = -5;
2519     
2520    }
2521    return returnCode;
2522  } else if (updateStatus==3) {
2523    // out of memory
2524    // increase space if not many iterations
2525    if (factorization_->pivots()<
2526        0.5*factorization_->maximumPivots()&&
2527        factorization_->pivots()<200)
2528      factorization_->areaFactor(
2529                                 factorization_->areaFactor() * 1.1);
2530    returnCode =-2; // factorize now
2531  } else if (updateStatus==5) {
2532    problemStatus_=-2; // factorize now
2533  }
2534 
2535  // update primal solution
2536 
2537  double objectiveChange=0.0;
2538  // after this rowArray_[1] is not empty - used to update djs
2539  // If pivot row >= numberRows then may be gub
2540  updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1);
2541 
2542  double oldValue = valueIn_;
2543  if (directionIn_==-1) {
2544    // as if from upper bound
2545    if (sequenceIn_!=sequenceOut_) {
2546      // variable becoming basic
2547      valueIn_ -= fabs(theta_);
2548    } else {
2549      valueIn_=lowerIn_;
2550    }
2551  } else {
2552    // as if from lower bound
2553    if (sequenceIn_!=sequenceOut_) {
2554      // variable becoming basic
2555      valueIn_ += fabs(theta_);
2556    } else {
2557      valueIn_=upperIn_;
2558    }
2559  }
2560  objectiveChange += dualIn_*(valueIn_-oldValue);
2561  // outgoing
2562  if (sequenceIn_!=sequenceOut_) {
2563    if (directionOut_>0) {
2564      valueOut_ = lowerOut_;
2565    } else {
2566      valueOut_ = upperOut_;
2567    }
2568    if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2569      valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2570    else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2571      valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2572    // may not be exactly at bound and bounds may have changed
2573    // Make sure outgoing looks feasible
2574    if (!isSuperBasic)
2575      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2576    solution_[sequenceOut_]=valueOut_;
2577  }
2578  // change cost and bounds on incoming if primal
2579  nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2580  int whatNext=housekeeping(objectiveChange);
2581  if (keepValue)
2582    solution_[sequenceOut_]=saveValue;
2583  if (isSuperBasic)
2584    setStatus(sequenceOut_,superBasic);
2585  //#define CLP_DEBUG
2586#if CLP_DEBUG > 1
2587  {
2588    int ninf= matrix_->checkFeasible(this);
2589    if (ninf)
2590      printf("infeas %d\n",ninf);
2591  }
2592#endif
2593  if (whatNext==1) {
2594    returnCode =-2; // refactorize
2595  } else if (whatNext==2) {
2596    // maximum iterations or equivalent
2597    returnCode=3;
2598  } else if(numberIterations_ == lastGoodIteration_
2599            + 2 * factorization_->maximumPivots()) {
2600    // done a lot of flips - be safe
2601    returnCode =-2; // refactorize
2602  }
2603  // Check event
2604  {
2605    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2606    if (status>=0) {
2607      problemStatus_=5;
2608      secondaryStatus_=ClpEventHandler::endOfIteration;
2609      returnCode=4;
2610    }
2611  }
2612  return returnCode;
2613}
2614// A sequential LP method
2615int 
2616ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
2617{
2618  // Are we minimizing or maximizing
2619  double whichWay=optimizationDirection();
2620  if (whichWay<0.0)
2621    whichWay=-1.0;
2622  else if (whichWay>0.0)
2623    whichWay=1.0;
2624
2625
2626  int numberColumns = this->numberColumns();
2627  int numberRows = this->numberRows();
2628  double * columnLower = this->columnLower();
2629  double * columnUpper = this->columnUpper();
2630  double * solution = this->primalColumnSolution();
2631 
2632  if (objective_->type()<2) {
2633    // no nonlinear part
2634    return ClpSimplex::primal(0);
2635  }
2636  // Get list of non linear columns
2637  char * markNonlinear = new char[numberColumns];
2638  memset(markNonlinear,0,numberColumns);
2639  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
2640 
2641  if (!numberNonLinearColumns) {
2642    delete [] markNonlinear;
2643    // no nonlinear part
2644    return ClpSimplex::primal(0);
2645  }
2646  int iColumn;
2647  int * listNonLinearColumn = new int[numberNonLinearColumns];
2648  numberNonLinearColumns=0;
2649  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2650    if(markNonlinear[iColumn])
2651      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
2652  }
2653  // Replace objective
2654  ClpObjective * trueObjective = objective_;
2655  objective_=new ClpLinearObjective(NULL,numberColumns);
2656  double * objective = this->objective();
2657 
2658  // get feasible
2659  if (this->status()<0||numberPrimalInfeasibilities())
2660    ClpSimplex::primal(1);
2661  // still infeasible
2662  if (numberPrimalInfeasibilities()) {
2663    delete [] listNonLinearColumn;
2664    return 0;
2665  }
2666  algorithm_=1;
2667  int jNon;
2668  int * last[3];
2669 
2670  double * trust = new double[numberNonLinearColumns];
2671  double * trueLower = new double[numberNonLinearColumns];
2672  double * trueUpper = new double[numberNonLinearColumns];
2673  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2674    iColumn=listNonLinearColumn[jNon];
2675    trust[jNon]=0.5;
2676    trueLower[jNon]=columnLower[iColumn];
2677    trueUpper[jNon]=columnUpper[iColumn];
2678    if (solution[iColumn]<trueLower[jNon])
2679      solution[iColumn]=trueLower[jNon];
2680    else if (solution[iColumn]>trueUpper[jNon])
2681      solution[iColumn]=trueUpper[jNon];
2682  }
2683  int saveLogLevel=logLevel();
2684  int iPass;
2685  double lastObjective=1.0e31;
2686  double * saveSolution = new double [numberColumns];
2687  double * saveRowSolution = new double [numberRows];
2688  memset(saveRowSolution,0,numberRows*sizeof(double));
2689  double * savePi = new double [numberRows];
2690  double * safeSolution = new double [numberColumns];
2691  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
2692#define MULTIPLE 0
2693#if MULTIPLE>2
2694  // Duplication but doesn't really matter
2695  double * saveSolutionM[MULTIPLE};
2696  for (jNon=0;jNon<MULTIPLE;jNon++) {
2697    saveSolutionM[jNon]=new double[numberColumns];
2698    memcpy(saveSolutionM,solution,numberColumns*sizeof(double));
2699  }
2700#endif
2701  double targetDrop=1.0e31;
2702  double objectiveOffset;
2703  getDblParam(ClpObjOffset,objectiveOffset);
2704  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
2705  for (iPass=0;iPass<3;iPass++) {
2706    last[iPass]=new int[numberNonLinearColumns];
2707    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
2708      last[iPass][jNon]=0;
2709  }
2710  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2711  int goodMove=-2;
2712  char * statusCheck = new char[numberColumns];
2713  double * changeRegion = new double [numberColumns];
2714  double offset=0.0;
2715  double objValue=0.0;
2716  int exitPass = 2*numberPasses+10;
2717  for (iPass=0;iPass<numberPasses;iPass++) {
2718    exitPass--;
2719    // redo objective
2720    offset=0.0;
2721    objValue=-objectiveOffset;
2722    // make sure x updated
2723    trueObjective->newXValues();
2724    double theta=-1.0;
2725    double maxTheta=COIN_DBL_MAX;
2726    //maxTheta=1.0;
2727    if (iPass) {
2728      int jNon=0;
2729      for (iColumn=0;iColumn<numberColumns;iColumn++) { 
2730        changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
2731        double alpha = changeRegion[iColumn];
2732        double oldValue = saveSolution[iColumn];
2733        if (markNonlinear[iColumn]==0) {
2734          // linear
2735          if (alpha<-1.0e-15) {
2736            // variable going towards lower bound
2737            double bound = columnLower[iColumn];
2738            oldValue -= bound;
2739            if (oldValue+maxTheta*alpha<0.0) {
2740              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2741            }
2742          } else if (alpha>1.0e-15) {
2743            // variable going towards upper bound
2744            double bound = columnUpper[iColumn];
2745            oldValue = bound-oldValue;
2746            if (oldValue-maxTheta*alpha<0.0) {
2747              maxTheta = CoinMax(0.0,oldValue/alpha);
2748            }
2749          }
2750        } else {
2751          // nonlinear
2752          if (alpha<-1.0e-15) {
2753            // variable going towards lower bound
2754            double bound = trueLower[jNon];
2755            oldValue -= bound;
2756            if (oldValue+maxTheta*alpha<0.0) {
2757              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2758            }
2759          } else if (alpha>1.0e-15) {
2760            // variable going towards upper bound
2761            double bound = trueUpper[jNon];
2762            oldValue = bound-oldValue;
2763            if (oldValue-maxTheta*alpha<0.0) {
2764              maxTheta = CoinMax(0.0,oldValue/alpha);
2765            }
2766          }
2767          jNon++;
2768        }
2769      }
2770      // make sure both accurate
2771      memset(rowActivity_,0,numberRows_*sizeof(double));
2772      times(1.0,solution,rowActivity_);
2773      memset(saveRowSolution,0,numberRows_*sizeof(double));
2774      times(1.0,saveSolution,saveRowSolution);
2775      for (int iRow=0;iRow<numberRows;iRow++) { 
2776        double alpha =rowActivity_[iRow]-saveRowSolution[iRow];
2777        double oldValue = saveRowSolution[iRow];
2778        if (alpha<-1.0e-15) {
2779          // variable going towards lower bound
2780          double bound = rowLower_[iRow];
2781          oldValue -= bound;
2782          if (oldValue+maxTheta*alpha<0.0) {
2783            maxTheta = CoinMax(0.0,oldValue/(-alpha));
2784          }
2785        } else if (alpha>1.0e-15) {
2786          // variable going towards upper bound
2787          double bound = rowUpper_[iRow];
2788          oldValue = bound-oldValue;
2789          if (oldValue-maxTheta*alpha<0.0) {
2790            maxTheta = CoinMax(0.0,oldValue/alpha);
2791          }
2792        }
2793      }
2794    } else {
2795      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2796        changeRegion[iColumn]=0.0;
2797        saveSolution[iColumn]=solution[iColumn];
2798      }
2799      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
2800    }
2801    // get current value anyway
2802    double predictedObj,thetaObj;
2803    double maxTheta2 = 2.0; // to work out a b c
2804    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
2805                                             objValue,predictedObj,thetaObj);
2806    int lastMoveStatus = goodMove;
2807    if (goodMove>=0) {
2808      theta = CoinMin(theta2,maxTheta);
2809#ifdef CLP_DEBUG
2810      if (handler_->logLevel()&32) 
2811        printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
2812               theta,objValue,thetaObj,predictedObj);
2813#endif
2814      if (theta>0.0&&theta<=1.0) {
2815        // update solution
2816        double lambda = 1.0-theta;
2817        for (iColumn=0;iColumn<numberColumns;iColumn++) 
2818          solution[iColumn] = lambda * saveSolution[iColumn] 
2819            + theta * solution[iColumn];
2820        memset(rowActivity_,0,numberRows_*sizeof(double));
2821        times(1.0,solution,rowActivity_);
2822        if (lambda>0.999) {
2823          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
2824          memcpy(status_,saveStatus,numberRows+numberColumns);
2825        }
2826        // Do local minimization
2827#define LOCAL
2828#ifdef LOCAL
2829        bool absolutelyOptimal=false;
2830        int saveScaling = scalingFlag();
2831        scaling(0);
2832        int savePerturbation=perturbation_;
2833        perturbation_=100;
2834        if (saveLogLevel==1)
2835          setLogLevel(0);
2836        int status=startup(1);
2837        if (!status) {
2838          int numberTotal = numberRows_+numberColumns_;
2839          // resize arrays
2840          for (int i=0;i<4;i++) {
2841            rowArray_[i]->reserve(CoinMax(numberRows_+numberColumns_,rowArray_[i]->capacity()));
2842          }
2843          CoinIndexedVector * longArray = rowArray_[3];
2844          CoinIndexedVector * rowArray = rowArray_[0];
2845          //CoinIndexedVector * columnArray = columnArray_[0];
2846          CoinIndexedVector * spare = rowArray_[1];
2847          double * work=longArray->denseVector();
2848          int * which=longArray->getIndices();
2849          int nPass=100;
2850          //bool conjugate=false;
2851          // Put back true bounds
2852          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2853            int iColumn=listNonLinearColumn[jNon];
2854            double value;
2855            value=trueLower[jNon];
2856            trueLower[jNon]=lower_[iColumn];
2857            lower_[iColumn]=value;
2858            value=trueUpper[jNon];
2859            trueUpper[jNon]=upper_[iColumn];
2860            upper_[iColumn]=value;
2861            switch(getStatus(iColumn)) {
2862             
2863            case basic:
2864            case isFree:
2865            case superBasic:
2866              break;
2867            case isFixed:
2868            case atUpperBound:
2869            case atLowerBound:
2870              if (solution_[iColumn]>lower_[iColumn]+primalTolerance_) {
2871                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2872                  setStatus(iColumn,superBasic);
2873                } else {
2874                  setStatus(iColumn,atUpperBound);
2875                }
2876              } else {
2877                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2878                  setStatus(iColumn,atLowerBound);
2879                } else {
2880                  setStatus(iColumn,isFixed);
2881                }
2882              }
2883              break;
2884            }
2885          }
2886          for (int jPass=0;jPass<nPass;jPass++) {
2887            trueObjective->reducedGradient(this,dj_,true);
2888            // get direction vector in longArray
2889            longArray->clear();
2890            double normFlagged=0.0;
2891            double normUnflagged=0.0;
2892            int numberNonBasic=0;
2893            directionVector(longArray,spare,rowArray,0,
2894                            normFlagged,normUnflagged,numberNonBasic);
2895            if (normFlagged+normUnflagged<1.0e-8) {
2896              absolutelyOptimal=true;
2897              break;  //looks optimal
2898            }
2899            double djNorm=0.0;
2900            int iIndex;
2901            for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
2902              int iSequence = which[iIndex];
2903              double alpha = work[iSequence];
2904              djNorm += alpha*alpha;
2905            }
2906            //if (!jPass)
2907            //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
2908            //int number=longArray->getNumElements();
2909            if (!numberNonBasic) {
2910              // looks optimal
2911              absolutelyOptimal=true;
2912              break;
2913            }
2914            int bestSequence=-1;
2915            double theta = 1.0e30;
2916            int iSequence;
2917            for (iSequence=0;iSequence<numberTotal;iSequence++) {
2918              double alpha = work[iSequence];
2919              double oldValue = solution_[iSequence];
2920              if (alpha<-1.0e-15) {
2921                // variable going towards lower bound
2922                double bound = lower_[iSequence];
2923                oldValue -= bound;
2924                if (oldValue+theta*alpha<0.0) {
2925                  bestSequence=iSequence;
2926                  theta = CoinMax(0.0,oldValue/(-alpha));
2927                }
2928              } else if (alpha>1.0e-15) {
2929                // variable going towards upper bound
2930                double bound = upper_[iSequence];
2931                oldValue = bound-oldValue;
2932                if (oldValue-theta*alpha<0.0) {
2933                  bestSequence=iSequence;
2934                  theta = CoinMax(0.0,oldValue/alpha);
2935                }
2936              }
2937            }
2938            // Now find minimum of function
2939            double currentObj;
2940            double predictedObj;
2941            double thetaObj;
2942            // need to use true objective
2943            double * saveCost = cost_;
2944            cost_=NULL;
2945            double objTheta =trueObjective->stepLength(this,solution_,work,theta,
2946                                                       currentObj,predictedObj,thetaObj);
2947            cost_=saveCost;
2948#ifdef CLP_DEBUG
2949            if (handler_->logLevel()&32) 
2950              printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2951#endif
2952            //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2953            //printf("objTheta %g theta %g\n",objTheta,theta);
2954            if (theta>objTheta) { 
2955              theta=objTheta;
2956              thetaObj=predictedObj;
2957            }
2958            // update one used outside
2959            objValue=currentObj;
2960            if (theta>1.0e-9&&
2961                (currentObj-thetaObj<-CoinMax(1.0e-8,1.0e-15*fabs(currentObj))||jPass<5)) {
2962              // Update solution
2963              for (iSequence=0;iSequence<numberTotal;iSequence++) {
2964                double alpha = work[iSequence];
2965                if (alpha) {
2966                  double value = solution_[iSequence]+theta*alpha;
2967                  solution_[iSequence]=value;
2968                  switch(getStatus(iSequence)) {
2969                   
2970                  case basic:
2971                  case isFixed:
2972                  case isFree:
2973                    break;
2974                  case atUpperBound:
2975                  case atLowerBound:
2976                  case superBasic:
2977                    if (value>lower_[iSequence]+primalTolerance_) {
2978                      if(value<upper_[iSequence]-primalTolerance_) {
2979                        setStatus(iSequence,superBasic);
2980                      } else {
2981                        setStatus(iSequence,atUpperBound);
2982                      }
2983                    } else {
2984                      if(value<upper_[iSequence]-primalTolerance_) {
2985                        setStatus(iSequence,atLowerBound);
2986                      } else {
2987                        setStatus(iSequence,isFixed);
2988                      }
2989                    }
2990                    break;
2991                  }
2992                }
2993              }
2994            } else {
2995              break;
2996            }
2997          }
2998          // Put back fake bounds
2999          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3000            int iColumn=listNonLinearColumn[jNon];
3001            double value;
3002            value=trueLower[jNon];
3003            trueLower[jNon]=lower_[iColumn];
3004            lower_[iColumn]=value;
3005            value=trueUpper[jNon];
3006            trueUpper[jNon]=upper_[iColumn];
3007            upper_[iColumn]=value;
3008          }
3009        }
3010        problemStatus_=0;
3011        finish();
3012        scaling(saveScaling);
3013        perturbation_=savePerturbation;
3014        setLogLevel(saveLogLevel);
3015#endif
3016        // redo rowActivity
3017        memset(rowActivity_,0,numberRows_*sizeof(double));
3018        times(1.0,solution,rowActivity_);
3019        if (theta>0.99999&&theta2<1.9&&!absolutelyOptimal) {
3020          // If big changes then tighten
3021          /*  thetaObj is objvalue + a*2*2 +b*2
3022              predictedObj is objvalue + a*theta2*theta2 +b*theta2
3023          */
3024          double rhs1 = thetaObj-objValue;
3025          double rhs2 = predictedObj-objValue;
3026          double subtractB = theta2*0.5;
3027          double a = (rhs2-subtractB*rhs1)/(theta2*theta2-4.0*subtractB);
3028          double b = 0.5*(rhs1 - 4.0*a);
3029          if (fabs(a+b)>1.0e-2) {
3030            // tighten all
3031            goodMove=-1;
3032          }
3033        }
3034      }
3035    }
3036    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3037           numberColumns*sizeof(double));
3038    //printf("offset comp %g orig %g\n",offset,objectiveOffset);
3039    // could do faster
3040    trueObjective->stepLength(this,solution,changeRegion,0.0,
3041                                             objValue,predictedObj,thetaObj);
3042    printf("offset comp %g orig %g - obj from stepLength %g\n",offset,objectiveOffset,objValue);
3043    setDblParam(ClpObjOffset,objectiveOffset+offset);
3044    int * temp=last[2];
3045    last[2]=last[1];
3046    last[1]=last[0];
3047    last[0]=temp;
3048    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3049      iColumn=listNonLinearColumn[jNon];
3050      double change = solution[iColumn]-saveSolution[iColumn];
3051      if (change<-1.0e-5) {
3052        if (fabs(change+trust[jNon])<1.0e-5) 
3053          temp[jNon]=-1;
3054        else
3055          temp[jNon]=-2;
3056      } else if(change>1.0e-5) {
3057        if (fabs(change-trust[jNon])<1.0e-5) 
3058          temp[jNon]=1;
3059        else
3060          temp[jNon]=2;
3061      } else {
3062        temp[jNon]=0;
3063      }
3064    } 
3065    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3066    double maxDelta=0.0;
3067    if (goodMove>=0) {
3068      if (objValue-lastObjective<=1.0e-15*fabs(lastObjective)) 
3069        goodMove=1;
3070      else
3071        goodMove=0;
3072    } else {
3073      maxDelta=1.0e10;
3074    }
3075    double maxGap=0.0;
3076    int numberSmaller=0;
3077    int numberSmaller2=0;
3078    int numberLarger=0;
3079    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3080      iColumn=listNonLinearColumn[jNon];
3081      maxDelta = CoinMax(maxDelta,
3082                     fabs(solution[iColumn]-saveSolution[iColumn]));
3083      if (goodMove>0) {
3084        if (last[0][jNon]*last[1][jNon]<0) {
3085          // halve
3086          trust[jNon] *= 0.5;
3087          numberSmaller2++;
3088        } else {
3089          if (last[0][jNon]==last[1][jNon]&&
3090              last[0][jNon]==last[2][jNon])
3091            trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6); 
3092          numberLarger++;
3093        }
3094      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
3095        trust[jNon] *= 0.2;
3096        numberSmaller++;
3097      }
3098      maxGap = CoinMax(maxGap,trust[jNon]);
3099    }
3100#ifdef CLP_DEBUG
3101          if (handler_->logLevel()&32) 
3102            std::cout<<"largest gap is "<<maxGap<<" "
3103                     <<numberSmaller+numberSmaller2<<" reduced ("
3104                     <<numberSmaller<<" badMove ), "
3105                     <<numberLarger<<" increased"<<std::endl;
3106#endif
3107    if (iPass>10000) {
3108      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3109        trust[jNon] *=0.0001;
3110    }
3111    if(lastMoveStatus==-1&&goodMove==-1)
3112      goodMove=1; // to force solve
3113    if (goodMove>0) {
3114      double drop = lastObjective-objValue;
3115      handler_->message(CLP_SLP_ITER,messages_)
3116        <<iPass<<objValue-objectiveOffset
3117        <<drop<<maxDelta
3118        <<CoinMessageEol;
3119      if (iPass>20&&drop<1.0e-12*fabs(objValue))
3120        drop=0.999e-4; // so will exit
3121      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999) {
3122        if (handler_->logLevel()>1) 
3123          std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
3124        break;
3125      }
3126    } else if (!numberSmaller&&iPass>1) {
3127      if (handler_->logLevel()>1) 
3128          std::cout<<"Exiting as all gaps small"<<std::endl;
3129        break;
3130    }
3131    if (!iPass)
3132      goodMove=1;
3133    targetDrop=0.0;
3134    double * r = this->dualColumnSolution();
3135    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3136      iColumn=listNonLinearColumn[jNon];
3137      columnLower[iColumn]=CoinMax(solution[iColumn]
3138                               -trust[jNon],
3139                               trueLower[jNon]);
3140      columnUpper[iColumn]=CoinMin(solution[iColumn]
3141                               +trust[jNon],
3142                               trueUpper[jNon]);
3143    }
3144    if (iPass) {
3145      // get reduced costs
3146      this->matrix()->transposeTimes(savePi,
3147                                     this->dualColumnSolution());
3148      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3149        iColumn=listNonLinearColumn[jNon];
3150        double dj = objective[iColumn]-r[iColumn];
3151        r[iColumn]=dj;
3152        if (dj<-dualTolerance_) 
3153          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
3154        else if (dj>dualTolerance_)
3155          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
3156      }
3157    } else {
3158      memset(r,0,numberColumns*sizeof(double));
3159    }
3160#if 0
3161    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3162      iColumn=listNonLinearColumn[jNon];
3163      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
3164        columnLower[iColumn]=CoinMax(solution[iColumn],
3165                                 trueLower[jNon]);
3166        columnUpper[iColumn]=CoinMin(solution[iColumn]
3167                                 +trust[jNon],
3168                                 trueUpper[jNon]);
3169      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
3170        columnLower[iColumn]=CoinMax(solution[iColumn]
3171                                 -trust[jNon],
3172                                 trueLower[jNon]);
3173        columnUpper[iColumn]=CoinMin(solution[iColumn],
3174                                 trueUpper[jNon]);
3175      } else {
3176        columnLower[iColumn]=CoinMax(solution[iColumn]
3177                                 -trust[jNon],
3178                                 trueLower[jNon]);
3179        columnUpper[iColumn]=CoinMin(solution[iColumn]
3180                                 +trust[jNon],
3181                                 trueUpper[jNon]);
3182      }
3183    }
3184#endif
3185    if (goodMove>0) {
3186      memcpy(saveSolution,solution,numberColumns*sizeof(double));
3187      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
3188      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
3189      memcpy(saveStatus,status_,numberRows+numberColumns);
3190#if MULTIPLE>2
3191      double * tempSol=saveSolutionM[0];
3192      for (jNon=0;jNon<MULTIPLE-1;jNon++) {
3193        saveSolutionM[jNon]=saveSolutionM[jNon+1];
3194      }
3195      saveSolutionM[MULTIPLE-1]=tempSol;
3196      memcpy(tempSol,solution,numberColumns*sizeof(double));
3197     
3198#endif
3199     
3200#ifdef CLP_DEBUG
3201      if (handler_->logLevel()&32) 
3202        std::cout<<"Pass - "<<iPass
3203                 <<", target drop is "<<targetDrop
3204                 <<std::endl;
3205#endif
3206      lastObjective = objValue;
3207      if (targetDrop<CoinMax(1.0e-8,CoinMin(1.0e-6,1.0e-6*fabs(objValue)))&&goodMove&&iPass>3) {
3208        if (handler_->logLevel()>1) 
3209        printf("Exiting on target drop %g\n",targetDrop);
3210        break;
3211      }
3212#ifdef CLP_DEBUG
3213      {
3214        double * r = this->dualColumnSolution();
3215        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3216          iColumn=listNonLinearColumn[jNon];
3217          if (handler_->logLevel()&32) 
3218            printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
3219                   jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
3220                   r[iColumn],statusCheck[iColumn],columnLower[iColumn],
3221                   columnUpper[iColumn]);
3222        }
3223      }
3224#endif
3225      //setLogLevel(63);
3226      this->scaling(false);
3227      if (saveLogLevel==1)
3228        setLogLevel(0);
3229      ClpSimplex::primal(1);
3230      algorithm_=1;
3231      setLogLevel(saveLogLevel);
3232#ifdef CLP_DEBUG
3233      if (this->status()) {
3234        writeMps("xx.mps");
3235      }
3236#endif
3237      if (this->status()==1) {
3238        // not feasible ! - backtrack and exit
3239        // use safe solution
3240        memcpy(solution,safeSolution,numberColumns*sizeof(double));
3241        memcpy(saveSolution,solution,numberColumns*sizeof(double));
3242        memset(rowActivity_,0,numberRows_*sizeof(double));
3243        times(1.0,solution,rowActivity_);
3244        memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
3245        memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
3246        memcpy(status_,saveStatus,numberRows+numberColumns);
3247        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3248          iColumn=listNonLinearColumn[jNon];
3249          columnLower[iColumn]=CoinMax(solution[iColumn]
3250                                       -trust[jNon],
3251                                       trueLower[jNon]);
3252          columnUpper[iColumn]=CoinMin(solution[iColumn]
3253                                       +trust[jNon],
3254                                       trueUpper[jNon]);
3255        }
3256        break;
3257      } else {
3258        // save in case problems
3259        memcpy(safeSolution,solution,numberColumns*sizeof(double));
3260      }
3261      if (problemStatus_==3)
3262        break; // probably user interrupt
3263      goodMove=1;
3264    } else {
3265      // bad pass - restore solution
3266#ifdef CLP_DEBUG
3267      if (handler_->logLevel()&32) 
3268        printf("Backtracking\n");
3269#endif
3270      memcpy(solution,saveSolution,numberColumns*sizeof(double));
3271      memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
3272      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
3273      memcpy(status_,saveStatus,numberRows+numberColumns);
3274      iPass--;
3275      assert (exitPass>0);
3276      goodMove=-1;
3277    }
3278  }
3279#if MULTIPLE>2
3280  for (jNon=0;jNon<MULTIPLE;jNon++) {
3281    delete [] saveSolutionM[jNon];
3282  }
3283#endif
3284  // restore solution
3285  memcpy(solution,saveSolution,numberColumns*sizeof(double));
3286  memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
3287  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3288    iColumn=listNonLinearColumn[jNon];
3289    columnLower[iColumn]=CoinMax(solution[iColumn],
3290                             trueLower[jNon]);
3291    columnUpper[iColumn]=CoinMin(solution[iColumn],
3292                             trueUpper[jNon]);
3293  }
3294  delete [] markNonlinear;
3295  delete [] statusCheck;
3296  delete [] savePi;
3297  delete [] saveStatus;
3298  // redo objective
3299  memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3300         numberColumns*sizeof(double));
3301  ClpSimplex::primal(1);
3302  delete objective_;
3303  objective_=trueObjective;
3304  // redo values
3305  setDblParam(ClpObjOffset,objectiveOffset);
3306  objectiveValue_ += whichWay*offset;
3307  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3308    iColumn=listNonLinearColumn[jNon];
3309    columnLower[iColumn]= trueLower[jNon];
3310    columnUpper[iColumn]= trueUpper[jNon];
3311  }
3312  delete [] saveSolution;
3313  delete [] safeSolution;
3314  delete [] saveRowSolution;
3315  for (iPass=0;iPass<3;iPass++) 
3316    delete [] last[iPass];
3317  delete [] trust;
3318  delete [] trueUpper;
3319  delete [] trueLower;
3320  delete [] listNonLinearColumn;
3321  delete [] changeRegion;
3322  // temp
3323  //setLogLevel(63);
3324  return 0;
3325}
3326/* Primal algorithm for nonlinear constraints
3327   Using a semi-trust region approach as for pooling problem
3328   This is in because I have it lying around
3329   
3330*/
3331int 
3332ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
3333                               int numberPasses, double deltaTolerance)
3334{
3335  if (!numberConstraints) {
3336    // no nonlinear constraints - may be nonlinear objective
3337    return primalSLP(numberPasses,deltaTolerance);
3338  }
3339  // check all matrix for odd rows is empty
3340  int iConstraint;
3341  int numberBad=0;
3342  CoinPackedMatrix * columnCopy = matrix();
3343  // Get a row copy in standard format
3344  CoinPackedMatrix copy;
3345  copy.reverseOrderedCopyOf(*columnCopy);
3346  // get matrix data pointers
3347  //const int * column = copy.getIndices();
3348  //const CoinBigIndex * rowStart = copy.getVectorStarts();
3349  const int * rowLength = copy.getVectorLengths(); 
3350  //const double * elementByRow = copy.getElements();
3351  int numberArtificials=0;
3352  // see how many extra we need
3353  CoinBigIndex numberExtra=0;
3354  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3355    ClpConstraint * constraint = constraints[iConstraint];
3356    int iRow = constraint->rowNumber();
3357    assert (iRow>=0);
3358    int n = constraint->numberCoefficients() -rowLength[iRow];
3359    numberExtra += n;
3360    if (iRow>=numberRows_)
3361      numberBad++;
3362    else if (rowLength[iRow]&&n)
3363      numberBad++;
3364    if (rowLower_[iRow]>-1.0e20)
3365      numberArtificials++;
3366    if (rowUpper_[iRow]<1.0e20)
3367      numberArtificials++;
3368  }
3369  if (numberBad)
3370    return numberBad;
3371  ClpObjective * trueObjective = NULL;
3372  if (objective_->type()>=2) {
3373    // Replace objective
3374    trueObjective = objective_;
3375    objective_=new ClpLinearObjective(NULL,numberColumns_);
3376  }
3377  ClpSimplex newModel(*this);
3378  // we can put back true objective
3379  if (trueObjective) {
3380    // Replace objective
3381    delete objective_;
3382    objective_=trueObjective;
3383  }
3384  int numberColumns2 = numberColumns_;
3385  double penalty=1.0e9;
3386  if (numberArtificials) {
3387    numberColumns2 += numberArtificials;
3388    int * addStarts = new int [numberArtificials+1];
3389    int * addRow = new int[numberArtificials];
3390    double * addElement = new double[numberArtificials];
3391    addStarts[0]=0;
3392    double * addCost = new double [numberArtificials];
3393    numberArtificials=0;
3394    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3395      ClpConstraint * constraint = constraints[iConstraint];
3396      int iRow = constraint->rowNumber();
3397      if (rowLower_[iRow]>-1.0e20) {
3398        addRow[numberArtificials]=iRow;
3399        addElement[numberArtificials]=1.0;
3400        addCost[numberArtificials]=penalty;
3401        numberArtificials++;
3402        addStarts[numberArtificials]=numberArtificials;
3403      }
3404      if (rowUpper_[iRow]<1.0e20) {
3405        addRow[numberArtificials]=iRow;
3406        addElement[numberArtificials]=-1.0;
3407        addCost[numberArtificials]=penalty;
3408        numberArtificials++;
3409        addStarts[numberArtificials]=numberArtificials;
3410      }
3411    }
3412    newModel.addColumns(numberArtificials,NULL,NULL,addCost,
3413                       addStarts,addRow,addElement);
3414    delete [] addStarts;
3415    delete [] addRow;
3416    delete [] addElement;
3417    delete [] addCost;
3418  }
3419  // find nonlinear columns
3420  int * listNonLinearColumn = new int [numberColumns_];
3421  char * mark = new char[numberColumns_];
3422  memset(mark,0,numberColumns_);
3423  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3424    ClpConstraint * constraint = constraints[iConstraint];
3425    constraint->markNonlinear(mark);
3426  }
3427  if (trueObjective)
3428    trueObjective->markNonlinear(mark);
3429  int iColumn;
3430  int numberNonLinearColumns=0;
3431  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3432    if (mark[iColumn])
3433      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
3434  }
3435  // build row copy of matrix with space for nonzeros
3436  // Get column copy
3437  columnCopy = newModel.matrix();
3438  copy.reverseOrderedCopyOf(*columnCopy);
3439  // get matrix data pointers
3440  const int * column = copy.getIndices();
3441  const CoinBigIndex * rowStart = copy.getVectorStarts();
3442  rowLength = copy.getVectorLengths(); 
3443  const double * elementByRow = copy.getElements();
3444  numberExtra +=copy.getNumElements();
3445  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
3446  int * newColumn = new int[numberExtra];
3447  double * newElement = new double[numberExtra];
3448  newStarts[0]=0;
3449  int * backRow = new int [numberRows_];
3450  int iRow;
3451  for (iRow=0;iRow<numberRows_;iRow++) 
3452    backRow[iRow]=-1;
3453  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3454    ClpConstraint * constraint = constraints[iConstraint];
3455    iRow = constraint->rowNumber();
3456    backRow[iRow]=iConstraint;
3457  }
3458  numberExtra=0;
3459  for (iRow=0;iRow<numberRows_;iRow++) {
3460    if (backRow[iRow]<0) {
3461      // copy normal
3462      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3463           j++) {
3464        newColumn[numberExtra]=column[j];
3465        newElement[numberExtra++] = elementByRow[j];
3466      }
3467    } else {
3468      ClpConstraint * constraint = constraints[backRow[iRow]];
3469      assert(iRow == constraint->rowNumber());
3470      int numberArtificials=0;
3471      if (rowLower_[iRow]>-1.0e20)
3472        numberArtificials++;
3473      if (rowUpper_[iRow]<1.0e20)
3474        numberArtificials++;
3475      if (numberArtificials==rowLength[iRow]) {
3476        // all possible
3477        memset(mark,0,numberColumns_);
3478        constraint->markNonzero(mark);
3479        for (int k=0;k<numberColumns_;k++) {
3480          if (mark[k]) {
3481            newColumn[numberExtra]=k;
3482            newElement[numberExtra++] = 1.0;
3483          }
3484        }
3485        // copy artificials
3486        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3487             j++) {
3488          newColumn[numberExtra]=column[j];
3489          newElement[numberExtra++] = elementByRow[j];
3490        }
3491      } else {
3492        // there already
3493        // copy
3494        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
3495             j++) {
3496          newColumn[numberExtra]=column[j];
3497          assert (elementByRow[j]);
3498          newElement[numberExtra++] = elementByRow[j];
3499        }
3500      }
3501    }
3502    newStarts[iRow+1]=numberExtra;
3503  }
3504  delete [] backRow;
3505  CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_,
3506                              numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0);
3507  delete [] newStarts;
3508  delete [] newColumn;
3509  delete [] newElement;
3510  delete [] mark;
3511  // get feasible
3512  newModel.primal(1);
3513  // still infeasible
3514  if (newModel.numberPrimalInfeasibilities()) {
3515    delete [] listNonLinearColumn;
3516    return 0;
3517  }
3518  int numberRows = newModel.numberRows();
3519  double * columnLower = newModel.columnLower();
3520  double * columnUpper = newModel.columnUpper();
3521  double * objective = newModel.objective();
3522  double * rowLower = newModel.rowLower();
3523  double * rowUpper = newModel.rowUpper(); 
3524  double * solution = newModel.primalColumnSolution();
3525  int jNon;
3526  int * last[3];
3527 
3528  double * trust = new double[numberNonLinearColumns];
3529  double * trueLower = new double[numberNonLinearColumns];
3530  double * trueUpper = new double[numberNonLinearColumns];
3531  double objectiveOffset;
3532  double objectiveOffset2;
3533  getDblParam(ClpObjOffset,objectiveOffset);
3534  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3535    iColumn=listNonLinearColumn[jNon];
3536    double upper = columnUpper[iColumn];
3537    double lower = columnLower[iColumn];
3538    if (solution[iColumn]<trueLower[jNon])
3539      solution[iColumn]=trueLower[jNon];
3540    else if (solution[iColumn]>trueUpper[jNon])
3541      solution[iColumn]=trueUpper[jNon];
3542#if 0
3543    double large = CoinMax(1000.0,10.0*fabs(solution[iColumn]));
3544    if (upper>1.0e10)
3545      upper = solution[iColumn]+large;
3546    if (lower<-1.0e10)
3547      lower = solution[iColumn]-large;
3548#else
3549    upper = solution[iColumn]+0.5;
3550    lower = solution[iColumn]-0.5;
3551#endif
3552    //columnUpper[iColumn]=upper;
3553    trust[jNon]=0.05*(1.0+upper-lower);
3554    trueLower[jNon]=columnLower[iColumn];
3555    //trueUpper[jNon]=upper;
3556    trueUpper[jNon]=columnUpper[iColumn];
3557  }
3558  bool tryFix=false;
3559  int iPass;
3560  double lastObjective=1.0e31;
3561  double * saveSolution = new double [numberColumns2+numberRows];
3562  char * saveStatus = new char [numberColumns2+numberRows];
3563  double targetDrop=1.0e31;
3564  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
3565  for (iPass=0;iPass<3;iPass++) {
3566    last[iPass]=new int[numberNonLinearColumns];
3567    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3568      last[iPass][jNon]=0;
3569  }
3570  int numberZeroPasses=0;
3571  bool zeroTargetDrop=false;
3572  double * gradient = new double [numberColumns_];
3573#define SMALL_FIX 0.0
3574  for (iPass=0;iPass<numberPasses;iPass++) {
3575    objectiveOffset2 = objectiveOffset;
3576    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3577      iColumn=listNonLinearColumn[jNon];
3578      if (solution[iColumn]<trueLower[jNon])
3579        solution[iColumn]=trueLower[jNon];
3580      else if (solution[iColumn]>trueUpper[jNon])
3581        solution[iColumn]=trueUpper[jNon];
3582      columnLower[iColumn]=CoinMax(solution[iColumn]
3583                                   -trust[jNon],
3584                                   trueLower[jNon]);
3585      if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
3586        columnLower[iColumn]=SMALL_FIX;
3587      columnUpper[iColumn]=CoinMin(solution[iColumn]
3588                                   +trust[jNon],
3589                                   trueUpper[jNon]);
3590      if (!trueLower[jNon]) 
3591        columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
3592                                       columnLower[iColumn]+SMALL_FIX);
3593      if (!trueLower[jNon]&&tryFix&&
3594          columnLower[iColumn]==SMALL_FIX&&
3595          columnUpper[iColumn]<3.0*SMALL_FIX) {
3596        columnLower[iColumn]=0.0;
3597        solution[iColumn]=0.0;
3598        columnUpper[iColumn]=0.0;
3599        printf("fixing %d\n",iColumn);
3600      }
3601    }
3602    // redo matrix
3603    double offset;
3604    CoinPackedMatrix newMatrix(saveMatrix);
3605    // get matrix data pointers
3606    column = newMatrix.getIndices();
3607    rowStart = newMatrix.getVectorStarts();
3608    rowLength = newMatrix.getVectorLengths(); 
3609    // make sure x updated
3610    if (numberConstraints)
3611      constraints[0]->newXValues();
3612    else
3613      trueObjective->newXValues();
3614    double * changeableElement = newMatrix.getMutableElements();
3615    if (trueObjective) {
3616      memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
3617             numberColumns_*sizeof(double));
3618    } else {
3619      memcpy(objective,objective_->gradient(this,solution,offset,true,2),
3620             numberColumns_*sizeof(double));
3621    }
3622    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3623      ClpConstraint * constraint = constraints[iConstraint];
3624      int iRow = constraint->rowNumber();
3625      double functionValue;
3626#ifndef NDEBUG
3627      int numberErrors = 
3628#endif
3629      constraint->gradient(&newModel,solution,gradient,functionValue,offset);
3630      assert (!numberErrors);
3631      double dualValue = newModel.dualRowSolution()[iRow];
3632      int numberCoefficients = constraint->numberCoefficients();
3633      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
3634        int iColumn = column[j];
3635        changeableElement[j] = gradient[iColumn];
3636        //objective[iColumn] -= dualValue*gradient[iColumn];
3637        gradient[iColumn]=0.0;
3638      }
3639      for (int k=0;k<numberColumns_;k++)
3640        assert (!gradient[k]);
3641      if (rowLower_[iRow]>-1.0e20)
3642        rowLower[iRow] = rowLower_[iRow] - offset;
3643      if (rowUpper_[iRow]<1.0e20)
3644        rowUpper[iRow] = rowUpper_[iRow] - offset;
3645    }
3646    // Replace matrix
3647    // Get a column copy in standard format
3648    CoinPackedMatrix * columnCopy = new CoinPackedMatrix();;
3649    columnCopy->reverseOrderedCopyOf(newMatrix);
3650    newModel.replaceMatrix(columnCopy,true);
3651    // solve
3652    newModel.primal(1);
3653    if (newModel.status()==1) { 
3654      // Infeasible!
3655      newModel.allSlackBasis();
3656      newModel.primal();
3657      newModel.writeMps("infeas.mps");
3658      assert(!newModel.status());
3659    }
3660    double sumInfeas=0.0;
3661    int numberInfeas=0;
3662    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
3663      if (solution[iColumn]>1.0e-8) {
3664        numberInfeas++;
3665        sumInfeas += solution[iColumn];
3666      }
3667    }
3668    printf("%d artificial infeasibilities - summing to %g\n",
3669           numberInfeas,sumInfeas);
3670    double infValue=0.0;
3671    double objValue=0.0;
3672    // make sure x updated
3673    if (numberConstraints)
3674      constraints[0]->newXValues();
3675    else
3676      trueObjective->newXValues();
3677    if (trueObjective) {
3678      objValue=trueObjective->objectiveValue(this,solution);
3679      printf("objective offset %g\n",offset);
3680      objectiveOffset2 =objectiveOffset+offset; // ? sign
3681      newModel.setObjectiveOffset(objectiveOffset2);
3682    } else {
3683      objValue=objective_->objectiveValue(this,solution);
3684    }
3685    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
3686      ClpConstraint * constraint = constraints[iConstraint];
3687      int iRow = constraint->rowNumber();
3688      double functionValue=constraint->functionValue(this,solution);
3689      double dualValue = newModel.dualRowSolution()[iRow];
3690      if (numberConstraints<50)
3691        printf("For row %d current value is %g (activity %g) , dual is %g - offset %g\n",iRow,functionValue,
3692               newModel.primalRowSolution()[iRow],
3693               dualValue,offset);
3694      if (functionValue<rowLower_[iRow]-1.0e-5) {
3695        double under = rowLower_[iRow]-functionValue;
3696        infValue += under;
3697      } else if (functionValue>rowUpper_[iRow]+1.0e-5) {
3698        double over = functionValue-rowUpper_[iRow];
3699        infValue += over;
3700      }
3701    }
3702    if (infValue)
3703      printf("Sum infeasibilities %g ",infValue);
3704    if (infValue<0.1)
3705      infValue=0.0;
3706    infValue *= penalty;
3707    if (infValue)
3708      printf("Infeasible obj %g ",infValue);
3709    if (objectiveOffset2)
3710      printf("offset2 %g\n",objectiveOffset2);
3711    objValue -= objectiveOffset2;
3712    printf("True objective %g\n",objValue);
3713    objValue += infValue;
3714    if (iPass) {
3715      double drop = lastObjective-objValue;
3716      std::cout<<"True drop was "<<drop<<std::endl;
3717      if (drop<-0.05*fabs(objValue)-1.0e-4) {
3718        // pretty bad - go back and halve
3719        memcpy(solution,saveSolution,numberColumns2*sizeof(double));
3720        memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2,
3721               numberRows*sizeof(double));
3722        memcpy(newModel.statusArray(),saveStatus,
3723               numberColumns2+numberRows);
3724        for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3725          if (trust[jNon]>0.1)
3726            trust[jNon] *= 0.5;
3727          else
3728            trust[jNon] *= 0.9;
3729       
3730        printf("** Halving trust\n");
3731        objValue=lastObjective;
3732        continue;
3733      } else if ((iPass%25)==-1) {
3734        for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3735          trust[jNon] *= 2.0;
3736        printf("** Doubling trust\n");
3737      }
3738      int * temp=last[2];
3739      last[2]=last[1];
3740      last[1]=last[0];
3741      last[0]=temp;
3742      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3743        iColumn=listNonLinearColumn[jNon];
3744        double change = solution[iColumn]-saveSolution[iColumn];
3745        if (change<-1.0e-5) {
3746          if (fabs(change+trust[jNon])<1.0e-5) 
3747            temp[jNon]=-1;
3748          else
3749            temp[jNon]=-2;
3750        } else if(change>1.0e-5) {
3751          if (fabs(change-trust[jNon])<1.0e-5) 
3752            temp[jNon]=1;
3753          else
3754            temp[jNon]=2;
3755        } else {
3756          temp[jNon]=0;
3757        }
3758      } 
3759      double maxDelta=0.0;
3760      double smallestTrust=1.0e31;
3761      double smallestNonLinearGap=1.0e31;
3762      double smallestGap=1.0e31;
3763      bool increasing=false;
3764      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3765        double gap = columnUpper[iColumn]-columnLower[iColumn];
3766        assert (gap>=0.0);
3767        if (gap)
3768          smallestGap = CoinMin(smallestGap,gap);
3769      }
3770      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3771        iColumn=listNonLinearColumn[jNon];
3772        double gap = columnUpper[iColumn]-columnLower[iColumn];
3773        assert (gap>=0.0);
3774        if (gap) {
3775          smallestNonLinearGap = CoinMin(smallestNonLinearGap,gap);
3776          if (gap<1.0e-7&&iPass==1) {
3777            printf("Small gap %d %d %g %g %g\n",
3778                   jNon,iColumn,columnLower[iColumn],columnUpper[iColumn],
3779                   gap);
3780            //trueUpper[jNon]=trueLower[jNon];
3781            //columnUpper[iColumn]=columnLower[iColumn];
3782          }
3783        }
3784        maxDelta = CoinMax(maxDelta,
3785                       fabs(solution[iColumn]-saveSolution[iColumn]));
3786        if (last[0][jNon]*last[1][jNon]<0) {
3787          // halve
3788          if (trust[jNon]>1.0)
3789            trust[jNon] *= 0.5;
3790          else
3791            trust[jNon] *= 0.7;
3792        } else {
3793          // ? only increase if +=1 ?
3794          if (last[0][jNon]==last[1][jNon]&&
3795              last[0][jNon]==last[2][jNon]&&
3796              last[0][jNon]) {
3797            trust[jNon] *= 1.8; 
3798            increasing=true;
3799          }
3800        }
3801        smallestTrust = CoinMin(smallestTrust,trust[jNon]);
3802      }
3803      std::cout<<"largest delta is "<<maxDelta
3804               <<", smallest trust is "<<smallestTrust
3805               <<", smallest gap is "<<smallestGap
3806               <<", smallest nonlinear gap is "<<smallestNonLinearGap
3807               <<std::endl;
3808      if (maxDelta<deltaTolerance&&!increasing&&iPass>100) {
3809        numberZeroPasses++;
3810        if (numberZeroPasses==3) {
3811          if (tryFix) {
3812            std::cout<<"Exiting"<<std::endl;
3813            break;
3814          } else {
3815            tryFix=true;
3816            for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
3817              trust[jNon] = CoinMax(trust[jNon],1.0e-1);
3818            numberZeroPasses=0;
3819          }
3820        }
3821      } else {
3822        numberZeroPasses=0;
3823      }
3824    }
3825    memcpy(saveSolution,solution,numberColumns2*sizeof(double));
3826    memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(),
3827           numberRows*sizeof(double));
3828    memcpy(saveStatus,newModel.statusArray(),
3829           numberColumns2+numberRows);
3830   
3831    targetDrop=0.0;
3832    if (iPass) {
3833      // get reduced costs
3834      const double * pi = newModel.dualRowSolution();
3835      newModel.matrix()->transposeTimes(pi,
3836                                        newModel.dualColumnSolution());
3837      double * r = newModel.dualColumnSolution();
3838      for (iColumn=0;iColumn<numberColumns_;iColumn++) 
3839        r[iColumn] = objective[iColumn]-r[iColumn];
3840      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3841        iColumn=listNonLinearColumn[jNon];
3842        double dj = r[iColumn];
3843        if (dj<-1.0e-6) {
3844          double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
3845          //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
3846          //double drop2 = -dj*(upper-solution[iColumn]);
3847#if 0
3848          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2&&iPass>100))
3849            printf("Big drop %d %g %g %g %g T %g\n",
3850                   iColumn,columnLower[iColumn],solution[iColumn],
3851                   columnUpper[iColumn],dj,trueUpper[jNon]);
3852#endif
3853          targetDrop += drop;
3854          if (dj<-1.0e-1&&trust[jNon]<1.0e-3
3855              &&trueUpper[jNon]-solution[iColumn]>1.0e-2) {
3856            trust[jNon] *= 1.5;
3857            //printf("Increasing trust on %d to %g\n",
3858            //     iColumn,trust[jNon]);
3859          }
3860        } else if (dj>1.0e-6) {
3861          double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
3862          //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
3863          //double drop2 = -dj*(lower-solution[iColumn]);
3864#if 0
3865          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2))
3866            printf("Big drop %d %g %g %g %g T %g\n",
3867                   iColumn,columnLower[iColumn],solution[iColumn],
3868                   columnUpper[iColumn],dj,trueLower[jNon]);
3869#endif
3870          targetDrop += drop;
3871          if (dj>1.0e-1&&trust[jNon]<1.0e-3
3872              &&solution[iColumn]-trueLower[jNon]>1.0e-2) {
3873            trust[jNon] *= 1.5;
3874            printf("Increasing trust on %d to %g\n",
3875                   iColumn,trust[jNon]);
3876          }
3877        }
3878      }
3879    }
3880    std::cout<<"Pass - "<<iPass
3881             <<", target drop is "<<targetDrop
3882             <<std::endl;
3883    if (iPass>1&&targetDrop<1.0e-5&&zeroTargetDrop)
3884      break;
3885    if (iPass>1&&targetDrop<1.0e-5)
3886      zeroTargetDrop = true;
3887    else
3888      zeroTargetDrop = false;
3889    //if (iPass==5)
3890    //newModel.setLogLevel(63);
3891    lastObjective = objValue;
3892    // take out when ClpPackedMatrix changed
3893    //newModel.scaling(false);
3894#if 0
3895    CoinMpsIO writer;
3896    writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
3897                      newModel.getColLower(), newModel.getColUpper(),
3898                      newModel.getObjCoefficients(),
3899                      (const char*) 0 /*integrality*/,
3900                      newModel.getRowLower(), newModel.getRowUpper(),
3901                      NULL,NULL);
3902    writer.writeMps("xx.mps");
3903#endif
3904  }
3905  delete [] saveSolution;
3906  delete [] saveStatus;
3907  for (iPass=0;iPass<3;iPass++) 
3908    delete [] last[iPass];
3909  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
3910    iColumn=listNonLinearColumn[jNon];
3911    columnLower[iColumn]=trueLower[jNon];
3912    columnUpper[iColumn]=trueUpper[jNon];
3913  }
3914  delete [] trust;
3915  delete [] trueUpper;
3916  delete [] trueLower;
3917  // Simplest way to get true row activity ?
3918  double * rowActivity = newModel.primalRowSolution();
3919  for (iRow=0;iRow<numberRows;iRow++) {
3920    double difference;
3921    if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow]))
3922      difference = rowLower_[iRow]-rowLower[iRow];
3923    else
3924      difference = rowUpper_[iRow]-rowUpper[iRow];
3925    rowLower[iRow]=rowLower_[iRow];
3926    rowUpper[iRow]=rowUpper_[iRow];
3927    if (difference) {
3928      if (numberRows<50)
3929        printf("For row %d activity changes from %g to %g\n",
3930               iRow,rowActivity[iRow],rowActivity[iRow]+difference);
3931      rowActivity[iRow]+= difference;
3932    }
3933  }
3934  delete [] listNonLinearColumn;
3935  delete [] gradient;
3936  printf("solution still in newModel - do objective etc!\n");
3937  objectiveValue_=newModel.objectiveValue();
3938  numberIterations_ =newModel.numberIterations();
3939  problemStatus_ =newModel.problemStatus();
3940  secondaryStatus_ =newModel.secondaryStatus();
3941  memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double));
3942  // should do status region
3943  CoinZeroN(rowActivity_,numberRows_);
3944  matrix_->times(1.0,columnActivity_,rowActivity_);
3945  return 0;
3946}
Note: See TracBrowser for help on using the repository browser.