source: trunk/ClpSimplexNonlinear.cpp @ 437

Last change on this file since 437 was 414, checked in by forrest, 16 years ago

barrier stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 85.8 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <math.h>
7#include "CoinHelperFunctions.hpp"
8#include "ClpHelperFunctions.hpp"
9#include "ClpSimplexNonlinear.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpNonLinearCost.hpp"
12#include "ClpLinearObjective.hpp"
13#include "CoinPackedMatrix.hpp"
14#include "CoinIndexedVector.hpp"
15#include "ClpPrimalColumnPivot.hpp"
16#include "ClpMessage.hpp"
17#include "ClpEventHandler.hpp"
18#include <cfloat>
19#include <cassert>
20#include <string>
21#include <stdio.h>
22#include <iostream>
23#ifndef NDEBUG
24#define CLP_DEBUG 1
25#endif
26// primal
27int ClpSimplexNonlinear::primal ()
28{
29
30  int ifValuesPass=1;
31  algorithm_ = +3;
32
33  // save data
34  ClpDataSave data = saveData();
35  matrix_->refresh(this); // make sure matrix okay
36
37  double bestObjectiveWhenFlagged=COIN_DBL_MAX;
38  int pivotMode=15;
39
40  // initialize - maybe values pass and algorithm_ is +1
41  // true does something (? perturbs)
42  if (!startup(true)) {
43   
44    // Set average theta
45    nonLinearCost_->setAverageTheta(1.0e3);
46    int lastCleaned=0; // last time objective or bounds cleaned up
47   
48    // Say no pivot has occurred (for steepest edge and updates)
49    pivotRow_=-2;
50   
51    // This says whether to restore things etc
52    int factorType=0;
53    // Start check for cycles
54    progress_->startCheck();
55    /*
56      Status of problem:
57      0 - optimal
58      1 - infeasible
59      2 - unbounded
60      -1 - iterating
61      -2 - factorization wanted
62      -3 - redo checking without factorization
63      -4 - looks infeasible
64      -5 - looks unbounded
65    */
66    while (problemStatus_<0) {
67      int iRow,iColumn;
68      // clear
69      for (iRow=0;iRow<4;iRow++) {
70        rowArray_[iRow]->clear();
71      }   
72     
73      for (iColumn=0;iColumn<2;iColumn++) {
74        columnArray_[iColumn]->clear();
75      }   
76     
77      // give matrix (and model costs and bounds a chance to be
78      // refreshed (normally null)
79      matrix_->refresh(this);
80      // If getting nowhere - why not give it a kick
81      // If we have done no iterations - special
82      if (lastGoodIteration_==numberIterations_&&factorType)
83        factorType=3;
84     
85      // may factorize, checks if problem finished
86      if (objective_->type()>1&&lastFlaggedIteration_>=0&&
87          numberIterations_>lastFlaggedIteration_+507) {
88        unflag();
89        lastFlaggedIteration_=numberIterations_;
90        if (pivotMode>=10) {
91          pivotMode--;
92#ifdef CLP_DEBUG
93          if (handler_->logLevel()&32) 
94            printf("pivot mode now %d\n",pivotMode);
95#endif
96          if (pivotMode==9) 
97            pivotMode=0;        // switch off fast attempt
98        }
99      }
100      statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,
101                              bestObjectiveWhenFlagged);
102     
103      // Say good factorization
104      factorType=1;
105     
106      // Say no pivot has occurred (for steepest edge and updates)
107      pivotRow_=-2;
108
109      // exit if victory declared
110      if (problemStatus_>=0)
111        break;
112     
113      // test for maximum iterations
114      if (hitMaximumIterations()||(ifValuesPass==2&&firstFree_<0)) {
115        problemStatus_=3;
116        break;
117      }
118
119      if (firstFree_<0) {
120        if (ifValuesPass) {
121          // end of values pass
122          ifValuesPass=0;
123          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
124          if (status>=0) {
125            problemStatus_=5;
126            secondaryStatus_=ClpEventHandler::endOfValuesPass;
127            break;
128          }
129        }
130      }
131      // Check event
132      {
133        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
134        if (status>=0) {
135          problemStatus_=5;
136          secondaryStatus_=ClpEventHandler::endOfFactorization;
137          break;
138        }
139      }
140      // Iterate
141      whileIterating(pivotMode);
142    }
143  }
144  // if infeasible get real values
145  if (problemStatus_==1) {
146    infeasibilityCost_=0.0;
147    createRim(7);
148    nonLinearCost_->checkInfeasibilities(0.0);
149    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
150    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
151    // and get good feasible duals
152    computeDuals(NULL);
153  }
154  // correct objective value
155  if (numberColumns_)
156    objectiveValue_ = nonLinearCost_->feasibleCost()+objective_->nonlinearOffset();
157  // clean up
158  unflag();
159  finish();
160  restoreData(data);
161  return problemStatus_;
162}
163/*  Refactorizes if necessary
164    Checks if finished.  Updates status.
165    lastCleaned refers to iteration at which some objective/feasibility
166    cleaning too place.
167   
168    type - 0 initial so set up save arrays etc
169    - 1 normal -if good update save
170    - 2 restoring from saved
171*/
172void 
173ClpSimplexNonlinear::statusOfProblemInPrimal(int & lastCleaned, int type,
174                                                 ClpSimplexProgress * progress,
175                                                 bool doFactorization,
176                                                 double & bestObjectiveWhenFlagged)
177{
178  int dummy; // for use in generalExpanded
179  if (type==2) {
180    // trouble - restore solution
181    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
182    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
183           numberRows_*sizeof(double));
184    memcpy(columnActivityWork_,savedSolution_ ,
185           numberColumns_*sizeof(double));
186    // restore extra stuff
187    matrix_->generalExpanded(this,6,dummy);
188    forceFactorization_=1; // a bit drastic but ..
189    pivotRow_=-1; // say no weights update
190    changeMade_++; // say change made
191  }
192  int saveThreshold = factorization_->sparseThreshold();
193  int tentativeStatus = problemStatus_;
194  int numberThrownOut=1; // to loop round on bad factorization in values pass
195  while (numberThrownOut) {
196    if (problemStatus_>-3||problemStatus_==-4) {
197      // factorize
198      // later on we will need to recover from singularities
199      // also we could skip if first time
200      // do weights
201      // This may save pivotRow_ for use
202      if (doFactorization)
203        primalColumnPivot_->saveWeights(this,1);
204     
205      if (type&&doFactorization) {
206        // is factorization okay?
207        int factorStatus = internalFactorize(1);
208        if (factorStatus) {
209          if (type!=1||largestPrimalError_>1.0e3
210              ||largestDualError_>1.0e3) {
211            // was ||largestDualError_>1.0e3||objective_->type()>1) {
212            // switch off dense
213            int saveDense = factorization_->denseThreshold();
214            factorization_->setDenseThreshold(0);
215            // make sure will do safe factorization
216            pivotVariable_[0]=-1;
217            internalFactorize(2);
218            factorization_->setDenseThreshold(saveDense);
219            // Go to safe
220            factorization_->pivotTolerance(0.99);
221            // restore extra stuff
222            matrix_->generalExpanded(this,6,dummy);
223          } else {
224            // no - restore previous basis
225            memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
226            memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
227                   numberRows_*sizeof(double));
228            memcpy(columnActivityWork_,savedSolution_ ,
229                   numberColumns_*sizeof(double));
230            // restore extra stuff
231            matrix_->generalExpanded(this,6,dummy);
232            matrix_->generalExpanded(this,5,dummy);
233            forceFactorization_=1; // a bit drastic but ..
234            type = 2;
235            // Go to safe
236            factorization_->pivotTolerance(0.99);
237            if (internalFactorize(1)!=0)
238               largestPrimalError_=1.0e4; // force other type
239          }
240          // flag incoming
241          if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
242            setFlagged(sequenceIn_);
243            saveStatus_[sequenceIn_]=status_[sequenceIn_];
244          }
245          changeMade_++; // say change made
246        }
247      }
248      if (problemStatus_!=-4)
249        problemStatus_=-3;
250    }
251    // at this stage status is -3 or -5 if looks unbounded
252    // get primal and dual solutions
253    // put back original costs and then check
254    createRim(4);
255    // May need to do more if column generation
256    dummy=4;
257    matrix_->generalExpanded(this,9,dummy);
258    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
259    if (numberThrownOut) {
260      problemStatus_=tentativeStatus;
261      doFactorization=true;
262    }
263  }
264  // Double check reduced costs if no action
265  if (progress->lastIterationNumber(0)==numberIterations_) {
266    if (primalColumnPivot_->looksOptimal()) {
267      numberDualInfeasibilities_ = 0;
268      sumDualInfeasibilities_ = 0.0;
269    }
270  }
271  // Check if looping
272  int loop;
273  if (type!=2) 
274    loop = progress->looping();
275  else
276    loop=-1;
277  if (loop>=0) {
278    if (!problemStatus_) {
279      // declaring victory
280      numberPrimalInfeasibilities_ = 0;
281      sumPrimalInfeasibilities_ = 0.0;
282    } else {
283      problemStatus_ = loop; //exit if in loop
284      problemStatus_ = 10; // instead - try other algorithm
285    }
286    problemStatus_ = 10; // instead - try other algorithm
287    return ;
288  } else if (loop<-1) {
289    // Is it time for drastic measures
290    if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
291        progress->oddState()<10&&progress->oddState()>=0) {
292      progress->newOddState();
293      nonLinearCost_->zapCosts();
294    }
295    // something may have changed
296    gutsOfSolution(NULL,NULL,true);
297  }
298  // If progress then reset costs
299  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
300    createRim(4,false); // costs back
301    delete nonLinearCost_;
302    nonLinearCost_ = new ClpNonLinearCost(this);
303    progress->endOddState();
304    gutsOfSolution(NULL,NULL,true);
305  }
306  // Flag to say whether to go to dual to clean up
307  bool goToDual=false;
308  // really for free variables in
309  //if((progressFlag_&2)!=0)
310  //problemStatus_=-1;;
311  progressFlag_ = 0; //reset progress flag
312
313  handler_->message(CLP_SIMPLEX_STATUS,messages_)
314    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
315  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
316    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
317  handler_->printing(sumDualInfeasibilities_>0.0)
318    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
319  handler_->printing(numberDualInfeasibilitiesWithoutFree_
320                     <numberDualInfeasibilities_)
321                       <<numberDualInfeasibilitiesWithoutFree_;
322  handler_->message()<<CoinMessageEol;
323  if (!primalFeasible()) {
324    nonLinearCost_->checkInfeasibilities(primalTolerance_);
325    gutsOfSolution(NULL,NULL,true);
326    nonLinearCost_->checkInfeasibilities(primalTolerance_);
327  }
328  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
329  if (trueInfeasibility>1.0) {
330    // If infeasibility going up may change weights
331    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
332    if(progress->lastInfeasibility()<testValue) {
333      if (infeasibilityCost_<1.0e14) {
334        infeasibilityCost_ *= 1.5;
335        printf("increasing weight to %g\n",infeasibilityCost_);
336        gutsOfSolution(NULL,NULL,true);
337      }
338    }
339  }
340  // we may wish to say it is optimal even if infeasible
341  bool alwaysOptimal = (specialOptions_&1)!=0;
342  // give code benefit of doubt
343  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
344      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
345    // say optimal (with these bounds etc)
346    numberDualInfeasibilities_ = 0;
347    sumDualInfeasibilities_ = 0.0;
348    numberPrimalInfeasibilities_ = 0;
349    sumPrimalInfeasibilities_ = 0.0;
350  }
351  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
352  if (dualFeasible()||problemStatus_==-4) {
353    // see if extra helps
354    if (nonLinearCost_->numberInfeasibilities()&&
355         (nonLinearCost_->sumInfeasibilities()>1.0e-3||sumOfRelaxedPrimalInfeasibilities_)
356        &&!alwaysOptimal) {
357      //may need infeasiblity cost changed
358      // we can see if we can construct a ray
359      // make up a new objective
360      double saveWeight = infeasibilityCost_;
361      // save nonlinear cost as we are going to switch off costs
362      ClpNonLinearCost * nonLinear = nonLinearCost_;
363      // do twice to make sure Primal solution has settled
364      // put non-basics to bounds in case tolerance moved
365      // put back original costs
366      createRim(4);
367      nonLinearCost_->checkInfeasibilities(0.0);
368      gutsOfSolution(NULL,NULL,true);
369
370      infeasibilityCost_=1.0e100;
371      // put back original costs
372      createRim(4);
373      nonLinearCost_->checkInfeasibilities(primalTolerance_);
374      // may have fixed infeasibilities - double check
375      if (nonLinearCost_->numberInfeasibilities()==0) {
376        // carry on
377        problemStatus_ = -1;
378        infeasibilityCost_=saveWeight;
379        nonLinearCost_->checkInfeasibilities(primalTolerance_);
380      } else {
381        nonLinearCost_=NULL;
382        // scale
383        int i;
384        for (i=0;i<numberRows_+numberColumns_;i++) 
385          cost_[i] *= 1.0e-95;
386        gutsOfSolution(NULL,NULL,false);
387        nonLinearCost_=nonLinear;
388        infeasibilityCost_=saveWeight;
389        if ((infeasibilityCost_>=1.0e18||
390             numberDualInfeasibilities_==0)&&perturbation_==101) {
391          goToDual=unPerturb(); // stop any further perturbation
392          if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
393            goToDual=false;
394          nonLinearCost_->checkInfeasibilities(primalTolerance_);
395          numberDualInfeasibilities_=1; // carry on
396          problemStatus_=-1;
397        }
398        if (infeasibilityCost_>=1.0e20||
399            numberDualInfeasibilities_==0) {
400          // we are infeasible - use as ray
401          delete [] ray_;
402          ray_ = new double [numberRows_];
403          memcpy(ray_,dual_,numberRows_*sizeof(double));
404          // and get feasible duals
405          infeasibilityCost_=0.0;
406          createRim(4);
407          nonLinearCost_->checkInfeasibilities(primalTolerance_);
408          gutsOfSolution(NULL,NULL,true);
409          // so will exit
410          infeasibilityCost_=1.0e30;
411          // reset infeasibilities
412          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
413          numberPrimalInfeasibilities_=
414            nonLinearCost_->numberInfeasibilities();
415        }
416        if (infeasibilityCost_<1.0e20) {
417          infeasibilityCost_ *= 5.0;
418          changeMade_++; // say change made
419          unflag();
420          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
421            <<infeasibilityCost_
422            <<CoinMessageEol;
423          // put back original costs and then check
424          createRim(4);
425          nonLinearCost_->checkInfeasibilities(0.0);
426          gutsOfSolution(NULL,NULL,true);
427          problemStatus_=-1; //continue
428          goToDual=false;
429        } else {
430          // say infeasible
431          problemStatus_ = 1;
432        }
433      }
434    } else {
435      // may be optimal
436      if (perturbation_==101) {
437        goToDual=unPerturb(); // stop any further perturbation
438        lastCleaned=-1; // carry on
439      }
440      bool unflagged = unflag()!=0;
441      if ( lastCleaned!=numberIterations_||unflagged) {
442        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
443          <<primalTolerance_
444          <<CoinMessageEol;
445        double current = nonLinearCost_->feasibleReportCost();
446        if (numberTimesOptimal_<4) {
447          if (bestObjectiveWhenFlagged<=current) {
448            numberTimesOptimal_++;
449#ifdef CLP_DEBUG
450            if (handler_->logLevel()&32) 
451              printf("%d times optimal, current %g best %g\n",numberTimesOptimal_,
452                     current,bestObjectiveWhenFlagged);
453#endif
454          } else {
455            bestObjectiveWhenFlagged=current; 
456          }
457          changeMade_++; // say change made
458          if (numberTimesOptimal_==1) {
459            // better to have small tolerance even if slower
460            factorization_->zeroTolerance(1.0e-15);
461          }
462          lastCleaned=numberIterations_;
463          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
464            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
465              <<CoinMessageEol;
466          double oldTolerance = primalTolerance_;
467          primalTolerance_=dblParam_[ClpPrimalTolerance];
468          // put back original costs and then check
469          createRim(4);
470          nonLinearCost_->checkInfeasibilities(oldTolerance);
471          gutsOfSolution(NULL,NULL,true);
472          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
473              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
474            // say optimal (with these bounds etc)
475            numberDualInfeasibilities_ = 0;
476            sumDualInfeasibilities_ = 0.0;
477            numberPrimalInfeasibilities_ = 0;
478            sumPrimalInfeasibilities_ = 0.0;
479          }
480          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
481            problemStatus_=0;
482          else
483            problemStatus_ = -1;
484        } else {
485          problemStatus_=0; // optimal
486          if (lastCleaned<numberIterations_) {
487            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
488              <<CoinMessageEol;
489          }
490        }
491      } else {
492        problemStatus_=0; // optimal
493      }
494    }
495  } else {
496    // see if looks unbounded
497    if (problemStatus_==-5) {
498      if (nonLinearCost_->numberInfeasibilities()) {
499        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
500          // back off weight
501          infeasibilityCost_ = 1.0e13;
502          unPerturb(); // stop any further perturbation
503        }
504        //we need infeasiblity cost changed
505        if (infeasibilityCost_<1.0e20) {
506          infeasibilityCost_ *= 5.0;
507          changeMade_++; // say change made
508          unflag();
509          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
510            <<infeasibilityCost_
511            <<CoinMessageEol;
512          // put back original costs and then check
513          createRim(4);
514          gutsOfSolution(NULL,NULL,true);
515          problemStatus_=-1; //continue
516        } else {
517          // say unbounded
518          problemStatus_ = 2;
519        }
520      } else {
521        // say unbounded
522        problemStatus_ = 2;
523      } 
524    } else {
525      if(type==3&&problemStatus_!=-5)
526        unflag(); // odd
527      // carry on
528      problemStatus_ = -1;
529    }
530  }
531  // save extra stuff
532  matrix_->generalExpanded(this,5,dummy);
533  if (type==0||type==1) {
534    if (type!=1||!saveStatus_) {
535      // create save arrays
536      delete [] saveStatus_;
537      delete [] savedSolution_;
538      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
539      savedSolution_ = new double [numberRows_+numberColumns_];
540    }
541    // save arrays
542    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
543    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
544           numberRows_*sizeof(double));
545    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
546  }
547  if (doFactorization) {
548    // restore weights (if saved) - also recompute infeasibility list
549    if (tentativeStatus>-3) 
550      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
551    else
552      primalColumnPivot_->saveWeights(this,3);
553    if (saveThreshold) {
554      // use default at present
555      factorization_->sparseThreshold(0);
556      factorization_->goSparse();
557    }
558  }
559  if (problemStatus_<0&&!changeMade_) {
560    problemStatus_=4; // unknown
561  }
562  lastGoodIteration_ = numberIterations_;
563  //if (goToDual)
564  //problemStatus_=10; // try dual
565}
566/*
567  Reasons to come out:
568  -1 iterations etc
569  -2 inaccuracy
570  -3 slight inaccuracy (and done iterations)
571  -4 end of values pass and done iterations
572  +0 looks optimal (might be infeasible - but we will investigate)
573  +2 looks unbounded
574  +3 max iterations
575*/
576int
577ClpSimplexNonlinear::whileIterating(int & pivotMode)
578{
579  // Say if values pass
580  //int ifValuesPass=(firstFree_>=0) ? 1 : 0;
581  int ifValuesPass=1;
582  int returnCode=-1;
583  // status stays at -1 while iterating, >=0 finished, -2 to invert
584  // status -3 to go to top without an invert
585  int numberInterior=0;
586  int nextUnflag=10;
587  int nextUnflagIteration=numberIterations_+10;
588  // get two arrays
589  double * array1 = new double[numberRows_+numberColumns_];
590  double solutionError=-1.0;
591  while (problemStatus_==-1) {
592    int result;
593    rowArray_[1]->clear();
594    //#define CLP_DEBUG
595#if CLP_DEBUG > 1
596    rowArray_[0]->checkClear();
597    rowArray_[1]->checkClear();
598    rowArray_[2]->checkClear();
599    rowArray_[3]->checkClear();
600    columnArray_[0]->checkClear();
601#endif
602    if (numberInterior>=5) {
603      // this can go when we have better minimization
604      if (pivotMode<10)
605        pivotMode=1;
606      unflag();
607#ifdef CLP_DEBUG
608          if (handler_->logLevel()&32) 
609            printf("interior unflag\n");
610#endif
611      numberInterior=0;
612      nextUnflag=10;
613      nextUnflagIteration=numberIterations_+10;
614    } else {
615      if (numberInterior>nextUnflag&&
616          numberIterations_>nextUnflagIteration) {
617        nextUnflagIteration=numberIterations_+10;
618        nextUnflag += 10;
619        unflag();
620#ifdef CLP_DEBUG
621        if (handler_->logLevel()&32) 
622          printf("unflagging as interior\n");
623#endif
624      }
625    } 
626    pivotRow_=-1;
627    result = pivotColumn(rowArray_[3],rowArray_[0],
628                         columnArray_[0],rowArray_[1],pivotMode,solutionError,
629                         array1);
630    if (result) {
631      if (result==3) 
632        break; // null vector not accurate
633#ifdef CLP_DEBUG
634      if (handler_->logLevel()&32) {
635        double currentObj;
636        double thetaObj;
637        double predictedObj;
638        objective_->stepLength(this,solution_,solution_,0.0,
639                               currentObj,predictedObj,thetaObj);
640        printf("obj %g after interior move\n",currentObj);
641      }
642#endif
643      // just move and try again
644      if (pivotMode<10) {
645        pivotMode=result-1;
646        numberInterior++;
647      }
648      continue;
649    } else {
650      if (pivotMode<10) {
651        if (theta_>0.001)
652          pivotMode=0;
653        else if (pivotMode==2)
654          pivotMode=1;
655      }
656      numberInterior=0;
657      nextUnflag=10;
658      nextUnflagIteration=numberIterations_+10;
659    }
660    sequenceOut_=-1;
661    rowArray_[1]->clear();
662    if (sequenceIn_>=0) {
663      // we found a pivot column
664      assert (!flagged(sequenceIn_));
665#ifdef CLP_DEBUG
666      if ((handler_->logLevel()&32)) {
667        char x = isColumn(sequenceIn_) ? 'C' :'R';
668        std::cout<<"pivot column "<<
669          x<<sequenceWithin(sequenceIn_)<<std::endl;
670      }
671#endif
672      // do second half of iteration
673      if (pivotRow_<0&&theta_<1.0e-8) {
674        assert (sequenceIn_>=0);
675        returnCode = pivotResult(ifValuesPass);
676      } else {
677        // specialized code
678        returnCode = pivotNonlinearResult();
679        //printf("odd pivrow %d\n",sequenceOut_);
680        if (sequenceOut_>=0&&theta_<1.0e-5) {
681          if (getStatus(sequenceOut_)!=isFixed) {
682            if (getStatus(sequenceOut_)==atUpperBound)
683              solution_[sequenceOut_]=upper_[sequenceOut_];
684            else if (getStatus(sequenceOut_)==atLowerBound)
685              solution_[sequenceOut_]=lower_[sequenceOut_];
686            setFlagged(sequenceOut_);
687          }
688        }
689      }
690      if (returnCode<-1&&returnCode>-5) {
691        problemStatus_=-2; //
692      } else if (returnCode==-5) {
693        // something flagged - continue;
694      } else if (returnCode==2) {
695        problemStatus_=-5; // looks unbounded
696      } else if (returnCode==4) {
697        problemStatus_=-2; // looks unbounded but has iterated
698      } else if (returnCode!=-1) {
699        assert(returnCode==3);
700        problemStatus_=3;
701      }
702    } else {
703      // no pivot column
704#ifdef CLP_DEBUG
705      if (handler_->logLevel()&32)
706        printf("** no column pivot\n");
707#endif
708      if (pivotMode<10) {
709        // looks optimal
710        primalColumnPivot_->setLooksOptimal(true);
711      } else {
712        pivotMode--;
713#ifdef CLP_DEBUG
714        if (handler_->logLevel()&32) 
715          printf("pivot mode now %d\n",pivotMode);
716#endif
717        if (pivotMode==9) 
718          pivotMode=0;  // switch off fast attempt
719        unflag();
720      }
721      if (nonLinearCost_->numberInfeasibilities())
722        problemStatus_=-4; // might be infeasible
723      returnCode=0;
724      break;
725    }
726  }
727  delete [] array1;
728  return returnCode;
729}
730// Creates direction vector
731void
732ClpSimplexNonlinear::directionVector (CoinIndexedVector * vectorArray,
733                                      CoinIndexedVector * spare1, CoinIndexedVector * spare2,
734                                      int pivotMode2,
735                                      double & normFlagged,double & normUnflagged,
736                                      int & numberNonBasic)
737{
738#if CLP_DEBUG > 1
739  vectorArray->checkClear();
740  spare1->checkClear();
741  spare2->checkClear();
742#endif
743  double *array = vectorArray->denseVector();
744  int * index = vectorArray->getIndices();
745  int number=0;
746  sequenceIn_=-1;
747  normFlagged=0.0;
748  normUnflagged=1.0;
749  double dualTolerance2 = CoinMin(1.0e-8,1.0e-2*dualTolerance_);
750  double dualTolerance3 = CoinMin(1.0e-3,1.0e3*dualTolerance_);
751  if (!numberNonBasic) {
752    //if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
753    //printf("infeasible\n");
754    if (!pivotMode2||pivotMode2>=10) {
755      normUnflagged=0.0;
756      double bestDj=0.0;
757      double bestSuper=0.0;
758      double sumSuper=0.0;
759      sequenceIn_=-1;
760      int nSuper=0;
761      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
762        array[iSequence]=0.0;
763        if (flagged(iSequence)) {
764          // accumulate  norm
765          switch(getStatus(iSequence)) {
766           
767          case basic:
768          case ClpSimplex::isFixed:
769            break;
770          case atUpperBound:
771            if (dj_[iSequence]>dualTolerance3) 
772              normFlagged += dj_[iSequence]*dj_[iSequence];
773            break;
774          case atLowerBound:
775            if (dj_[iSequence]<-dualTolerance3) 
776              normFlagged += dj_[iSequence]*dj_[iSequence];
777            break;
778          case isFree:
779          case superBasic:
780            if (fabs(dj_[iSequence])>dualTolerance3) 
781              normFlagged += dj_[iSequence]*dj_[iSequence];
782            break;
783          }
784          continue;
785        }
786        switch(getStatus(iSequence)) {
787         
788        case basic:
789        case ClpSimplex::isFixed:
790          break;
791        case atUpperBound:
792          if (dj_[iSequence]>dualTolerance_) {
793            if (dj_[iSequence]>dualTolerance3) 
794              normUnflagged += dj_[iSequence]*dj_[iSequence];
795            if (pivotMode2<10) {
796              array[iSequence]=-dj_[iSequence];
797              index[number++]=iSequence;
798            } else {
799              if (dj_[iSequence]>bestDj) {
800                bestDj=dj_[iSequence];
801                sequenceIn_=iSequence;
802              }
803            }
804          }
805          break;
806        case atLowerBound:
807          if (dj_[iSequence]<-dualTolerance_) {
808            if (dj_[iSequence]<-dualTolerance3) 
809              normUnflagged += dj_[iSequence]*dj_[iSequence];
810            if (pivotMode2<10) {
811              array[iSequence]=-dj_[iSequence];
812              index[number++]=iSequence;
813            } else {
814              if (-dj_[iSequence]>bestDj) {
815                bestDj=-dj_[iSequence];
816                sequenceIn_=iSequence;
817              }
818            }
819          }
820          break;
821        case isFree:
822        case superBasic:
823          //#define ALLSUPER
824#define NOSUPER
825#ifndef ALLSUPER
826          if (fabs(dj_[iSequence])>dualTolerance_) {
827            if (fabs(dj_[iSequence])>dualTolerance3) 
828              normUnflagged += dj_[iSequence]*dj_[iSequence];
829            nSuper++;
830            bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
831            sumSuper += fabs(dj_[iSequence]);
832          }
833          if (fabs(dj_[iSequence])>dualTolerance2) {
834            array[iSequence]=-dj_[iSequence];
835            index[number++]=iSequence;
836          }
837#else
838          array[iSequence]=-dj_[iSequence];
839          index[number++]=iSequence;
840          if (pivotMode2>=10)
841            bestSuper = CoinMax(fabs(dj_[iSequence]),bestSuper);
842#endif
843          break;
844        }
845      }
846#ifdef NOSUPER
847      // redo
848      bestSuper=sumSuper;
849      if(sequenceIn_>=0&&bestDj>bestSuper) {
850        int j;
851        // get rid of superbasics
852        for (j=0;j<number;j++) {
853          int iSequence = index[j];
854          array[iSequence]=0.0;
855        }
856        number=0;
857        array[sequenceIn_]=-dj_[sequenceIn_];
858        index[number++]=sequenceIn_;
859      } else {
860        sequenceIn_=-1;
861      }
862#else
863      if (pivotMode2>=10||!nSuper) {
864        bool takeBest=true;
865        if (pivotMode2==100&&nSuper>1)
866          takeBest=false;
867        if(sequenceIn_>=0&&takeBest) {
868          if (fabs(dj_[sequenceIn_])>bestSuper) {
869            array[sequenceIn_]=-dj_[sequenceIn_];
870            index[number++]=sequenceIn_;
871          } else {
872            sequenceIn_=-1;
873          }
874        } else {
875          sequenceIn_=-1;
876        }
877      }
878#endif
879#ifdef CLP_DEBUG
880      if (handler_->logLevel()&32) {
881        if (sequenceIn_>=0)
882          printf("%d superBasic, chosen %d - dj %g\n",nSuper,sequenceIn_,
883                 dj_[sequenceIn_]);
884        else
885          printf("%d superBasic - none chosen\n",nSuper);
886      }
887#endif
888    } else {
889      double bestDj=0.0;
890      double saveDj=0.0;
891      if (sequenceOut_>=0) {
892        saveDj=dj_[sequenceOut_];
893        dj_[sequenceOut_]=0.0;
894        switch(getStatus(sequenceOut_)) {
895         
896        case basic:
897          sequenceOut_=-1;
898        case ClpSimplex::isFixed:
899          break;
900        case atUpperBound:
901          if (dj_[sequenceOut_]>dualTolerance_) {
902#ifdef CLP_DEBUG
903            if (handler_->logLevel()&32) 
904              printf("after pivot out %d values %g %g %g, dj %g\n",
905                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
906                     upper_[sequenceOut_],dj_[sequenceOut_]);
907#endif
908          }
909          break;
910        case atLowerBound:
911          if (dj_[sequenceOut_]<-dualTolerance_) {
912#ifdef CLP_DEBUG
913            if (handler_->logLevel()&32) 
914              printf("after pivot out %d values %g %g %g, dj %g\n",
915                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
916                     upper_[sequenceOut_],dj_[sequenceOut_]);
917#endif
918          }
919          break;
920        case isFree:
921        case superBasic:
922          if (dj_[sequenceOut_]>dualTolerance_) {
923#ifdef CLP_DEBUG
924            if (handler_->logLevel()&32) 
925              printf("after pivot out %d values %g %g %g, dj %g\n",
926                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
927                     upper_[sequenceOut_],dj_[sequenceOut_]);
928#endif
929          } else if (dj_[sequenceOut_]<-dualTolerance_) {
930#ifdef CLP_DEBUG
931            if (handler_->logLevel()&32) 
932              printf("after pivot out %d values %g %g %g, dj %g\n",
933                     sequenceOut_,lower_[sequenceOut_],solution_[sequenceOut_],
934                     upper_[sequenceOut_],dj_[sequenceOut_]);
935#endif
936          }
937          break;
938        }
939      }
940      // Go for dj
941      pivotMode2=3;
942      for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
943        array[iSequence]=0.0;
944        if (flagged(iSequence))
945          continue;
946        switch(getStatus(iSequence)) {
947         
948        case basic:
949        case ClpSimplex::isFixed:
950          break;
951        case atUpperBound:
952          if (dj_[iSequence]>dualTolerance_) {
953            double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
954            double merit = distance*dj_[iSequence];
955            if (pivotMode2==1)
956              merit *= 1.0e-20; // discourage
957            if (pivotMode2==3)
958              merit = fabs(dj_[iSequence]);
959            if (merit>bestDj) {
960              sequenceIn_=iSequence;
961              bestDj=merit;
962            }
963          }
964          break;
965        case atLowerBound:
966          if (dj_[iSequence]<-dualTolerance_) {
967            double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
968            double merit = -distance*dj_[iSequence];
969            if (pivotMode2==1)
970              merit *= 1.0e-20; // discourage
971            if (pivotMode2==3)
972              merit = fabs(dj_[iSequence]);
973            if (merit>bestDj) {
974              sequenceIn_=iSequence;
975              bestDj=merit;
976            }
977          }
978          break;
979        case isFree:
980        case superBasic:
981          if (dj_[iSequence]>dualTolerance_) {
982            double distance = CoinMin(1.0e-2,solution_[iSequence]-lower_[iSequence]);
983            distance = CoinMin(solution_[iSequence]-lower_[iSequence],
984                           upper_[iSequence]-solution_[iSequence]);
985            double merit = distance*dj_[iSequence];
986            if (pivotMode2==1)
987              merit=distance;
988            if (pivotMode2==3)
989              merit = fabs(dj_[iSequence]);
990            if (merit>bestDj) {
991              sequenceIn_=iSequence;
992              bestDj=merit;
993            }
994          } else if (dj_[iSequence]<-dualTolerance_) {
995            double distance = CoinMin(1.0e-2,upper_[iSequence]-solution_[iSequence]);
996            distance = CoinMin(solution_[iSequence]-lower_[iSequence],
997                           upper_[iSequence]-solution_[iSequence]);
998            double merit = -distance*dj_[iSequence];
999            if (pivotMode2==1)
1000              merit=distance;
1001            if (pivotMode2==3)
1002              merit = fabs(dj_[iSequence]);
1003            if (merit>bestDj) {
1004              sequenceIn_=iSequence;
1005              bestDj=merit;
1006            }
1007          }
1008          break;
1009        }
1010      }
1011      if (sequenceOut_>=0) {
1012        dj_[sequenceOut_]=saveDj;
1013        sequenceOut_=-1;
1014      }
1015      if (sequenceIn_>=0) {
1016        array[sequenceIn_]=-dj_[sequenceIn_];
1017        index[number++]=sequenceIn_;
1018      }
1019    }
1020    numberNonBasic = number;
1021  } else {
1022    // compute norm of flagged
1023    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
1024      if (flagged(iSequence)) {
1025        // accumulate  norm
1026        switch(getStatus(iSequence)) {
1027         
1028        case basic:
1029        case ClpSimplex::isFixed:
1030          break;
1031        case atUpperBound:
1032          if (dj_[iSequence]>dualTolerance_) 
1033            normFlagged += dj_[iSequence]*dj_[iSequence];
1034          break;
1035        case atLowerBound:
1036          if (dj_[iSequence]<-dualTolerance_) 
1037            normFlagged += dj_[iSequence]*dj_[iSequence];
1038          break;
1039        case isFree:
1040        case superBasic:
1041          if (fabs(dj_[iSequence])>dualTolerance_) 
1042            normFlagged += dj_[iSequence]*dj_[iSequence];
1043          break;
1044        }
1045      }
1046    }
1047    // re-use list
1048    number=0;
1049    int j;
1050    for (j=0;j<numberNonBasic;j++) {
1051      int iSequence = index[j];
1052      if (flagged(iSequence))
1053        continue;
1054      switch(getStatus(iSequence)) {
1055         
1056      case basic:
1057      case ClpSimplex::isFixed:
1058        abort();
1059        break;
1060      case atUpperBound:
1061        if (dj_[iSequence]>dualTolerance_) 
1062          number++;
1063        break;
1064      case atLowerBound:
1065        if (dj_[iSequence]<-dualTolerance_) 
1066          number++;
1067        break;
1068      case isFree:
1069      case superBasic:
1070        if (fabs(dj_[iSequence])>dualTolerance_) 
1071          number++;
1072        break;
1073      }
1074      array[iSequence]=-dj_[iSequence];
1075    }
1076    if (!number) {
1077      for (j=0;j<numberNonBasic;j++) {
1078        int iSequence = index[j];
1079        array[iSequence]=0.0;
1080      }
1081      numberNonBasic=0;
1082    }
1083    number=numberNonBasic;
1084  }
1085  if (number) {
1086    int j;
1087    // Now do basic ones - how do I compensate for small basic infeasibilities?
1088    int iRow;
1089    for (iRow=0;iRow<numberRows_;iRow++) {
1090      int iPivot=pivotVariable_[iRow];
1091      double value=0.0;
1092      if (solution_[iPivot]>upper_[iPivot]) {
1093        value = upper_[iPivot]-solution_[iPivot];
1094      } else if (solution_[iPivot]<lower_[iPivot]) {
1095        value = lower_[iPivot]-solution_[iPivot];
1096      }
1097      //if (value)
1098      //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
1099      //     upper_[iPivot]);
1100      //value=0.0;
1101      value *= -1.0e0;
1102      if (value) {
1103        array[iPivot]=value;
1104        index[number++]=iPivot;
1105      }
1106    }
1107    double * array2=spare1->denseVector();
1108    int * index2 = spare1->getIndices();
1109    int number2=0;
1110    times(-1.0,array,array2);
1111    array = array+numberColumns_;
1112    for (iRow=0;iRow<numberRows_;iRow++) {
1113      double value = array2[iRow] + array[iRow];
1114      if (value) {
1115        array2[iRow]=value;
1116        index2[number2++]=iRow;
1117      } else {
1118        array2[iRow]=0.0;
1119      }
1120    }
1121    array -= numberColumns_;
1122    spare1->setNumElements(number2);
1123    // Ftran
1124    factorization_->updateColumn(spare2,spare1);
1125    number2=spare1->getNumElements();
1126    for (j=0;j<number2;j++) {
1127      int iSequence = index2[j];
1128      double value = array2[iSequence];
1129      array2[iSequence]=0.0;
1130      if (value) {
1131        int iPivot=pivotVariable_[iSequence];
1132        double oldValue = array[iPivot];
1133        if (!oldValue) {
1134          array[iPivot] = value;
1135          index[number++]=iPivot;
1136        } else {
1137          // something already there
1138          array[iPivot] = value+oldValue;
1139        }
1140      }
1141    }
1142    spare1->setNumElements(0);
1143  }
1144  vectorArray->setNumElements(number);
1145}
1146/*
1147   Row array and column array have direction
1148   Returns 0 - can do normal iteration (basis change)
1149   1 - no basis change
1150*/
1151int 
1152ClpSimplexNonlinear::pivotColumn(CoinIndexedVector * longArray,
1153                                 CoinIndexedVector * rowArray,
1154                                 CoinIndexedVector * columnArray,
1155                                 CoinIndexedVector * spare,
1156                                 int & pivotMode,
1157                                 double & solutionError,
1158                                 double * dArray)
1159{
1160  // say not optimal
1161  primalColumnPivot_->setLooksOptimal(false);
1162  double acceptablePivot=1.0e-10;
1163  int lastSequenceIn=-1;
1164  if (pivotMode&&pivotMode<10) {
1165    acceptablePivot=1.0e-6;
1166    if (factorization_->pivots())
1167      acceptablePivot=1.0e-5; // if we have iterated be more strict
1168  }
1169  double acceptableBasic=1.0e-7;
1170 
1171  int number=longArray->getNumElements();
1172  int numberTotal = numberRows_+numberColumns_;
1173  int bestSequence=-1;
1174  int bestBasicSequence=-1;
1175  double eps= 1.0e-1;
1176  eps=1.0e-6;
1177  double basicTheta = 1.0e30;
1178  double objTheta=0.0;
1179  bool finished = false;
1180  sequenceIn_=-1;
1181  int nPasses=0;
1182  int nTotalPasses=0;
1183  int nBigPasses=0;
1184  double djNorm0=0.0;
1185  double djNorm=0.0;
1186  double normFlagged=0.0;
1187  double normUnflagged=0.0;
1188  int localPivotMode=pivotMode;
1189  bool allFinished=false;
1190  bool justOne=false;
1191  int returnCode=1;
1192  double currentObj;
1193  double predictedObj;
1194  double thetaObj;
1195  objective_->stepLength(this,solution_,solution_,0.0,
1196                                               currentObj,predictedObj,thetaObj);
1197  double saveObj=currentObj;
1198#define MINTYPE 1
1199#if MINTYPE ==2
1200  // try Shanno's method
1201  // Copy code from ....cpp.q
1202#endif
1203  // big big loop
1204  while (!allFinished) {
1205    double * work=longArray->denseVector();
1206    int * which=longArray->getIndices();
1207    allFinished=true;
1208    // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
1209    //#define SMALLTHETA1 1.0e-25
1210    //#define SMALLTHETA2 1.0e-25
1211#define SMALLTHETA1 1.0e-10
1212#define SMALLTHETA2 1.0e-10
1213#define CONJUGATE 2
1214#if CONJUGATE == 0
1215    int conjugate = 0;
1216#elif CONJUGATE == 1
1217    int conjugate = (pivotMode<10) ? MINTYPE : 0;
1218#elif CONJUGATE == 2
1219    int conjugate = MINTYPE;
1220#else
1221    int conjugate = MINTYPE;
1222#endif
1223   
1224    //if (pivotMode==1)
1225    //localPivotMode=11;;
1226#if CLP_DEBUG > 1
1227    longArray->checkClear();
1228#endif
1229    int numberNonBasic=0;
1230    int startLocalMode=-1;
1231    while (!finished) {
1232      double simpleObjective=COIN_DBL_MAX;
1233      returnCode=1;
1234      int iSequence;
1235      objective_->reducedGradient(this,dj_,false);
1236      // get direction vector in longArray
1237      longArray->clear();
1238      // take out comment to try slightly different idea
1239      if (nPasses+nTotalPasses>3000)
1240        break;
1241      if (!nPasses) {
1242        numberNonBasic=0;
1243        nBigPasses++;
1244      }
1245      // just do superbasic if in cleanup mode
1246      int local=localPivotMode;
1247      if (!local&&pivotMode>=10&&nBigPasses<10) {
1248        local=100;
1249      } else if (local==-1||nBigPasses>=10) {
1250        local=0;
1251        localPivotMode=0;
1252      } 
1253      if (justOne) {
1254        local=2;
1255        //local=100;
1256        justOne=false;
1257      }
1258      if (!nPasses)
1259        startLocalMode=local;
1260      directionVector(longArray,spare,rowArray,local,
1261                      normFlagged,normUnflagged,numberNonBasic);
1262      {
1263        // check null vector
1264        double * rhs = spare->denseVector();
1265        int iRow;
1266        multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,rhs,0.0);
1267        matrix_->times(1.0,solution_,rhs,rowScale_,columnScale_);
1268        double largest=0.0;
1269        int iLargest=-1;
1270        for (iRow=0;iRow<numberRows_;iRow++) {
1271          double value=fabs(rhs[iRow]);
1272          rhs[iRow]=0.0;
1273          if (value>largest) {
1274            largest=value;
1275            iLargest=iRow;
1276          }
1277        }
1278#if CLP_DEBUG > 0
1279        if ((handler_->logLevel()&32)&&largest>1.0e-8) 
1280          printf("largest rhs error %g on row %d\n",largest,iLargest);
1281#endif
1282        if (solutionError<0.0) {
1283          solutionError=largest;
1284        } else if (largest>max(1.0e-8,1.0e2*solutionError)&&
1285                   factorization_->pivots()) {
1286          longArray->clear();
1287          pivotRow_ = -1;
1288          theta_=0.0;
1289          return 3;
1290        }
1291      }
1292      if (sequenceIn_>=0)
1293        lastSequenceIn=sequenceIn_;
1294      double djNormSave = djNorm;
1295      djNorm=0.0;
1296      int iIndex;
1297      for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
1298        int iSequence = which[iIndex];
1299        double alpha = work[iSequence];
1300        djNorm += alpha*alpha;
1301      }
1302      // go to conjugate gradient if necessary
1303      if (numberNonBasic&&localPivotMode>=10&&(nPasses>4||sequenceIn_<0)) {
1304        localPivotMode=0;
1305        nTotalPasses+= nPasses;
1306        nPasses=0;
1307      }
1308#if CONJUGATE == 2
1309      conjugate = (localPivotMode<10) ? MINTYPE : 0;
1310#endif
1311      if (!nPasses) {
1312        djNorm0 = CoinMax(djNorm,1.0e-20);
1313        memcpy(dArray,work,numberTotal*sizeof(double));
1314        if (sequenceIn_>=0&&numberNonBasic==1) {
1315          // see if simple move
1316          double objTheta2 =objective_->stepLength(this,solution_,work,1.0e30,
1317                                               currentObj,predictedObj,thetaObj);
1318          rowArray->clear();
1319         
1320          // update the incoming column
1321          unpackPacked(rowArray);
1322          factorization_->updateColumnFT(spare,rowArray);
1323          theta_=0.0;
1324          double * work2=rowArray->denseVector();
1325          int number=rowArray->getNumElements();
1326          int * which2=rowArray->getIndices();
1327          int iIndex;
1328          bool easyMove=false;
1329          double way;
1330          if (dj_[sequenceIn_]>0.0)
1331            way=-1.0;
1332          else
1333            way=1.0;
1334          double largest=COIN_DBL_MAX;
1335          int kPivot=-1;
1336          for (iIndex=0;iIndex<number;iIndex++) {
1337            int iRow = which2[iIndex];
1338            double alpha = way*work2[iIndex];
1339            int iPivot = pivotVariable_[iRow];
1340            if (alpha<-1.0e-5) {
1341              double distance = upper_[iPivot]-solution_[iPivot];
1342              if (distance<-largest*alpha) {
1343                kPivot=iPivot;
1344                largest=CoinMax(0.0,-distance/alpha);
1345              }
1346              if (distance<-1.0e-12*alpha) {
1347                easyMove=true;
1348                break;
1349              }
1350            } else if (alpha>1.0e-5) {
1351              double distance = solution_[iPivot]-lower_[iPivot];
1352              if (distance<largest*alpha) {
1353                kPivot=iPivot;
1354                largest=CoinMax(0.0,distance/alpha);
1355              }
1356              if (distance<1.0e-12*alpha) {
1357                easyMove=true;
1358                break;
1359              }
1360            }
1361          }
1362          if (largest<objTheta2) {
1363            easyMove=true;
1364          } else if (!easyMove) {
1365            double objDrop = currentObj-predictedObj;
1366            double th = objective_->stepLength(this,solution_,work,largest,
1367                                   currentObj,predictedObj,simpleObjective);
1368            simpleObjective = CoinMax(simpleObjective,predictedObj);
1369            double easyDrop = currentObj-simpleObjective;
1370            if (easyDrop>1.0e-8&&easyDrop>0.5*objDrop) {
1371              easyMove=true;
1372#ifdef CLP_DEBUG
1373              if (handler_->logLevel()&32) 
1374                printf("easy - obj drop %g, easy drop %g\n",objDrop,easyDrop);
1375#endif
1376              if (easyDrop>objDrop) {
1377                // debug
1378                printf("****** th %g simple %g\n",th,simpleObjective);
1379                objective_->stepLength(this,solution_,work,1.0e30,
1380                                   currentObj,predictedObj,simpleObjective);
1381                objective_->stepLength(this,solution_,work,largest,
1382                                   currentObj,predictedObj,simpleObjective);
1383              }
1384            }
1385          }
1386          rowArray->clear();
1387#ifdef CLP_DEBUG
1388          if (handler_->logLevel()&32)
1389            printf("largest %g piv %d\n",largest,kPivot);
1390#endif
1391          if (easyMove) {
1392            valueIn_=solution_[sequenceIn_];
1393            dualIn_=dj_[sequenceIn_];
1394            lowerIn_=lower_[sequenceIn_];
1395            upperIn_=upper_[sequenceIn_];
1396            if (dualIn_>0.0)
1397              directionIn_ = -1;
1398            else 
1399              directionIn_ = 1;
1400            longArray->clear();
1401            pivotRow_ = -1;
1402            theta_=0.0;
1403            return 0;
1404          }
1405        }
1406      } else {
1407        double beta = djNorm/djNormSave;
1408        if (conjugate) {
1409          //printf("current norm %g, old %g - beta %g\n",
1410          //     djNorm,djNormSave,beta);
1411          for (iSequence=0;iSequence<numberTotal;iSequence++)
1412            dArray[iSequence] = work[iSequence] + beta *dArray[iSequence];
1413        } else {
1414          for (iSequence=0;iSequence<numberTotal;iSequence++)
1415            dArray[iSequence] = work[iSequence];
1416        }
1417      }
1418      if (djNorm<eps*djNorm0||(nPasses>100&&djNorm<CoinMin(1.0e-1*djNorm0,1.0e-12))) {
1419#ifdef CLP_DEBUG
1420        if (handler_->logLevel()&32) 
1421          printf("dj norm reduced from %g to %g\n",djNorm0,djNorm);
1422#endif
1423        if (pivotMode<10||!numberNonBasic) {
1424          finished=true;
1425        } else {
1426          localPivotMode=pivotMode;
1427          nTotalPasses+= nPasses;
1428          nPasses=0;
1429          continue;
1430        }
1431      }
1432      //if (nPasses==100||nPasses==50)
1433      //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
1434      //         normFlagged);
1435      if (nPasses>100&&djNorm<1.0e-2*normFlagged&&!startLocalMode) {
1436#ifdef CLP_DEBUG
1437        if (handler_->logLevel()&32) 
1438          printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n",djNorm0,djNorm,
1439                 normFlagged);
1440#endif
1441        unflag();
1442        localPivotMode=0;
1443        nTotalPasses+= nPasses;
1444        nPasses=0;
1445        continue;
1446      }
1447      if (djNorm>0.99*djNorm0&&nPasses>1500) {
1448        finished=true;
1449#ifdef CLP_DEBUG
1450        if (handler_->logLevel()&32) 
1451          printf("dj norm NOT reduced from %g to %g\n",djNorm0,djNorm);
1452#endif
1453        djNorm=1.2345e-20;
1454      }
1455      number=longArray->getNumElements();
1456      if (!numberNonBasic) {
1457        // looks optimal
1458        // check if any dj look plausible
1459        int nSuper=0;
1460        int nFlagged=0;
1461        int nNormal=0;
1462        for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
1463          if (flagged(iSequence)) {
1464            switch(getStatus(iSequence)) {
1465             
1466            case basic:
1467            case ClpSimplex::isFixed:
1468              break;
1469            case atUpperBound:
1470              if (dj_[iSequence]>dualTolerance_) 
1471                nFlagged++;
1472              break;
1473            case atLowerBound:
1474              if (dj_[iSequence]<-dualTolerance_) 
1475                nFlagged++;
1476              break;
1477            case isFree:
1478            case superBasic:
1479              if (fabs(dj_[iSequence])>dualTolerance_) 
1480                nFlagged++;
1481              break;
1482            }
1483            continue;
1484          }
1485          switch(getStatus(iSequence)) {
1486           
1487          case basic:
1488          case ClpSimplex::isFixed:
1489            break;
1490          case atUpperBound:
1491            if (dj_[iSequence]>dualTolerance_) 
1492              nNormal++;
1493            break;
1494          case atLowerBound:
1495            if (dj_[iSequence]<-dualTolerance_) 
1496              nNormal++;
1497            break;
1498          case isFree:
1499          case superBasic:
1500            if (fabs(dj_[iSequence])>dualTolerance_) 
1501              nSuper++;
1502            break;
1503          }
1504        }
1505#ifdef CLP_DEBUG
1506        if (handler_->logLevel()&32) 
1507          printf("%d super, %d normal, %d flagged\n",
1508                 nSuper,nNormal,nFlagged);
1509#endif
1510       
1511        int nFlagged2=1;
1512        if (lastSequenceIn<0&&!nNormal&&!nSuper) {
1513          nFlagged2=unflag();
1514          if (pivotMode>=10) {
1515            pivotMode--;
1516#ifdef CLP_DEBUG
1517            if (handler_->logLevel()&32) 
1518              printf("pivot mode now %d\n",pivotMode);
1519#endif
1520            if (pivotMode==9) 
1521              pivotMode=0;      // switch off fast attempt
1522          }
1523        } else {
1524          lastSequenceIn=-1;
1525        }
1526        if (!nFlagged2||(normFlagged+normUnflagged<1.0e-8)) {
1527          primalColumnPivot_->setLooksOptimal(true);
1528          return 0;
1529        } else {
1530          localPivotMode=-1;
1531          nTotalPasses+= nPasses;
1532          nPasses=0;
1533          finished=false;
1534          continue;
1535        }
1536      }
1537      bestSequence=-1;
1538      bestBasicSequence=-1;
1539     
1540      // temp
1541      nPasses++;
1542      if (nPasses>2000)
1543        finished=true;
1544      double theta = 1.0e30;
1545      basicTheta = 1.0e30;
1546      theta_=1.0e30;
1547      double basicTolerance = 1.0e-4*primalTolerance_;
1548      for (iSequence=0;iSequence<numberTotal;iSequence++) {
1549        //if (flagged(iSequence)
1550        //  continue;
1551        double alpha = dArray[iSequence];
1552        Status thisStatus = getStatus(iSequence);
1553        double oldValue = solution_[iSequence];
1554        if (thisStatus!=basic) {
1555          if (fabs(alpha)>=acceptablePivot) {
1556            if (alpha<0.0) {
1557              // variable going towards lower bound
1558              double bound = lower_[iSequence];
1559              oldValue -= bound;
1560              if (oldValue+theta*alpha<0.0) {
1561                bestSequence=iSequence;
1562                theta = CoinMax(0.0,oldValue/(-alpha));
1563              }
1564            } else {
1565              // variable going towards upper bound
1566              double bound = upper_[iSequence];
1567              oldValue = bound-oldValue;
1568              if (oldValue-theta*alpha<0.0) {
1569                bestSequence=iSequence;
1570                theta = CoinMax(0.0,oldValue/alpha);
1571              }
1572            }
1573          }
1574        } else {
1575          if (fabs(alpha)>=acceptableBasic) {
1576            if (alpha<0.0) {
1577              // variable going towards lower bound
1578              double bound = lower_[iSequence];
1579              oldValue -= bound;
1580              if (oldValue+basicTheta*alpha<-basicTolerance) {
1581                bestBasicSequence=iSequence;
1582                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/(-alpha));
1583              }
1584            } else {
1585              // variable going towards upper bound
1586              double bound = upper_[iSequence];
1587              oldValue = bound-oldValue;
1588              if (oldValue-basicTheta*alpha<-basicTolerance) {
1589                bestBasicSequence=iSequence;
1590                basicTheta = CoinMax(0.0,(oldValue+basicTolerance)/alpha);
1591              }
1592            }
1593          }
1594        }
1595      }
1596      theta_ = CoinMin(theta,basicTheta);
1597      // Now find minimum of function
1598      double objTheta2 =objective_->stepLength(this,solution_,dArray,CoinMin(theta,basicTheta),
1599                                               currentObj,predictedObj,thetaObj);
1600#ifdef CLP_DEBUG
1601      if (handler_->logLevel()&32) 
1602        printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
1603#endif
1604      if (conjugate) {
1605        double offset;
1606        const double * gradient = objective_->gradient(this,
1607                                                       dArray, offset,
1608                                                       true,0);
1609        double product=0.0;
1610        for (iSequence = 0;iSequence<numberColumns_;iSequence++) {
1611          double alpha = dArray[iSequence];
1612          double value = alpha*gradient[iSequence];
1613          product += value;
1614        }
1615        //#define INCLUDESLACK
1616#ifdef INCLUDESLACK
1617        for (;iSequence<numberColumns_+numberRows_;iSequence++) {
1618          double alpha = dArray[iSequence];
1619          double value = alpha*cost_[iSequence];
1620          product += value;
1621        }
1622#endif
1623        if (product>0.0)
1624          objTheta = djNorm/product;
1625        else
1626          objTheta=COIN_DBL_MAX;
1627        assert (product>-1.0e-5);
1628      } else {
1629        objTheta =objTheta2;
1630      }
1631      // if very small difference then take pivot (depends on djNorm?)
1632      // distinguish between basic and non-basic
1633      bool chooseObjTheta=objTheta<theta_;
1634      if (chooseObjTheta) {
1635        if (thetaObj<currentObj-1.0e-12&&objTheta+1.0e-10>theta_)
1636          chooseObjTheta=false;
1637        //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
1638        //chooseObjTheta=false;
1639      }
1640      //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) {
1641      if (chooseObjTheta) {
1642        theta_ = objTheta;
1643      } else {
1644        objTheta = CoinMax(objTheta,1.00000001*theta_+1.0e-12);
1645        //if (theta+1.0e-13>basicTheta) {
1646        //theta = CoinMax(theta,1.00000001*basicTheta);
1647        //theta_ = basicTheta;
1648        //}
1649      }
1650      // always out if one variable in and zero move
1651      if (theta_==basicTheta||(sequenceIn_>=0&&!theta_))
1652        finished=true; // come out
1653#ifdef CLP_DEBUG
1654      if (handler_->logLevel()&32) {
1655        printf("%d pass,",nPasses);
1656        if (sequenceIn_>=0)
1657          printf (" Sin %d,",sequenceIn_);
1658        if (basicTheta==theta_)
1659          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1660        else
1661          printf(" basicTheta %g",basicTheta);
1662        if (theta==theta_)
1663          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1664        else
1665          printf(" non-basicTheta %g",theta);
1666        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1667        printf(" djNorm %g\n",djNorm);
1668      }
1669#endif
1670      if (handler_->logLevel()>3&&objTheta!=theta_) {
1671        printf("%d pass obj %g,",nPasses,currentObj);
1672        if (sequenceIn_>=0)
1673          printf (" Sin %d,",sequenceIn_);
1674        if (basicTheta==theta_)
1675          printf(" X(%d) basicTheta %g",bestBasicSequence,basicTheta);
1676        else
1677          printf(" basicTheta %g",basicTheta);
1678        if (theta==theta_)
1679          printf(" X(%d) non-basicTheta %g",bestSequence,theta);
1680        else
1681          printf(" non-basicTheta %g",theta);
1682        printf(" %s objTheta %g",objTheta==theta_ ? "X" : "",objTheta);
1683        printf(" djNorm %g\n",djNorm);
1684      }
1685      if (objTheta!=theta_) {
1686        //printf("hit boundary after %d passes\n",nPasses);
1687        nTotalPasses+= nPasses;
1688        nPasses=0; // start again
1689      }
1690      if (localPivotMode<10||lastSequenceIn==bestSequence) {
1691        if (theta_==theta&&theta_<basicTheta&&theta_<1.0e-5)
1692          setFlagged(bestSequence); // out of active set
1693      }
1694      // Update solution
1695      for (iSequence=0;iSequence<numberTotal;iSequence++) {
1696        //for (iIndex=0;iIndex<number;iIndex++) {
1697       
1698        //int iSequence = which[iIndex];
1699        double alpha = dArray[iSequence];
1700        if (alpha) {
1701          double value = solution_[iSequence]+theta_*alpha;
1702          solution_[iSequence]=value;
1703          switch(getStatus(iSequence)) {
1704           
1705          case basic:
1706          case isFixed:
1707          case isFree:
1708          case atUpperBound:
1709          case atLowerBound:
1710            nonLinearCost_->setOne(iSequence,value);
1711            break;
1712          case superBasic:
1713            // To get correct action
1714            setStatus(iSequence,isFixed);
1715            nonLinearCost_->setOne(iSequence,value);
1716            //assert (getStatus(iSequence)!=isFixed);
1717            break;
1718          }
1719        }
1720      }
1721      if (objTheta<SMALLTHETA2&&objTheta==theta_) {
1722        if (sequenceIn_>=0&&basicTheta<1.0e-9) {
1723          // try just one
1724          localPivotMode=0;
1725          sequenceIn_=-1;
1726          nTotalPasses+= nPasses;
1727          nPasses=0;
1728          //finished=true;
1729          //objTheta=0.0;
1730          //theta_=0.0;
1731        }
1732        //if (objTheta<1.0e-12)
1733        //finished=true;
1734        //printf("zero move\n");
1735      }
1736#ifdef CLP_DEBUG
1737      if (handler_->logLevel()&32) {
1738        objective_->stepLength(this,solution_,work,0.0,
1739                               currentObj,predictedObj,thetaObj);
1740        printf("current obj %g after update - simple was %g\n",currentObj,simpleObjective);
1741      }
1742#endif
1743      if (sequenceIn_>=0&&!finished&&objTheta>1.0e-4) {
1744        // we made some progress - back to normal
1745        if (localPivotMode<10) {
1746          localPivotMode=0;
1747          sequenceIn_=-1;
1748          nTotalPasses+= nPasses;
1749          nPasses=0;
1750        }
1751#ifdef CLP_DEBUG
1752        if (handler_->logLevel()&32) 
1753          printf("some progress?\n");
1754#endif
1755      }
1756#if CLP_DEBUG > 1
1757      longArray->checkClean();
1758#endif
1759    }
1760#ifdef CLP_DEBUG
1761    if (handler_->logLevel()&32) 
1762      printf("out of loop after %d passes\n",nPasses);
1763#endif
1764    bool ordinaryDj=false;
1765    //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
1766    //printf("could try ordinary iteration %g\n",theta_);
1767    if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-15) {
1768      //printf("trying ordinary iteration\n");
1769      ordinaryDj=true;
1770    }
1771    if (!basicTheta&&!ordinaryDj) {
1772      //returnCode=2;
1773      //objTheta=-1.0; // so we fall through
1774    }
1775    assert (theta_<1.0e30); // for now
1776    // See if we need to pivot
1777    if (theta_==basicTheta||ordinaryDj) {
1778      if (!ordinaryDj) {
1779        // find an incoming column which will force pivot
1780        int iRow;
1781        pivotRow_=-1;
1782        for (iRow=0;iRow<numberRows_;iRow++) {
1783          if (pivotVariable_[iRow]==bestBasicSequence) {
1784            pivotRow_=iRow;
1785            break;
1786          }
1787        }
1788        assert (pivotRow_>=0);
1789        // Get good size for pivot
1790        double acceptablePivot=1.0e-7;
1791        if (factorization_->pivots()>10)
1792          acceptablePivot=1.0e-5; // if we have iterated be more strict
1793        else if (factorization_->pivots()>5)
1794          acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
1795        // should be dArray but seems better this way!
1796        double direction = work[bestBasicSequence]>0.0 ? -1.0 : 1.0;
1797        // create as packed
1798        rowArray->createPacked(1,&pivotRow_,&direction);
1799        factorization_->updateColumnTranspose(spare,rowArray);
1800        // put row of tableau in rowArray and columnArray
1801        matrix_->transposeTimes(this,-1.0,
1802                                rowArray,spare,columnArray);
1803        // choose one futhest away from bound which has reasonable pivot
1804        // If increasing we want negative alpha
1805       
1806        double * work2;
1807        int iSection;
1808       
1809        sequenceIn_=-1;
1810        double bestValue=-1.0;
1811        double bestDirection=0.0;
1812        // First pass we take correct direction and large pivots
1813        // then correct direction
1814        // then any
1815        double check[]={1.0e-8,-1.0e-12,-1.0e30};
1816        double mult[]={100.0,1.0,1.0};
1817        for (int iPass=0;iPass<3;iPass++) {
1818          //if (!bestValue&&iPass==2)
1819          //bestValue=-1.0;
1820          double acceptable=acceptablePivot*mult[iPass];
1821          double checkValue=check[iPass];
1822          for (iSection=0;iSection<2;iSection++) {
1823           
1824            int addSequence;
1825           
1826            if (!iSection) {
1827              work2 = rowArray->denseVector();
1828              number = rowArray->getNumElements();
1829              which = rowArray->getIndices();
1830              addSequence = numberColumns_;
1831            } else {
1832              work2 = columnArray->denseVector();
1833              number = columnArray->getNumElements();
1834              which = columnArray->getIndices();
1835              addSequence = 0;
1836            }
1837            int i;
1838           
1839            for (i=0;i<number;i++) {
1840              int iSequence = which[i]+addSequence;
1841              if (flagged(iSequence))
1842                continue;
1843              //double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
1844              //                  upper_[iSequence]-solution_[iSequence]);
1845              double alpha=work2[i];
1846              // should be dArray but seems better this way!
1847              double change=work[iSequence];
1848              Status thisStatus = getStatus(iSequence);
1849              double direction=0;;
1850              switch(thisStatus) {
1851               
1852              case basic:
1853              case ClpSimplex::isFixed:
1854                break;
1855              case isFree:
1856              case superBasic:
1857                if (alpha<-acceptable&&change>checkValue)
1858                  direction=1.0;
1859                else if (alpha>acceptable&&change<-checkValue)
1860                  direction=-1.0;
1861                break;
1862              case atUpperBound:
1863                if (alpha>acceptable&&change<-checkValue)
1864                  direction=-1.0;
1865                else if (iPass==2&&alpha<-acceptable&&change<-checkValue)
1866                  direction=1.0;
1867                break;
1868              case atLowerBound:
1869                if (alpha<-acceptable&&change>checkValue)
1870                  direction=1.0;
1871                else if (iPass==2&&alpha>acceptable&&change>checkValue)
1872                  direction=-1.0;
1873                break;
1874              }
1875              if (direction) {
1876                if (sequenceIn_!=lastSequenceIn||localPivotMode<10) { 
1877                  if (CoinMin(solution_[iSequence]-lower_[iSequence],
1878                          upper_[iSequence]-solution_[iSequence])>bestValue) {
1879                    bestValue=CoinMin(solution_[iSequence]-lower_[iSequence],
1880                                  upper_[iSequence]-solution_[iSequence]);
1881                    sequenceIn_=iSequence;
1882                    bestDirection=direction;
1883                  }
1884                } else {
1885                  // choose
1886                  bestValue=COIN_DBL_MAX;
1887                  sequenceIn_=iSequence;
1888                  bestDirection=direction;
1889                }
1890              }
1891            }
1892          }
1893          if (sequenceIn_>=0&&bestValue>0.0)
1894            break;
1895        }
1896        if (sequenceIn_>=0) {
1897          valueIn_=solution_[sequenceIn_];
1898          dualIn_=dj_[sequenceIn_];
1899          if (bestDirection<0.0) {
1900            // we want positive dj
1901            if (dualIn_<=0.0) {
1902              //printf("bad dj - xx %g\n",dualIn_);
1903              dualIn_=1.0;
1904            }
1905          } else {
1906            // we want negative dj
1907            if (dualIn_>=0.0) {
1908              //printf("bad dj - xx %g\n",dualIn_);
1909              dualIn_=-1.0;
1910            }
1911          }
1912          lowerIn_=lower_[sequenceIn_];
1913          upperIn_=upper_[sequenceIn_];
1914          if (dualIn_>0.0)
1915            directionIn_ = -1;
1916          else 
1917            directionIn_ = 1;
1918        } else {
1919          //ordinaryDj=true;
1920#ifdef CLP_DEBUG
1921          if (handler_->logLevel()&32) {
1922            printf("no easy pivot - norm %g mode %d\n",djNorm,localPivotMode);
1923            if (rowArray->getNumElements()+columnArray->getNumElements()<12) {
1924              for (iSection=0;iSection<2;iSection++) {
1925               
1926                int addSequence;
1927               
1928                if (!iSection) {
1929                  work2 = rowArray->denseVector();
1930                  number = rowArray->getNumElements();
1931                  which = rowArray->getIndices();
1932                  addSequence = numberColumns_;
1933                } else {
1934                  work2 = columnArray->denseVector();
1935                  number = columnArray->getNumElements();
1936                  which = columnArray->getIndices();
1937                  addSequence = 0;
1938                }
1939                int i;
1940                char section[]={'R','C'};
1941                for (i=0;i<number;i++) {
1942                  int iSequence = which[i]+addSequence;
1943                  if (flagged(iSequence)) {
1944                    printf("%c%d flagged\n",section[iSection],which[i]);
1945                    continue;
1946                  } else {
1947                    printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
1948                           section[iSection],which[i],status_[iSequence],
1949                           lower_[iSequence],solution_[iSequence],upper_[iSequence],
1950                           work2[i],work[iSequence]);
1951                  }
1952                }
1953              }
1954            }
1955          }
1956#endif
1957          assert (sequenceIn_<0);
1958          justOne=true;
1959          allFinished=false; // Round again
1960          finished=false;
1961          nTotalPasses+= nPasses;
1962          nPasses=0;
1963          if (djNorm<0.9*djNorm0&&djNorm<1.0e-3&&!localPivotMode) {
1964#ifdef CLP_DEBUG
1965            if (handler_->logLevel()&32) 
1966              printf("no pivot - mode %d norms %g %g - unflagging\n",
1967                     localPivotMode,djNorm0,djNorm);
1968#endif
1969            unflag(); //unflagging
1970            returnCode=1;
1971          } else {
1972            returnCode=2; // do single incoming
1973            //returnCode=1;
1974          }
1975        }
1976      }
1977      rowArray->clear();
1978      columnArray->clear();
1979      longArray->clear();
1980      if (ordinaryDj) {
1981        valueIn_=solution_[sequenceIn_];
1982        dualIn_=dj_[sequenceIn_];
1983        lowerIn_=lower_[sequenceIn_];
1984        upperIn_=upper_[sequenceIn_];
1985        if (dualIn_>0.0)
1986          directionIn_ = -1;
1987        else 
1988          directionIn_ = 1;
1989      }
1990      if (returnCode==1)
1991        returnCode=0;
1992    } else {
1993      // round again
1994      longArray->clear();
1995      if (djNorm<1.0e-3&&!localPivotMode) {
1996        if (djNorm==1.2345e-20&&djNorm0>1.0e-4) {
1997#ifdef CLP_DEBUG
1998          if (handler_->logLevel()&32) 
1999            printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n",djNorm0,nPasses,
2000                   localPivotMode,returnCode);
2001#endif
2002          //if (!localPivotMode)
2003          //returnCode=2; // force singleton
2004        } else {
2005#ifdef CLP_DEBUG
2006          if (handler_->logLevel()&32) 
2007            printf("unflagging as djNorm %g %g, %d passes\n",djNorm,djNorm0,nPasses);
2008#endif
2009          if (pivotMode>=10) {
2010            pivotMode--;
2011#ifdef CLP_DEBUG
2012            if (handler_->logLevel()&32) 
2013              printf("pivot mode now %d\n",pivotMode);
2014#endif
2015            if (pivotMode==9) 
2016              pivotMode=0;      // switch off fast attempt
2017          }
2018          bool unflagged = unflag()!=0;
2019          if (!unflagged&&djNorm<1.0e-10) {
2020            // ? declare victory
2021            sequenceIn_=-1;
2022            returnCode=0;
2023          }
2024        }
2025      }
2026    }
2027  }
2028  if (djNorm0<1.0e-12*normFlagged) {
2029#ifdef CLP_DEBUG
2030    if (handler_->logLevel()&32) 
2031      printf("unflagging as djNorm %g %g and flagged norm %g\n",djNorm,djNorm0,normFlagged);
2032#endif
2033    unflag();
2034  }
2035  if (saveObj-currentObj<1.0e-5&&nTotalPasses>2000) {
2036    normUnflagged=0.0;
2037    double dualTolerance3 = CoinMin(1.0e-3,1.0e3*dualTolerance_);
2038    for (int iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
2039      switch(getStatus(iSequence)) {
2040       
2041      case basic:
2042      case ClpSimplex::isFixed:
2043        break;
2044      case atUpperBound:
2045        if (dj_[iSequence]>dualTolerance3) 
2046          normFlagged += dj_[iSequence]*dj_[iSequence];
2047        break;
2048      case atLowerBound:
2049        if (dj_[iSequence]<-dualTolerance3) 
2050          normFlagged += dj_[iSequence]*dj_[iSequence];
2051        break;
2052      case isFree:
2053      case superBasic:
2054        if (fabs(dj_[iSequence])>dualTolerance3) 
2055          normFlagged += dj_[iSequence]*dj_[iSequence];
2056        break;
2057      }
2058    }
2059    if (handler_->logLevel()>2)
2060      printf("possible optimal  %d %d %g %g\n",
2061             nBigPasses,nTotalPasses,saveObj-currentObj,normFlagged);
2062    if (normFlagged<1.0e-5) {
2063      unflag();
2064      primalColumnPivot_->setLooksOptimal(true);
2065      returnCode=0;
2066    }
2067  }
2068  return returnCode;
2069}
2070/* Do last half of an iteration.
2071   Return codes
2072   Reasons to come out normal mode
2073   -1 normal
2074   -2 factorize now - good iteration
2075   -3 slight inaccuracy - refactorize - iteration done
2076   -4 inaccuracy - refactorize - no iteration
2077   -5 something flagged - go round again
2078   +2 looks unbounded
2079   +3 max iterations (iteration done)
2080   
2081*/
2082int 
2083ClpSimplexNonlinear::pivotNonlinearResult()
2084{
2085
2086  int returnCode=-1;
2087
2088  rowArray_[1]->clear();
2089 
2090  // we found a pivot column
2091  // update the incoming column
2092  unpackPacked(rowArray_[1]);
2093  factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2094  theta_=0.0;
2095  double * work=rowArray_[1]->denseVector();
2096  int number=rowArray_[1]->getNumElements();
2097  int * which=rowArray_[1]->getIndices();
2098  bool keepValue=false;
2099  double saveValue=0.0;
2100  if (pivotRow_>=0) {
2101    sequenceOut_=pivotVariable_[pivotRow_];;
2102    valueOut_ = solution(sequenceOut_);
2103    keepValue=true;
2104    saveValue=valueOut_;
2105    lowerOut_=lower_[sequenceOut_];
2106    upperOut_=upper_[sequenceOut_];
2107    int iIndex;
2108    for (iIndex=0;iIndex<number;iIndex++) {
2109     
2110      int iRow = which[iIndex];
2111      if (iRow==pivotRow_) {
2112        alpha_ = work[iIndex];
2113        break;
2114      }
2115    }
2116  } else {
2117    int iIndex;
2118    double smallest=COIN_DBL_MAX;
2119    for (iIndex=0;iIndex<number;iIndex++) {
2120      int iRow = which[iIndex];
2121      double alpha = work[iIndex];
2122      if (fabs(alpha)>1.0e-6) {
2123        int iPivot = pivotVariable_[iRow];
2124        double distance = CoinMin(upper_[iPivot]-solution_[iPivot],
2125                              solution_[iPivot]-lower_[iPivot]);
2126        if (distance<smallest) {
2127          pivotRow_=iRow;
2128          alpha_=alpha;
2129          smallest=distance;
2130        }
2131      }
2132    }
2133    if (smallest>primalTolerance_) {
2134      smallest=COIN_DBL_MAX;
2135      for (iIndex=0;iIndex<number;iIndex++) {
2136        int iRow = which[iIndex];
2137        double alpha = work[iIndex];
2138        if (fabs(alpha)>1.0e-6) {
2139          double distance = CoinDrand48();
2140          if (distance<smallest) {
2141            pivotRow_=iRow;
2142            alpha_=alpha;
2143            smallest=distance;
2144          }
2145        }
2146      }
2147    }
2148    assert (pivotRow_>=0);
2149    sequenceOut_=pivotVariable_[pivotRow_];;
2150    valueOut_ = solution(sequenceOut_);
2151    lowerOut_=lower_[sequenceOut_];
2152    upperOut_=upper_[sequenceOut_];
2153  }
2154  double newValue = valueOut_ - theta_*alpha_;
2155  bool isSuperBasic=false;
2156  if (valueOut_>=upperOut_-primalTolerance_) {
2157    directionOut_=-1;      // to upper bound
2158    upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2159    upperOut_ = newValue;
2160  } else if (valueOut_<=lowerOut_+primalTolerance_) {
2161    directionOut_=1;      // to lower bound
2162    lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2163  } else {
2164    lowerOut_=valueOut_;
2165    upperOut_=valueOut_;
2166    isSuperBasic=true;
2167    //printf("XX superbasic out\n");
2168  }
2169  dualOut_ = dj_[sequenceOut_];
2170  double checkValue=1.0e-2;
2171  if (largestDualError_>1.0e-5)
2172    checkValue=1.0e-1;
2173  // if stable replace in basis
2174 
2175  int updateStatus = factorization_->replaceColumn(this,
2176                                                   rowArray_[2],
2177                                                   rowArray_[1],
2178                                                   pivotRow_,
2179                                                   alpha_);
2180 
2181  // if no pivots, bad update but reasonable alpha - take and invert
2182  if (updateStatus==2&&
2183      lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2184    updateStatus=4;
2185  if (updateStatus==1||updateStatus==4) {
2186    // slight error
2187    if (factorization_->pivots()>5||updateStatus==4) {
2188      returnCode=-3;
2189    }
2190  } else if (updateStatus==2) {
2191    // major error
2192    // better to have small tolerance even if slower
2193    factorization_->zeroTolerance(1.0e-15);
2194    int maxFactor = factorization_->maximumPivots();
2195    if (maxFactor>10) {
2196      if (forceFactorization_<0)
2197        forceFactorization_= maxFactor;
2198      forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2199    } 
2200    // later we may need to unwind more e.g. fake bounds
2201    if(lastGoodIteration_ != numberIterations_) {
2202      clearAll();
2203      pivotRow_=-1;
2204      returnCode=-4;
2205    } else {
2206      // need to reject something
2207      char x = isColumn(sequenceIn_) ? 'C' :'R';
2208      handler_->message(CLP_SIMPLEX_FLAG,messages_)
2209        <<x<<sequenceWithin(sequenceIn_)
2210        <<CoinMessageEol;
2211      setFlagged(sequenceIn_);
2212      progress_->clearBadTimes();
2213      lastBadIteration_ = numberIterations_; // say be more cautious
2214      clearAll();
2215      pivotRow_=-1;
2216      sequenceOut_=-1;
2217      returnCode = -5;
2218     
2219    }
2220    return returnCode;
2221  } else if (updateStatus==3) {
2222    // out of memory
2223    // increase space if not many iterations
2224    if (factorization_->pivots()<
2225        0.5*factorization_->maximumPivots()&&
2226        factorization_->pivots()<200)
2227      factorization_->areaFactor(
2228                                 factorization_->areaFactor() * 1.1);
2229    returnCode =-2; // factorize now
2230  } else if (updateStatus==5) {
2231    problemStatus_=-2; // factorize now
2232  }
2233 
2234  // update primal solution
2235 
2236  double objectiveChange=0.0;
2237  // after this rowArray_[1] is not empty - used to update djs
2238  // If pivot row >= numberRows then may be gub
2239  updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1);
2240 
2241  double oldValue = valueIn_;
2242  if (directionIn_==-1) {
2243    // as if from upper bound
2244    if (sequenceIn_!=sequenceOut_) {
2245      // variable becoming basic
2246      valueIn_ -= fabs(theta_);
2247    } else {
2248      valueIn_=lowerIn_;
2249    }
2250  } else {
2251    // as if from lower bound
2252    if (sequenceIn_!=sequenceOut_) {
2253      // variable becoming basic
2254      valueIn_ += fabs(theta_);
2255    } else {
2256      valueIn_=upperIn_;
2257    }
2258  }
2259  objectiveChange += dualIn_*(valueIn_-oldValue);
2260  // outgoing
2261  if (sequenceIn_!=sequenceOut_) {
2262    if (directionOut_>0) {
2263      valueOut_ = lowerOut_;
2264    } else {
2265      valueOut_ = upperOut_;
2266    }
2267    if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2268      valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2269    else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2270      valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2271    // may not be exactly at bound and bounds may have changed
2272    // Make sure outgoing looks feasible
2273    if (!isSuperBasic)
2274      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2275    solution_[sequenceOut_]=valueOut_;
2276  }
2277  // change cost and bounds on incoming if primal
2278  nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2279  int whatNext=housekeeping(objectiveChange);
2280  if (keepValue)
2281    solution_[sequenceOut_]=saveValue;
2282  if (isSuperBasic)
2283    setStatus(sequenceOut_,superBasic);
2284  //#define CLP_DEBUG
2285#if CLP_DEBUG > 1
2286  {
2287    int ninf= matrix_->checkFeasible(this);
2288    if (ninf)
2289      printf("infeas %d\n",ninf);
2290  }
2291#endif
2292  if (whatNext==1) {
2293    returnCode =-2; // refactorize
2294  } else if (whatNext==2) {
2295    // maximum iterations or equivalent
2296    returnCode=3;
2297  } else if(numberIterations_ == lastGoodIteration_
2298            + 2 * factorization_->maximumPivots()) {
2299    // done a lot of flips - be safe
2300    returnCode =-2; // refactorize
2301  }
2302  // Check event
2303  {
2304    int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2305    if (status>=0) {
2306      problemStatus_=5;
2307      secondaryStatus_=ClpEventHandler::endOfIteration;
2308      returnCode=4;
2309    }
2310  }
2311  return returnCode;
2312}
2313// A sequential LP method
2314int 
2315ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
2316{
2317  // Are we minimizing or maximizing
2318  double whichWay=optimizationDirection();
2319  if (whichWay<0.0)
2320    whichWay=-1.0;
2321  else if (whichWay>0.0)
2322    whichWay=1.0;
2323
2324
2325  int numberColumns = this->numberColumns();
2326  int numberRows = this->numberRows();
2327  double * columnLower = this->columnLower();
2328  double * columnUpper = this->columnUpper();
2329  double * solution = this->primalColumnSolution();
2330 
2331  if (objective_->type()<2) {
2332    // no nonlinear part
2333    return ClpSimplex::primal(0);
2334  }
2335  // Get list of non linear columns
2336  char * markNonlinear = new char[numberColumns];
2337  memset(markNonlinear,0,numberColumns);
2338  int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
2339 
2340  if (!numberNonLinearColumns) {
2341    delete [] markNonlinear;
2342    // no nonlinear part
2343    return ClpSimplex::primal(0);
2344  }
2345  int iColumn;
2346  int * listNonLinearColumn = new int[numberNonLinearColumns];
2347  numberNonLinearColumns=0;
2348  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2349    if(markNonlinear[iColumn])
2350      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
2351  }
2352  // Replace objective
2353  ClpObjective * trueObjective = objective_;
2354  objective_=new ClpLinearObjective(NULL,numberColumns);
2355  double * objective = this->objective();
2356 
2357  // get feasible
2358  if (this->status()<0||numberPrimalInfeasibilities())
2359    ClpSimplex::primal(1);
2360  // still infeasible
2361  if (numberPrimalInfeasibilities()) {
2362    delete [] listNonLinearColumn;
2363    return 0;
2364  }
2365
2366  int jNon;
2367  int * last[3];
2368 
2369  double * trust = new double[numberNonLinearColumns];
2370  double * trueLower = new double[numberNonLinearColumns];
2371  double * trueUpper = new double[numberNonLinearColumns];
2372  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2373    iColumn=listNonLinearColumn[jNon];
2374    trust[jNon]=0.5;
2375    trueLower[jNon]=columnLower[iColumn];
2376    trueUpper[jNon]=columnUpper[iColumn];
2377    if (solution[iColumn]<trueLower[jNon])
2378      solution[iColumn]=trueLower[jNon];
2379    else if (solution[iColumn]>trueUpper[jNon])
2380      solution[iColumn]=trueUpper[jNon];
2381  }
2382  int saveLogLevel=logLevel();
2383  int iPass;
2384  double lastObjective=1.0e31;
2385  double * saveSolution = new double [numberColumns];
2386  double * saveRowSolution = new double [numberRows];
2387  memset(saveRowSolution,0,numberRows*sizeof(double));
2388  double * savePi = new double [numberRows];
2389  double * safeSolution = new double [numberColumns];
2390  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
2391#define MULTIPLE 0
2392#if MULTIPLE>2
2393  // Duplication but doesn't really matter
2394  double * saveSolutionM[MULTIPLE};
2395  for (jNon=0;jNon<MULTIPLE;jNon++) {
2396    saveSolutionM[jNon]=new double[numberColumns];
2397    memcpy(saveSolutionM,solution,numberColumns*sizeof(double));
2398  }
2399#endif
2400  double targetDrop=1.0e31;
2401  double objectiveOffset;
2402  getDblParam(ClpObjOffset,objectiveOffset);
2403  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
2404  for (iPass=0;iPass<3;iPass++) {
2405    last[iPass]=new int[numberNonLinearColumns];
2406    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
2407      last[iPass][jNon]=0;
2408  }
2409  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2410  int goodMove=-2;
2411  char * statusCheck = new char[numberColumns];
2412  double * changeRegion = new double [numberColumns];
2413  double offset=0.0;
2414  double objValue=0.0;
2415  for (iPass=0;iPass<numberPasses;iPass++) {
2416    // redo objective
2417    offset=0.0;
2418    objValue=-objectiveOffset;
2419    double theta=-1.0;
2420    double maxTheta=COIN_DBL_MAX;
2421    //maxTheta=1.0;
2422    if (iPass) {
2423      int jNon=0;
2424      for (iColumn=0;iColumn<numberColumns;iColumn++) { 
2425        changeRegion[iColumn]=solution[iColumn]-saveSolution[iColumn];
2426        double alpha = changeRegion[iColumn];
2427        double oldValue = saveSolution[iColumn];
2428        if (markNonlinear[iColumn]==0) {
2429          // linear
2430          if (alpha<-1.0e-15) {
2431            // variable going towards lower bound
2432            double bound = columnLower[iColumn];
2433            oldValue -= bound;
2434            if (oldValue+maxTheta*alpha<0.0) {
2435              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2436            }
2437          } else if (alpha>1.0e-15) {
2438            // variable going towards upper bound
2439            double bound = columnUpper[iColumn];
2440            oldValue = bound-oldValue;
2441            if (oldValue-maxTheta*alpha<0.0) {
2442              maxTheta = CoinMax(0.0,oldValue/alpha);
2443            }
2444          }
2445        } else {
2446          // nonlinear
2447          if (alpha<-1.0e-15) {
2448            // variable going towards lower bound
2449            double bound = trueLower[jNon];
2450            oldValue -= bound;
2451            if (oldValue+maxTheta*alpha<0.0) {
2452              maxTheta = CoinMax(0.0,oldValue/(-alpha));
2453            }
2454          } else if (alpha>1.0e-15) {
2455            // variable going towards upper bound
2456            double bound = trueUpper[jNon];
2457            oldValue = bound-oldValue;
2458            if (oldValue-maxTheta*alpha<0.0) {
2459              maxTheta = CoinMax(0.0,oldValue/alpha);
2460            }
2461          }
2462          jNon++;
2463        }
2464      }
2465      // make sure both accurate
2466      memset(rowActivity_,0,numberRows_*sizeof(double));
2467      times(1.0,solution,rowActivity_);
2468      memset(saveRowSolution,0,numberRows_*sizeof(double));
2469      times(1.0,saveSolution,saveRowSolution);
2470      for (int iRow=0;iRow<numberRows;iRow++) { 
2471        double alpha =rowActivity_[iRow]-saveRowSolution[iRow];
2472        double oldValue = saveRowSolution[iRow];
2473        if (alpha<-1.0e-15) {
2474          // variable going towards lower bound
2475          double bound = rowLower_[iRow];
2476          oldValue -= bound;
2477          if (oldValue+maxTheta*alpha<0.0) {
2478            maxTheta = CoinMax(0.0,oldValue/(-alpha));
2479          }
2480        } else if (alpha>1.0e-15) {
2481          // variable going towards upper bound
2482          double bound = rowUpper_[iRow];
2483          oldValue = bound-oldValue;
2484          if (oldValue-maxTheta*alpha<0.0) {
2485            maxTheta = CoinMax(0.0,oldValue/alpha);
2486          }
2487        }
2488      }
2489    } else {
2490      for (iColumn=0;iColumn<numberColumns;iColumn++) {
2491        changeRegion[iColumn]=0.0;
2492        saveSolution[iColumn]=solution[iColumn];
2493      }
2494      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
2495    }
2496    // get current value anyway
2497    double predictedObj,thetaObj;
2498    double maxTheta2 = 2.0; // to work out a b c
2499    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
2500                                             objValue,predictedObj,thetaObj);
2501   
2502    if (goodMove>=0) {
2503      theta = CoinMin(theta2,maxTheta);
2504#ifdef CLP_DEBUG
2505      if (handler_->logLevel()&32) 
2506        printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
2507               theta,objValue,thetaObj,predictedObj);
2508#endif
2509      if (theta>0.0&&theta<=1.0) {
2510        // update solution
2511        double lambda = 1.0-theta;
2512        for (iColumn=0;iColumn<numberColumns;iColumn++) 
2513          solution[iColumn] = lambda * saveSolution[iColumn] 
2514            + theta * solution[iColumn];
2515        memset(rowActivity_,0,numberRows_*sizeof(double));
2516        times(1.0,solution,rowActivity_);
2517        if (lambda>0.999) {
2518          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
2519          memcpy(status_,saveStatus,numberRows+numberColumns);
2520        }
2521        // Do local minimization
2522#define LOCAL
2523#ifdef LOCAL
2524        int saveScaling = scalingFlag();
2525        scaling(0);
2526        int savePerturbation=perturbation_;
2527        perturbation_=100;
2528        if (saveLogLevel==1)
2529          setLogLevel(0);
2530        int status=startup(1);
2531        if (!status) {
2532          int numberTotal = numberRows_+numberColumns_;
2533          // resize arrays
2534          for (int i=0;i<4;i++) {
2535            rowArray_[i]->reserve(CoinMax(numberRows_+numberColumns_,rowArray_[i]->capacity()));
2536          }
2537          CoinIndexedVector * longArray = rowArray_[3];
2538          CoinIndexedVector * rowArray = rowArray_[0];
2539          //CoinIndexedVector * columnArray = columnArray_[0];
2540          CoinIndexedVector * spare = rowArray_[1];
2541          double * work=longArray->denseVector();
2542          int * which=longArray->getIndices();
2543          int nPass=100;
2544          //bool conjugate=false;
2545          // Put back true bounds
2546          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2547            int iColumn=listNonLinearColumn[jNon];
2548            double value;
2549            value=trueLower[jNon];
2550            trueLower[jNon]=lower_[iColumn];
2551            lower_[iColumn]=value;
2552            value=trueUpper[jNon];
2553            trueUpper[jNon]=upper_[iColumn];
2554            upper_[iColumn]=value;
2555            switch(getStatus(iColumn)) {
2556             
2557            case basic:
2558            case isFree:
2559            case superBasic:
2560              break;
2561            case isFixed:
2562            case atUpperBound:
2563            case atLowerBound:
2564              if (solution_[iColumn]>lower_[iColumn]+primalTolerance_) {
2565                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2566                  setStatus(iColumn,superBasic);
2567                } else {
2568                  setStatus(iColumn,atUpperBound);
2569                }
2570              } else {
2571                if(solution_[iColumn]<upper_[iColumn]-primalTolerance_) {
2572                  setStatus(iColumn,atLowerBound);
2573                } else {
2574                  setStatus(iColumn,isFixed);
2575                }
2576              }
2577              break;
2578            }
2579          }
2580          for (int jPass=0;jPass<nPass;jPass++) {
2581            trueObjective->reducedGradient(this,dj_,true);
2582            // get direction vector in longArray
2583            longArray->clear();
2584            double normFlagged=0.0;
2585            double normUnflagged=0.0;
2586            int numberNonBasic=0;
2587            directionVector(longArray,spare,rowArray,0,
2588                            normFlagged,normUnflagged,numberNonBasic);
2589            if (normFlagged+normUnflagged<1.0e-8)
2590              break;  //looks optimal
2591            double djNorm=0.0;
2592            int iIndex;
2593            for (iIndex=0;iIndex<numberNonBasic;iIndex++) {
2594              int iSequence = which[iIndex];
2595              double alpha = work[iSequence];
2596              djNorm += alpha*alpha;
2597            }
2598            //if (!jPass)
2599            //printf("dj norm %g - %d \n",djNorm,numberNonBasic);
2600            //int number=longArray->getNumElements();
2601            if (!numberNonBasic) {
2602              // looks optimal
2603              break;
2604            }
2605            int bestSequence=-1;
2606            double theta = 1.0e30;
2607            int iSequence;
2608            for (iSequence=0;iSequence<numberTotal;iSequence++) {
2609              double alpha = work[iSequence];
2610              double oldValue = solution_[iSequence];
2611              if (alpha<-1.0e-15) {
2612                // variable going towards lower bound
2613                double bound = lower_[iSequence];
2614                oldValue -= bound;
2615                if (oldValue+theta*alpha<0.0) {
2616                  bestSequence=iSequence;
2617                  theta = CoinMax(0.0,oldValue/(-alpha));
2618                }
2619              } else if (alpha>1.0e-15) {
2620                // variable going towards upper bound
2621                double bound = upper_[iSequence];
2622                oldValue = bound-oldValue;
2623                if (oldValue-theta*alpha<0.0) {
2624                  bestSequence=iSequence;
2625                  theta = CoinMax(0.0,oldValue/alpha);
2626                }
2627              }
2628            }
2629            // Now find minimum of function
2630            double currentObj;
2631            double predictedObj;
2632            double thetaObj;
2633            // need to use true objective
2634            double * saveCost = cost_;
2635            cost_=NULL;
2636            double objTheta =trueObjective->stepLength(this,solution_,work,theta,
2637                                                       currentObj,predictedObj,thetaObj);
2638            cost_=saveCost;
2639#ifdef CLP_DEBUG
2640            if (handler_->logLevel()&32) 
2641              printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2642#endif
2643            //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2644            //printf("objTheta %g theta %g\n",objTheta,theta);
2645            if (theta>objTheta) { 
2646              theta=objTheta;
2647              thetaObj=predictedObj;
2648            }
2649            // update one used outside
2650            objValue=currentObj;
2651            if (theta>1.0e-9&&
2652                (currentObj-thetaObj<-CoinMax(1.0e-8,1.0e-15*fabs(currentObj))||jPass<5)) {
2653              // Update solution
2654              for (iSequence=0;iSequence<numberTotal;iSequence++) {
2655                double alpha = work[iSequence];
2656                if (alpha) {
2657                  double value = solution_[iSequence]+theta*alpha;
2658                  solution_[iSequence]=value;
2659                  switch(getStatus(iSequence)) {
2660                   
2661                  case basic:
2662                  case isFixed:
2663                  case isFree:
2664                    break;
2665                  case atUpperBound:
2666                  case atLowerBound:
2667                  case superBasic:
2668                    if (value>lower_[iSequence]+primalTolerance_) {
2669                      if(value<upper_[iSequence]-primalTolerance_) {
2670                        setStatus(iSequence,superBasic);
2671                      } else {
2672                        setStatus(iSequence,atUpperBound);
2673                      }
2674                    } else {
2675                      if(value<upper_[iSequence]-primalTolerance_) {
2676                        setStatus(iSequence,atLowerBound);
2677                      } else {
2678                        setStatus(iSequence,isFixed);
2679                      }
2680                    }
2681                    break;
2682                  }
2683                }
2684              }
2685            } else {
2686              break;
2687            }
2688          }
2689          // Put back fake bounds
2690          for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2691            int iColumn=listNonLinearColumn[jNon];
2692            double value;
2693            value=trueLower[jNon];
2694            trueLower[jNon]=lower_[iColumn];
2695            lower_[iColumn]=value;
2696            value=trueUpper[jNon];
2697            trueUpper[jNon]=upper_[iColumn];
2698            upper_[iColumn]=value;
2699          }
2700        }
2701        problemStatus_=0;
2702        finish();
2703        scaling(saveScaling);
2704        perturbation_=savePerturbation;
2705        setLogLevel(saveLogLevel);
2706#endif
2707        // redo rowActivity
2708        memset(rowActivity_,0,numberRows_*sizeof(double));
2709        times(1.0,solution,rowActivity_);
2710        if (theta>0.99999&&theta2<1.9) {
2711          // If big changes then tighten
2712          /*  thetaObj is objvalue + a*2*2 +b*2
2713              predictedObj is objvalue + a*theta2*theta2 +b*theta2
2714          */
2715          double rhs1 = thetaObj-objValue;
2716          double rhs2 = predictedObj-objValue;
2717          double subtractB = theta2*0.5;
2718          double a = (rhs2-subtractB*rhs1)/(theta2*theta2-4.0*subtractB);
2719          double b = 0.5*(rhs1 - 4.0*a);
2720          if (fabs(a+b)>1.0e-2) {
2721            // tighten all
2722            goodMove=-1;
2723          }
2724        }
2725      }
2726    }
2727    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
2728           numberColumns*sizeof(double));
2729    // could do faster
2730    trueObjective->stepLength(this,solution,changeRegion,0.0,
2731                                             objValue,predictedObj,thetaObj);
2732    setDblParam(ClpObjOffset,objectiveOffset+offset);
2733    int * temp=last[2];
2734    last[2]=last[1];
2735    last[1]=last[0];
2736    last[0]=temp;
2737    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2738      iColumn=listNonLinearColumn[jNon];
2739      double change = solution[iColumn]-saveSolution[iColumn];
2740      if (change<-1.0e-5) {
2741        if (fabs(change+trust[jNon])<1.0e-5) 
2742          temp[jNon]=-1;
2743        else
2744          temp[jNon]=-2;
2745      } else if(change>1.0e-5) {
2746        if (fabs(change-trust[jNon])<1.0e-5) 
2747          temp[jNon]=1;
2748        else
2749          temp[jNon]=2;
2750      } else {
2751        temp[jNon]=0;
2752      }
2753    } 
2754    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2755    double maxDelta=0.0;
2756    if (goodMove>=0) {
2757      if (objValue<=lastObjective+1.0e-15*fabs(lastObjective)) 
2758        goodMove=1;
2759      else
2760        goodMove=0;
2761    } else {
2762      maxDelta=1.0e10;
2763    }
2764    double maxGap=0.0;
2765    int numberSmaller=0;
2766    int numberSmaller2=0;
2767    int numberLarger=0;
2768    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2769      iColumn=listNonLinearColumn[jNon];
2770      maxDelta = CoinMax(maxDelta,
2771                     fabs(solution[iColumn]-saveSolution[iColumn]));
2772      if (goodMove>0) {
2773        if (last[0][jNon]*last[1][jNon]<0) {
2774          // halve
2775          trust[jNon] *= 0.5;
2776          numberSmaller2++;
2777        } else {
2778          if (last[0][jNon]==last[1][jNon]&&
2779              last[0][jNon]==last[2][jNon])
2780            trust[jNon] = CoinMin(1.5*trust[jNon],1.0e6); 
2781          numberLarger++;
2782        }
2783      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
2784        trust[jNon] *= 0.2;
2785        numberSmaller++;
2786      }
2787      maxGap = CoinMax(maxGap,trust[jNon]);
2788    }
2789#ifdef CLP_DEBUG
2790          if (handler_->logLevel()&32) 
2791            std::cout<<"largest gap is "<<maxGap<<" "
2792                     <<numberSmaller+numberSmaller2<<" reduced ("
2793                     <<numberSmaller<<" badMove ), "
2794                     <<numberLarger<<" increased"<<std::endl;
2795#endif
2796    if (iPass>10000) {
2797      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
2798        trust[jNon] *=0.0001;
2799    }
2800    if (goodMove>0) {
2801      double drop = lastObjective-objValue;
2802      handler_->message(CLP_SLP_ITER,messages_)
2803        <<iPass<<objValue-objectiveOffset
2804        <<drop<<maxDelta
2805        <<CoinMessageEol;
2806      if (iPass>20&&drop<1.0e-12*fabs(objValue))
2807        drop=0.999e-4; // so will exit
2808      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&theta<0.99999) {
2809        if (handler_->logLevel()>1) 
2810          std::cout<<"Exiting as maxDelta < tolerance and small drop"<<std::endl;
2811        break;
2812      }
2813    } else if (!numberSmaller&&iPass>1) {
2814      if (handler_->logLevel()>1) 
2815          std::cout<<"Exiting as all gaps small"<<std::endl;
2816        break;
2817    }
2818    if (!iPass)
2819      goodMove=1;
2820    targetDrop=0.0;
2821    double * r = this->dualColumnSolution();
2822    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2823      iColumn=listNonLinearColumn[jNon];
2824      columnLower[iColumn]=CoinMax(solution[iColumn]
2825                               -trust[jNon],
2826                               trueLower[jNon]);
2827      columnUpper[iColumn]=CoinMin(solution[iColumn]
2828                               +trust[jNon],
2829                               trueUpper[jNon]);
2830    }
2831    if (iPass) {
2832      // get reduced costs
2833      this->matrix()->transposeTimes(savePi,
2834                                     this->dualColumnSolution());
2835      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2836        iColumn=listNonLinearColumn[jNon];
2837        double dj = objective[iColumn]-r[iColumn];
2838        r[iColumn]=dj;
2839        if (dj<-dualTolerance_) 
2840          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
2841        else if (dj>dualTolerance_)
2842          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
2843      }
2844    } else {
2845      memset(r,0,numberColumns*sizeof(double));
2846    }
2847#if 0
2848    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2849      iColumn=listNonLinearColumn[jNon];
2850      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
2851        columnLower[iColumn]=CoinMax(solution[iColumn],
2852                                 trueLower[jNon]);
2853        columnUpper[iColumn]=CoinMin(solution[iColumn]
2854                                 +trust[jNon],
2855                                 trueUpper[jNon]);
2856      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
2857        columnLower[iColumn]=CoinMax(solution[iColumn]
2858                                 -trust[jNon],
2859                                 trueLower[jNon]);
2860        columnUpper[iColumn]=CoinMin(solution[iColumn],
2861                                 trueUpper[jNon]);
2862      } else {
2863        columnLower[iColumn]=CoinMax(solution[iColumn]
2864                                 -trust[jNon],
2865                                 trueLower[jNon]);
2866        columnUpper[iColumn]=CoinMin(solution[iColumn]
2867                                 +trust[jNon],
2868                                 trueUpper[jNon]);
2869      }
2870    }
2871#endif
2872    if (goodMove) {
2873      memcpy(saveSolution,solution,numberColumns*sizeof(double));
2874      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
2875      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
2876      memcpy(saveStatus,status_,numberRows+numberColumns);
2877#if MULTIPLE>2
2878      double * tempSol=saveSolutionM[0];
2879      for (jNon=0;jNon<MULTIPLE-1;jNon++) {
2880        saveSolutionM[jNon]=saveSolutionM[jNon+1];
2881      }
2882      saveSolutionM[MULTIPLE-1]=tempSol;
2883      memcpy(tempSol,solution,numberColumns*sizeof(double));
2884     
2885#endif
2886     
2887#ifdef CLP_DEBUG
2888      if (handler_->logLevel()&32) 
2889        std::cout<<"Pass - "<<iPass
2890                 <<", target drop is "<<targetDrop
2891                 <<std::endl;
2892#endif
2893      lastObjective = objValue;
2894      if (targetDrop<1.0e-5&&goodMove&&iPass) {
2895        if (handler_->logLevel()>1) 
2896        printf("Exiting on target drop %g\n",targetDrop);
2897        break;
2898      }
2899#ifdef CLP_DEBUG
2900      {
2901        double * r = this->dualColumnSolution();
2902        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2903          iColumn=listNonLinearColumn[jNon];
2904          if (handler_->logLevel()&32) 
2905            printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
2906                   jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
2907                   r[iColumn],statusCheck[iColumn],columnLower[iColumn],
2908                   columnUpper[iColumn]);
2909        }
2910      }
2911#endif
2912      //setLogLevel(63);
2913      this->scaling(false);
2914      if (saveLogLevel==1)
2915        setLogLevel(0);
2916      ClpSimplex::primal(1);
2917      setLogLevel(saveLogLevel);
2918#ifdef CLP_DEBUG
2919      if (this->status()) {
2920        writeMps("xx.mps");
2921      }
2922#endif
2923      if (this->status()==1) {
2924        // not feasible ! - backtrack and exit
2925        // use safe solution
2926        memcpy(solution,safeSolution,numberColumns*sizeof(double));
2927        memcpy(saveSolution,solution,numberColumns*sizeof(double));
2928        memset(rowActivity_,0,numberRows_*sizeof(double));
2929        times(1.0,solution,rowActivity_);
2930        memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
2931        memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
2932        memcpy(status_,saveStatus,numberRows+numberColumns);
2933        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2934          iColumn=listNonLinearColumn[jNon];
2935          columnLower[iColumn]=CoinMax(solution[iColumn]
2936                                       -trust[jNon],
2937                                       trueLower[jNon]);
2938          columnUpper[iColumn]=CoinMin(solution[iColumn]
2939                                       +trust[jNon],
2940                                       trueUpper[jNon]);
2941        }
2942        break;
2943      } else {
2944        // save in case problems
2945        memcpy(safeSolution,solution,numberColumns*sizeof(double));
2946      }
2947      if (problemStatus_==3)
2948        break; // probably user interrupt
2949      goodMove=1;
2950    } else {
2951      // bad pass - restore solution
2952#ifdef CLP_DEBUG
2953      if (handler_->logLevel()&32) 
2954        printf("Backtracking\n");
2955#endif
2956      memcpy(solution,saveSolution,numberColumns*sizeof(double));
2957      memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
2958      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
2959      memcpy(status_,saveStatus,numberRows+numberColumns);
2960      iPass--;
2961      goodMove=-1;
2962    }
2963  }
2964#if MULTIPLE>2
2965  for (jNon=0;jNon<MULTIPLE;jNon++) {
2966    delete [] saveSolutionM[jNon];
2967  }
2968#endif
2969  // restore solution
2970  memcpy(solution,saveSolution,numberColumns*sizeof(double));
2971  memcpy(rowActivity_,saveRowSolution,numberRows*sizeof(double));
2972  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2973    iColumn=listNonLinearColumn[jNon];
2974    columnLower[iColumn]=CoinMax(solution[iColumn],
2975                             trueLower[jNon]);
2976    columnUpper[iColumn]=CoinMin(solution[iColumn],
2977                             trueUpper[jNon]);
2978  }
2979  delete [] markNonlinear;
2980  delete [] statusCheck;
2981  delete [] savePi;
2982  delete [] saveStatus;
2983  // redo objective
2984  memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
2985         numberColumns*sizeof(double));
2986  ClpSimplex::primal(1);
2987  delete objective_;
2988  objective_=trueObjective;
2989  // redo values
2990  setDblParam(ClpObjOffset,objectiveOffset);
2991  objectiveValue_ += whichWay*offset;
2992  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
2993    iColumn=listNonLinearColumn[jNon];
2994    columnLower[iColumn]= trueLower[jNon];
2995    columnUpper[iColumn]= trueUpper[jNon];
2996  }
2997  delete [] saveSolution;
2998  delete [] safeSolution;
2999  delete [] saveRowSolution;
3000  for (iPass=0;iPass<3;iPass++) 
3001    delete [] last[iPass];
3002  delete [] trust;
3003  delete [] trueUpper;
3004  delete [] trueLower;
3005  delete [] listNonLinearColumn;
3006  delete [] changeRegion;
3007  // temp
3008  //setLogLevel(63);
3009  return 0;
3010}
Note: See TracBrowser for help on using the repository browser.