source: tags/move-to-subversion/ClpSimplexNonlinear.cpp @ 1355

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

for quadratic

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