source: trunk/Clp/src/ClpSimplexPrimal.cpp @ 1304

Last change on this file since 1304 was 1303, checked in by forrest, 12 years ago

fix saying feasible for Scip

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 109.2 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/* Notes on implementation of primal simplex algorithm.
5
6   When primal feasible(A):
7
8   If dual feasible, we are optimal.  Otherwise choose an infeasible
9   basic variable to enter basis from a bound (B).  We now need to find an
10   outgoing variable which will leave problem primal feasible so we get
11   the column of the tableau corresponding to the incoming variable
12   (with the correct sign depending if variable will go up or down).
13
14   We now perform a ratio test to determine which outgoing variable will
15   preserve primal feasibility (C).  If no variable found then problem
16   is unbounded (in primal sense).  If there is a variable, we then
17   perform pivot and repeat.  Trivial?
18
19   -------------------------------------------
20
21   A) How do we get primal feasible?  All variables have fake costs
22   outside their feasible region so it is trivial to declare problem
23   feasible.  OSL did not have a phase 1/phase 2 approach but
24   instead effectively put an extra cost on infeasible basic variables
25   I am taking the same approach here, although it is generalized
26   to allow for non-linear costs and dual information.
27
28   In OSL, this weight was changed heuristically, here at present
29   it is only increased if problem looks finished.  If problem is
30   feasible I check for unboundedness.  If not unbounded we
31   could play with going into dual.  As long as weights increase
32   any algorithm would be finite.
33   
34   B) Which incoming variable to choose is a virtual base class.
35   For difficult problems steepest edge is preferred while for
36   very easy (large) problems we will need partial scan.
37
38   C) Sounds easy, but this is hardest part of algorithm.
39      1) Instead of stopping at first choice, we may be able
40      to allow that variable to go through bound and if objective
41      still improving choose again.  These mini iterations can
42      increase speed by orders of magnitude but we may need to
43      go to more of a bucket choice of variable rather than looking
44      at them one by one (for speed).
45      2) Accuracy.  Basic infeasibilities may be less than
46      tolerance.  Pivoting on these makes objective go backwards.
47      OSL modified cost so a zero move was made, Gill et al
48      modified so a strictly positive move was made.
49      The two problems are that re-factorizations can
50      change rinfeasibilities above and below tolerances and that when
51      finished we need to reset costs and try again.
52      3) Degeneracy.  Gill et al helps but may not be enough.  We
53      may need more.  Also it can improve speed a lot if we perturb
54      the rhs and bounds significantly. 
55
56  References:
57     Forrest and Goldfarb, Steepest-edge simplex algorithms for
58       linear programming - Mathematical Programming 1992
59     Forrest and Tomlin, Implementing the simplex method for
60       the Optimization Subroutine Library - IBM Systems Journal 1992
61     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
62       Procedure for Linear and Nonlinear Programming SOL report 1988
63
64
65  TODO:
66 
67  a) Better recovery procedures.  At present I never check on forward
68     progress.  There is checkpoint/restart with reducing
69     re-factorization frequency, but this is only on singular
70     factorizations.
71  b) Fast methods for large easy problems (and also the option for
72     the code to automatically choose which method).
73  c) We need to be able to stop in various ways for OSI - this
74     is fairly easy.
75
76 */
77
78
79#include "CoinPragma.hpp"
80
81#include <math.h>
82
83#include "CoinHelperFunctions.hpp"
84#include "ClpSimplexPrimal.hpp"
85#include "ClpFactorization.hpp"
86#include "ClpNonLinearCost.hpp"
87#include "CoinPackedMatrix.hpp"
88#include "CoinIndexedVector.hpp"
89#include "ClpPrimalColumnPivot.hpp"
90#include "ClpMessage.hpp"
91#include "ClpEventHandler.hpp"
92#include <cfloat>
93#include <cassert>
94#include <string>
95#include <stdio.h>
96#include <iostream>
97// primal
98int ClpSimplexPrimal::primal (int ifValuesPass , int startFinishOptions)
99{
100
101  /*
102      Method
103
104     It tries to be a single phase approach with a weight of 1.0 being
105     given to getting optimal and a weight of infeasibilityCost_ being
106     given to getting primal feasible.  In this version I have tried to
107     be clever in a stupid way.  The idea of fake bounds in dual
108     seems to work so the primal analogue would be that of getting
109     bounds on reduced costs (by a presolve approach) and using
110     these for being above or below feasible region.  I decided to waste
111     memory and keep these explicitly.  This allows for non-linear
112     costs!
113
114     The code is designed to take advantage of sparsity so arrays are
115     seldom zeroed out from scratch or gone over in their entirety.
116     The only exception is a full scan to find incoming variable for
117     Dantzig row choice.  For steepest edge we keep an updated list
118     of dual infeasibilities (actually squares). 
119     On easy problems we don't need full scan - just
120     pick first reasonable.
121
122     One problem is how to tackle degeneracy and accuracy.  At present
123     I am using the modification of costs which I put in OSL and which was
124     extended by Gill et al.  I am still not sure of the exact details.
125
126     The flow of primal is three while loops as follows:
127
128     while (not finished) {
129
130       while (not clean solution) {
131
132          Factorize and/or clean up solution by changing bounds so
133          primal feasible.  If looks finished check fake primal bounds.
134          Repeat until status is iterating (-1) or finished (0,1,2)
135
136       }
137
138       while (status==-1) {
139
140         Iterate until no pivot in or out or time to re-factorize.
141
142         Flow is:
143
144         choose pivot column (incoming variable).  if none then
145         we are primal feasible so looks as if done but we need to
146         break and check bounds etc.
147
148         Get pivot column in tableau
149
150         Choose outgoing row.  If we don't find one then we look
151         primal unbounded so break and check bounds etc.  (Also the
152         pivot tolerance is larger after any iterations so that may be
153         reason)
154
155         If we do find outgoing row, we may have to adjust costs to
156         keep going forwards (anti-degeneracy).  Check pivot will be stable
157         and if unstable throw away iteration and break to re-factorize.
158         If minor error re-factorize after iteration.
159
160         Update everything (this may involve changing bounds on
161         variables to stay primal feasible.
162
163       }
164
165     }
166
167     At present we never check we are going forwards.  I overdid that in
168     OSL so will try and make a last resort.
169
170     Needs partial scan pivot in option.
171
172     May need other anti-degeneracy measures, especially if we try and use
173     loose tolerances as a way to solve in fewer iterations.
174
175     I like idea of dynamic scaling.  This gives opportunity to decouple
176     different implications of scaling for accuracy, iteration count and
177     feasibility tolerance.
178
179  */
180
181  algorithm_ = +1;
182  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
183
184  // save data
185  ClpDataSave data = saveData();
186  matrix_->refresh(this); // make sure matrix okay
187 
188  // Save so can see if doing after dual
189  int initialStatus=problemStatus_;
190  int initialIterations = numberIterations_;
191  int initialNegDjs=-1;
192  // initialize - maybe values pass and algorithm_ is +1
193#if 0
194  // if so - put in any superbasic costed slacks
195  if (ifValuesPass&&specialOptions_<0x01000000) {
196    // Get column copy
197    const CoinPackedMatrix * columnCopy = matrix();
198    const int * row = columnCopy->getIndices();
199    const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
200    const int * columnLength = columnCopy->getVectorLengths();
201    //const double * element = columnCopy->getElements();
202    int n=0;
203    for (int iColumn = 0;iColumn<numberColumns_;iColumn++) {
204      if (columnLength[iColumn]==1) {
205        Status status = getColumnStatus(iColumn);
206        if (status!=basic&&status!=isFree) {
207          double value = columnActivity_[iColumn];
208          if (fabs(value-columnLower_[iColumn])>primalTolerance_&&
209              fabs(value-columnUpper_[iColumn])>primalTolerance_) {
210            int iRow = row[columnStart[iColumn]];
211            if (getRowStatus(iRow)==basic) {
212              setRowStatus(iRow,superBasic);
213              setColumnStatus(iColumn,basic);
214              n++;
215            }
216          }
217        }   
218      }
219    }
220    printf("%d costed slacks put in basis\n",n);
221  }
222#endif
223  if (!startup(ifValuesPass,startFinishOptions)) {
224   
225    // Set average theta
226    nonLinearCost_->setAverageTheta(1.0e3);
227    int lastCleaned=0; // last time objective or bounds cleaned up
228   
229    // Say no pivot has occurred (for steepest edge and updates)
230    pivotRow_=-2;
231   
232    // This says whether to restore things etc
233    int factorType=0;
234    if (problemStatus_<0&&perturbation_<100) {
235      perturb(0);
236      // Can't get here if values pass
237      assert (!ifValuesPass);
238      gutsOfSolution(NULL,NULL);
239      if (handler_->logLevel()>2) {
240        handler_->message(CLP_SIMPLEX_STATUS,messages_)
241          <<numberIterations_<<objectiveValue();
242        handler_->printing(sumPrimalInfeasibilities_>0.0)
243          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
244        handler_->printing(sumDualInfeasibilities_>0.0)
245          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
246        handler_->printing(numberDualInfeasibilitiesWithoutFree_
247                           <numberDualInfeasibilities_)
248                             <<numberDualInfeasibilitiesWithoutFree_;
249        handler_->message()<<CoinMessageEol;
250      }
251    }
252    ClpSimplex * saveModel=NULL;
253    int stopSprint=-1;
254    int sprintPass=0;
255    int reasonableSprintIteration=0;
256    int lastSprintIteration=0;
257    double lastObjectiveValue=COIN_DBL_MAX;
258    // Start check for cycles
259    progress_.fillFromModel(this);
260    progress_.startCheck();
261    /*
262      Status of problem:
263      0 - optimal
264      1 - infeasible
265      2 - unbounded
266      -1 - iterating
267      -2 - factorization wanted
268      -3 - redo checking without factorization
269      -4 - looks infeasible
270      -5 - looks unbounded
271    */
272    while (problemStatus_<0) {
273      int iRow,iColumn;
274      // clear
275      for (iRow=0;iRow<4;iRow++) {
276        rowArray_[iRow]->clear();
277      }   
278     
279      for (iColumn=0;iColumn<2;iColumn++) {
280        columnArray_[iColumn]->clear();
281      }   
282     
283      // give matrix (and model costs and bounds a chance to be
284      // refreshed (normally null)
285      matrix_->refresh(this);
286      // If getting nowhere - why not give it a kick
287#if 1
288      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)&&(specialOptions_&4)==0
289          &&initialStatus!=10) {
290        perturb(1);
291        matrix_->rhsOffset(this,true,false);
292      }
293#endif
294      // If we have done no iterations - special
295      if (lastGoodIteration_==numberIterations_&&factorType)
296        factorType=3;
297      if (saveModel) {
298        // Doing sprint
299        if (sequenceIn_<0||numberIterations_>=stopSprint) {
300          problemStatus_=-1;
301          originalModel(saveModel);
302          saveModel=NULL;
303          if (sequenceIn_<0&&numberIterations_<reasonableSprintIteration&&
304              sprintPass>100)
305            primalColumnPivot_->switchOffSprint();
306          //lastSprintIteration=numberIterations_;
307          printf("End small model\n");
308        }
309      }
310     
311      // may factorize, checks if problem finished
312      statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
313      if (initialStatus==10) {
314        // cleanup phase
315        if(initialIterations != numberIterations_) {
316          if (numberDualInfeasibilities_>10000&&numberDualInfeasibilities_>10*initialNegDjs) {
317            // getting worse - try perturbing
318            if (perturbation_<101&&(specialOptions_&4)==0) {
319              perturb(1);
320              matrix_->rhsOffset(this,true,false);
321              statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
322            }
323          }
324        } else {
325          // save number of negative djs
326          if (!numberPrimalInfeasibilities_)
327            initialNegDjs=numberDualInfeasibilities_;
328          // make sure weight won't be changed
329          if (infeasibilityCost_==1.0e10)
330            infeasibilityCost_=1.000001e10;
331        }
332      }
333      // See if sprint says redo because of problems
334      if (numberDualInfeasibilities_==-776) {
335        // Need new set of variables
336        problemStatus_=-1;
337        originalModel(saveModel);
338        saveModel=NULL;
339        //lastSprintIteration=numberIterations_;
340        printf("End small model after\n");
341        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
342      } 
343      int numberSprintIterations=0;
344      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
345      if (problemStatus_==777) {
346        // problems so do one pass with normal
347        problemStatus_=-1;
348        originalModel(saveModel);
349        saveModel=NULL;
350        // Skip factorization
351        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
352        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
353      } else if (problemStatus_<0&&!saveModel&&numberSprintColumns&&firstFree_<0) {
354        int numberSort=0;
355        int numberFixed=0;
356        int numberBasic=0;
357        reasonableSprintIteration = numberIterations_ + 100;
358        int * whichColumns = new int[numberColumns_];
359        double * weight = new double[numberColumns_];
360        int numberNegative=0;
361        double sumNegative = 0.0;
362        // now massage weight so all basic in plus good djs
363        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
364          double dj = dj_[iColumn];
365          switch(getColumnStatus(iColumn)) {
366           
367          case basic:
368            dj = -1.0e50;
369            numberBasic++;
370            break;
371          case atUpperBound:
372            dj = -dj;
373            break;
374          case isFixed:
375            dj=1.0e50;
376            numberFixed++;
377            break;
378          case atLowerBound:
379            dj = dj;
380            break;
381          case isFree:
382            dj = -100.0*fabs(dj);
383              break;
384          case superBasic:
385            dj = -100.0*fabs(dj);
386            break;
387          }
388          if (dj<-dualTolerance_&&dj>-1.0e50) {
389            numberNegative++;
390            sumNegative -= dj;
391          }
392          weight[iColumn]=dj;
393          whichColumns[iColumn] = iColumn;
394        }
395        handler_->message(CLP_SPRINT,messages_)
396          <<sprintPass<<numberIterations_-lastSprintIteration<<objectiveValue()<<sumNegative
397          <<numberNegative
398          <<CoinMessageEol;
399        sprintPass++;
400        lastSprintIteration=numberIterations_;
401        if (objectiveValue()*optimizationDirection_>lastObjectiveValue-1.0e-7&&sprintPass>5) {
402          // switch off
403          printf("Switching off sprint\n");
404          primalColumnPivot_->switchOffSprint();
405        } else {
406          lastObjectiveValue = objectiveValue()*optimizationDirection_;
407          // sort
408          CoinSort_2(weight,weight+numberColumns_,whichColumns);
409          numberSort = CoinMin(numberColumns_-numberFixed,numberBasic+numberSprintColumns);
410          // Sort to make consistent ?
411          std::sort(whichColumns,whichColumns+numberSort);
412          saveModel = new ClpSimplex(this,numberSort,whichColumns);
413          delete [] whichColumns;
414          delete [] weight;
415          // Skip factorization
416          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
417          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
418          stopSprint = numberIterations_+numberSprintIterations;
419          printf("Sprint with %d columns for %d iterations\n",
420                 numberSprintColumns,numberSprintIterations);
421        }
422      }
423     
424      // Say good factorization
425      factorType=1;
426     
427      // Say no pivot has occurred (for steepest edge and updates)
428      pivotRow_=-2;
429
430      // exit if victory declared
431      if (problemStatus_>=0)
432        break;
433     
434      // test for maximum iterations
435      if (hitMaximumIterations()||(ifValuesPass==2&&firstFree_<0)) {
436        problemStatus_=3;
437        break;
438      }
439
440      if (firstFree_<0) {
441        if (ifValuesPass) {
442          // end of values pass
443          ifValuesPass=0;
444          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
445          if (status>=0) {
446            problemStatus_=5;
447            secondaryStatus_=ClpEventHandler::endOfValuesPass;
448            break;
449          }
450        }
451      }
452      // Check event
453      {
454        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
455        if (status>=0) {
456          problemStatus_=5;
457          secondaryStatus_=ClpEventHandler::endOfFactorization;
458          break;
459        }
460      }
461      // Iterate
462      whileIterating(ifValuesPass ? 1 : 0);
463    }
464  }
465  // if infeasible get real values
466  //printf("XXXXY final cost %g\n",infeasibilityCost_);
467  progress_.initialWeight_=0.0;
468  if (problemStatus_==1&&secondaryStatus_!=6) {
469    infeasibilityCost_=0.0;
470    createRim(1+4);
471    nonLinearCost_->checkInfeasibilities(0.0);
472    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
473    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
474    // and get good feasible duals
475    computeDuals(NULL);
476  }
477  // clean up
478  unflag();
479  finish(startFinishOptions);
480  restoreData(data);
481  return problemStatus_;
482}
483/*
484  Reasons to come out:
485  -1 iterations etc
486  -2 inaccuracy
487  -3 slight inaccuracy (and done iterations)
488  -4 end of values pass and done iterations
489  +0 looks optimal (might be infeasible - but we will investigate)
490  +2 looks unbounded
491  +3 max iterations
492*/
493int
494ClpSimplexPrimal::whileIterating(int valuesOption)
495{
496  // Say if values pass
497  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
498  int returnCode=-1;
499  int superBasicType=1;
500  if (valuesOption>1)
501    superBasicType=3;
502  // status stays at -1 while iterating, >=0 finished, -2 to invert
503  // status -3 to go to top without an invert
504  while (problemStatus_==-1) {
505    //#define CLP_DEBUG 1
506#ifdef CLP_DEBUG
507    {
508      int i;
509      // not [1] as has information
510      for (i=0;i<4;i++) {
511        if (i!=1)
512          rowArray_[i]->checkClear();
513      }   
514      for (i=0;i<2;i++) {
515        columnArray_[i]->checkClear();
516      }
517    }     
518#endif
519#if 0
520    {
521      int iPivot;
522      double * array = rowArray_[3]->denseVector();
523      int * index = rowArray_[3]->getIndices();
524      int i;
525      for (iPivot=0;iPivot<numberRows_;iPivot++) {
526        int iSequence = pivotVariable_[iPivot];
527        unpackPacked(rowArray_[3],iSequence);
528        factorization_->updateColumn(rowArray_[2],rowArray_[3]);
529        int number = rowArray_[3]->getNumElements();
530        for (i=0;i<number;i++) {
531          int iRow = index[i];
532          if (iRow==iPivot)
533            assert (fabs(array[i]-1.0)<1.0e-4);
534          else
535            assert (fabs(array[i])<1.0e-4);
536        }
537        rowArray_[3]->clear();
538      }
539    }
540#endif
541#if 0
542    nonLinearCost_->checkInfeasibilities(primalTolerance_);
543    printf("suminf %g number %d\n",nonLinearCost_->sumInfeasibilities(),
544           nonLinearCost_->numberInfeasibilities());
545#endif
546#if CLP_DEBUG>2
547    // very expensive
548    if (numberIterations_>0&&numberIterations_<100&&!ifValuesPass) {
549      handler_->setLogLevel(63);
550      double saveValue = objectiveValue_;
551      double * saveRow1 = new double[numberRows_];
552      double * saveRow2 = new double[numberRows_];
553      CoinMemcpyN(rowReducedCost_,numberRows_,saveRow1);
554      CoinMemcpyN(rowActivityWork_,numberRows_,saveRow2);
555      double * saveColumn1 = new double[numberColumns_];
556      double * saveColumn2 = new double[numberColumns_];
557      CoinMemcpyN(reducedCostWork_,numberColumns_,saveColumn1);
558      CoinMemcpyN(columnActivityWork_,numberColumns_,saveColumn2);
559      gutsOfSolution(NULL,NULL,false);
560      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
561             numberIterations_,
562             saveValue,objectiveValue_,sumPrimalInfeasibilities_);
563      CoinMemcpyN(saveRow1,numberRows_,rowReducedCost_);
564      CoinMemcpyN(saveRow2,numberRows_,rowActivityWork_);
565      CoinMemcpyN(saveColumn1,numberColumns_,reducedCostWork_);
566      CoinMemcpyN(saveColumn2,numberColumns_,columnActivityWork_);
567      delete [] saveRow1;
568      delete [] saveRow2;
569      delete [] saveColumn1;
570      delete [] saveColumn2;
571      objectiveValue_=saveValue;
572    }
573#endif
574    if (!ifValuesPass) {
575      // choose column to come in
576      // can use pivotRow_ to update weights
577      // pass in list of cost changes so can do row updates (rowArray_[1])
578      // NOTE rowArray_[0] is used by computeDuals which is a
579      // slow way of getting duals but might be used
580      primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
581                   columnArray_[0],columnArray_[1]);
582    } else {
583      // in values pass
584      int sequenceIn=nextSuperBasic(superBasicType,columnArray_[0]);
585      if (valuesOption>1)
586        superBasicType=2;
587      if (sequenceIn<0) {
588        // end of values pass - initialize weights etc
589        handler_->message(CLP_END_VALUES_PASS,messages_)
590          <<numberIterations_;
591        primalColumnPivot_->saveWeights(this,5);
592        problemStatus_=-2; // factorize now
593        pivotRow_=-1; // say no weights update
594        returnCode=-4;
595        // Clean up
596        int i;
597        for (i=0;i<numberRows_+numberColumns_;i++) {
598          if (getColumnStatus(i)==atLowerBound||getColumnStatus(i)==isFixed)
599            solution_[i]=lower_[i];
600          else if (getColumnStatus(i)==atUpperBound)
601            solution_[i]=upper_[i];
602        }
603        break;
604      } else {
605        // normal
606        sequenceIn_ = sequenceIn;
607        valueIn_=solution_[sequenceIn_];
608        lowerIn_=lower_[sequenceIn_];
609        upperIn_=upper_[sequenceIn_];
610        dualIn_=dj_[sequenceIn_];
611      }
612    }
613    pivotRow_=-1;
614    sequenceOut_=-1;
615    rowArray_[1]->clear();
616    if (sequenceIn_>=0) {
617      // we found a pivot column
618      assert (!flagged(sequenceIn_));
619#ifdef CLP_DEBUG
620      if ((handler_->logLevel()&32)) {
621        char x = isColumn(sequenceIn_) ? 'C' :'R';
622        std::cout<<"pivot column "<<
623          x<<sequenceWithin(sequenceIn_)<<std::endl;
624      }
625#endif
626#ifdef CLP_DEBUG
627    {
628      int checkSequence=-2077;
629      if (checkSequence>=0&&checkSequence<numberRows_+numberColumns_&&!ifValuesPass) {
630        rowArray_[2]->checkClear();
631        rowArray_[3]->checkClear();
632        double * array = rowArray_[3]->denseVector();
633        int * index = rowArray_[3]->getIndices();
634        unpackPacked(rowArray_[3],checkSequence);
635        factorization_->updateColumnForDebug(rowArray_[2],rowArray_[3]);
636        int number = rowArray_[3]->getNumElements();
637        double dualIn = cost_[checkSequence];
638        int i;
639        for (i=0;i<number;i++) {
640          int iRow = index[i];
641          int iPivot = pivotVariable_[iRow];
642          double alpha = array[i];
643          dualIn -= alpha*cost_[iPivot];
644        }
645        printf("old dj for %d was %g, recomputed %g\n",checkSequence,
646               dj_[checkSequence],dualIn);
647        rowArray_[3]->clear();
648        if (numberIterations_>2000)
649          exit(1);
650      }
651    }
652#endif
653      // do second half of iteration
654      returnCode = pivotResult(ifValuesPass);
655      if (returnCode<-1&&returnCode>-5) {
656        problemStatus_=-2; //
657      } else if (returnCode==-5) {
658        if ((moreSpecialOptions_&16)==0&&factorization_->pivots()) {
659          moreSpecialOptions_ |= 16;
660          problemStatus_=-2;
661        }
662        // otherwise something flagged - continue;
663      } else if (returnCode==2) {
664        problemStatus_=-5; // looks unbounded
665      } else if (returnCode==4) {
666        problemStatus_=-2; // looks unbounded but has iterated
667      } else if (returnCode!=-1) {
668        assert(returnCode==3);
669        if (problemStatus_!=5)
670          problemStatus_=3;
671      }
672    } else {
673      // no pivot column
674#ifdef CLP_DEBUG
675      if (handler_->logLevel()&32)
676        printf("** no column pivot\n");
677#endif
678      if (nonLinearCost_->numberInfeasibilities())
679        problemStatus_=-4; // might be infeasible
680      // Force to re-factorize early next time
681      int numberPivots = factorization_->pivots();
682      forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
683      returnCode=0;
684      break;
685    }
686  }
687  if (valuesOption>1) 
688    columnArray_[0]->setNumElements(0);
689  return returnCode;
690}
691/* Checks if finished.  Updates status */
692void 
693ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type,
694                                          ClpSimplexProgress * progress,
695                                          bool doFactorization,
696                                          int ifValuesPass,
697                                          ClpSimplex * originalModel)
698{
699  int dummy; // for use in generalExpanded
700  int saveFirstFree=firstFree_;
701  // number of pivots done
702  int numberPivots = factorization_->pivots();
703  if (type==2) {
704    // trouble - restore solution
705    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
706    CoinMemcpyN(savedSolution_+numberColumns_ ,
707           numberRows_,rowActivityWork_);
708    CoinMemcpyN(savedSolution_ ,
709           numberColumns_,columnActivityWork_);
710    // restore extra stuff
711    matrix_->generalExpanded(this,6,dummy);
712    forceFactorization_=1; // a bit drastic but ..
713    pivotRow_=-1; // say no weights update
714    changeMade_++; // say change made
715  }
716  int saveThreshold = factorization_->sparseThreshold();
717  int tentativeStatus = problemStatus_;
718  int numberThrownOut=1; // to loop round on bad factorization in values pass
719  double lastSumInfeasibility=COIN_DBL_MAX;
720  if (numberIterations_)
721    lastSumInfeasibility=nonLinearCost_->sumInfeasibilities();
722  int nPass=0;
723  while (numberThrownOut) {
724    int nSlackBasic=0;
725    if (nPass) {
726      for (int i=0;i<numberRows_;i++) {
727        if (getRowStatus(i)==basic)
728          nSlackBasic++;
729      }
730    }
731    nPass++;
732    if (problemStatus_>-3||problemStatus_==-4) {
733      // factorize
734      // later on we will need to recover from singularities
735      // also we could skip if first time
736      // do weights
737      // This may save pivotRow_ for use
738      if (doFactorization)
739        primalColumnPivot_->saveWeights(this,1);
740     
741      if ((type&&doFactorization)||nSlackBasic==numberRows_) {
742        // is factorization okay?
743        int factorStatus = internalFactorize(1);
744        if (factorStatus) {
745          if (solveType_==2+8) {
746            // say odd
747            problemStatus_=5;
748            return;
749          }
750          if (type!=1||largestPrimalError_>1.0e3
751              ||largestDualError_>1.0e3) {
752            // switch off dense
753            int saveDense = factorization_->denseThreshold();
754            factorization_->setDenseThreshold(0);
755            // Go to safe
756            factorization_->pivotTolerance(0.99);
757            // make sure will do safe factorization
758            pivotVariable_[0]=-1;
759            internalFactorize(2);
760            factorization_->setDenseThreshold(saveDense);
761            // restore extra stuff
762            matrix_->generalExpanded(this,6,dummy);
763          } else {
764            // no - restore previous basis
765            // Keep any flagged variables
766            int i;
767            for (i=0;i<numberRows_+numberColumns_;i++) {
768              if (flagged(i))
769                saveStatus_[i] |= 64; //say flagged
770            }
771            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
772            if (numberPivots<=1) {
773              // throw out something
774              if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
775                setFlagged(sequenceIn_);
776              } else if (sequenceOut_>=0&&getStatus(sequenceOut_)!=basic) {
777                setFlagged(sequenceOut_);
778              }
779              double newTolerance = CoinMax(0.5 + 0.499*randomNumberGenerator_.randomDouble(),factorization_->pivotTolerance());
780              factorization_->pivotTolerance(newTolerance);
781            } else {
782              // Go to safe
783              factorization_->pivotTolerance(0.99);
784            }
785            CoinMemcpyN(savedSolution_+numberColumns_ ,
786                   numberRows_,rowActivityWork_);
787            CoinMemcpyN(savedSolution_ ,
788                   numberColumns_,columnActivityWork_);
789            // restore extra stuff
790            matrix_->generalExpanded(this,6,dummy);
791            matrix_->generalExpanded(this,5,dummy);
792            forceFactorization_=1; // a bit drastic but ..
793            type = 2;
794            if (internalFactorize(2)!=0) {
795              largestPrimalError_=1.0e4; // force other type
796            }
797          }
798          changeMade_++; // say change made
799        }
800      }
801      if (problemStatus_!=-4)
802        problemStatus_=-3;
803    }
804    // at this stage status is -3 or -5 if looks unbounded
805    // get primal and dual solutions
806    // put back original costs and then check
807    // createRim(4); // costs do not change
808    // May need to do more if column generation
809    dummy=4;
810    matrix_->generalExpanded(this,9,dummy);
811    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
812    double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
813    if (numberThrownOut||
814        (sumInfeasibility>1.0e7&&sumInfeasibility>100.0*lastSumInfeasibility
815         &&factorization_->pivotTolerance()<0.11)||(largestPrimalError_>1.0e10&&largestDualError_>1.0e10)) {
816      problemStatus_=tentativeStatus;
817      doFactorization=true;
818      if (numberPivots) {
819        // go back
820        numberThrownOut=-1;
821        // trouble - restore solution
822        CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
823        CoinMemcpyN(savedSolution_+numberColumns_ ,
824                    numberRows_,rowActivityWork_);
825        CoinMemcpyN(savedSolution_ ,
826           numberColumns_,columnActivityWork_);
827        // restore extra stuff
828        matrix_->generalExpanded(this,6,dummy);
829        forceFactorization_=1; // a bit drastic but ..
830        // Go to safe
831        factorization_->pivotTolerance(0.99);
832        pivotRow_=-1; // say no weights update
833        changeMade_++; // say change made
834        if (numberPivots==1) {
835          // throw out something
836          if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
837            setFlagged(sequenceIn_);
838          } else if (sequenceOut_>=0&&getStatus(sequenceOut_)!=basic) {
839            setFlagged(sequenceOut_);
840          }
841        }
842        numberPivots=0;
843      }
844    }
845  }
846  // Double check reduced costs if no action
847  if (progress->lastIterationNumber(0)==numberIterations_) {
848    if (primalColumnPivot_->looksOptimal()) {
849      numberDualInfeasibilities_ = 0;
850      sumDualInfeasibilities_ = 0.0;
851    }
852  }
853  // If in primal and small dj give up
854  if ((specialOptions_&1024)!=0&&!numberPrimalInfeasibilities_&&numberDualInfeasibilities_) {
855    double average = sumDualInfeasibilities_/((double) numberDualInfeasibilities_);
856    if (numberIterations_>300&&average<1.0e-4) {
857      numberDualInfeasibilities_ = 0;
858      sumDualInfeasibilities_ = 0.0;
859    }
860  }
861  // Check if looping
862  int loop;
863  if (type!=2&&!ifValuesPass) 
864    loop = progress->looping();
865  else
866    loop=-1;
867  if (loop>=0) {
868    if (!problemStatus_) {
869      // declaring victory
870      numberPrimalInfeasibilities_ = 0;
871      sumPrimalInfeasibilities_ = 0.0;
872    } else {
873      problemStatus_ = loop; //exit if in loop
874      problemStatus_ = 10; // instead - try other algorithm
875      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
876    }
877    problemStatus_ = 10; // instead - try other algorithm
878    return ;
879  } else if (loop<-1) {
880    // Is it time for drastic measures
881    if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
882        progress->oddState()<10&&progress->oddState()>=0) {
883      progress->newOddState();
884      nonLinearCost_->zapCosts();
885    }
886    // something may have changed
887    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
888  }
889  // If progress then reset costs
890  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
891    createRim(4,false); // costs back
892    delete nonLinearCost_;
893    nonLinearCost_ = new ClpNonLinearCost(this);
894    progress->endOddState();
895    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
896  }
897  // Flag to say whether to go to dual to clean up
898  bool goToDual=false;
899  // really for free variables in
900  //if((progressFlag_&2)!=0)
901  //problemStatus_=-1;;
902  progressFlag_ = 0; //reset progress flag
903
904  handler_->message(CLP_SIMPLEX_STATUS,messages_)
905    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
906  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
907    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
908  handler_->printing(sumDualInfeasibilities_>0.0)
909    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
910  handler_->printing(numberDualInfeasibilitiesWithoutFree_
911                     <numberDualInfeasibilities_)
912                       <<numberDualInfeasibilitiesWithoutFree_;
913  handler_->message()<<CoinMessageEol;
914  if (!primalFeasible()) {
915    nonLinearCost_->checkInfeasibilities(primalTolerance_);
916    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
917    nonLinearCost_->checkInfeasibilities(primalTolerance_);
918  }
919  if (nonLinearCost_->numberInfeasibilities()>0&&!progress->initialWeight_&&!ifValuesPass&&infeasibilityCost_==1.0e10) {
920    // first time infeasible - start up weight computation
921    double * oldDj = dj_;
922    double * oldCost = cost_;
923    int numberRows2 = numberRows_+numberExtraRows_;
924    int numberTotal = numberRows2+numberColumns_;
925    dj_ = new double[numberTotal];
926    cost_ = new double[numberTotal];
927    reducedCostWork_ = dj_;
928    rowReducedCost_ = dj_+numberColumns_;
929    objectiveWork_ = cost_;
930    rowObjectiveWork_ = cost_+numberColumns_;
931    double direction = optimizationDirection_*objectiveScale_;
932    const double * obj = objective();
933    memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
934    int iSequence;
935    if (columnScale_)
936      for (iSequence=0;iSequence<numberColumns_;iSequence++) 
937        cost_[iSequence] = obj[iSequence]*direction*columnScale_[iSequence];
938    else
939      for (iSequence=0;iSequence<numberColumns_;iSequence++) 
940        cost_[iSequence] = obj[iSequence]*direction;
941    computeDuals(NULL);
942    int numberSame=0;
943    int numberDifferent=0;
944    int numberZero=0;
945    int numberFreeSame=0;
946    int numberFreeDifferent=0;
947    int numberFreeZero=0;
948    int n=0;
949    for (iSequence=0;iSequence<numberTotal;iSequence++) {
950      if (getStatus(iSequence) != basic&&!flagged(iSequence)) {
951        // not basic
952        double distanceUp = upper_[iSequence]-solution_[iSequence];
953        double distanceDown = solution_[iSequence]-lower_[iSequence];
954        double feasibleDj = dj_[iSequence];
955        double infeasibleDj = oldDj[iSequence]-feasibleDj;
956        double value = feasibleDj*infeasibleDj;
957        if (distanceUp>primalTolerance_) {
958          // Check if "free"
959          if (distanceDown>primalTolerance_) {
960            // free
961            if (value>dualTolerance_) {
962              numberFreeSame++;
963            } else if(value<-dualTolerance_) {
964              numberFreeDifferent++;
965              dj_[n++] = feasibleDj/infeasibleDj;
966            } else {
967              numberFreeZero++;
968            }
969          } else {
970            // should not be negative
971            if (value>dualTolerance_) {
972              numberSame++;
973            } else if(value<-dualTolerance_) {
974              numberDifferent++;
975              dj_[n++] = feasibleDj/infeasibleDj;
976            } else {
977              numberZero++;
978            }
979          }
980        } else if (distanceDown>primalTolerance_) {
981          // should not be positive
982          if (value>dualTolerance_) {
983              numberSame++;
984            } else if(value<-dualTolerance_) {
985              numberDifferent++;
986              dj_[n++] = feasibleDj/infeasibleDj;
987            } else {
988              numberZero++;
989            }
990        }
991      }
992      progress->initialWeight_=-1.0;
993    }
994    //printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
995    //   numberSame,numberDifferent,numberZero,
996    //   numberFreeSame,numberFreeDifferent,numberFreeZero);
997    // we want most to be same
998    if (n) {
999      double most = 0.95;
1000      std::sort(dj_,dj_+n);
1001      int which = (int) ((1.0-most)*((double) n));
1002      double take = -dj_[which]*infeasibilityCost_;
1003      //printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
1004      take = -dj_[0]*infeasibilityCost_;
1005      infeasibilityCost_ = CoinMin(CoinMax(1000.0*take,1.0e8),1.0000001e10);;
1006      //printf("XXXX increasing weight to %g\n",infeasibilityCost_);
1007    }
1008    delete [] dj_;
1009    delete [] cost_;
1010    dj_= oldDj;
1011    cost_ = oldCost;
1012    reducedCostWork_ = dj_;
1013    rowReducedCost_ = dj_+numberColumns_;
1014    objectiveWork_ = cost_;
1015    rowObjectiveWork_ = cost_+numberColumns_;
1016    if (n)
1017      gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1018  }
1019  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
1020  if (!nonLinearCost_->numberInfeasibilities()&&infeasibilityCost_==1.0e10&&!ifValuesPass&&true) {
1021    // relax if default
1022    infeasibilityCost_ = CoinMin(CoinMax(100.0*sumDualInfeasibilities_,1.0e8),1.00000001e10);
1023    // reset looping criterion
1024    progress->reset();
1025    trueInfeasibility = 1.123456e10;
1026  }
1027  if (trueInfeasibility>1.0) {
1028    // If infeasibility going up may change weights
1029    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
1030    double lastInf = progress->lastInfeasibility(1);
1031    double lastInf3 = progress->lastInfeasibility(3);
1032    double thisObj = progress->lastObjective(0);
1033    double thisInf = progress->lastInfeasibility(0);
1034    thisObj += infeasibilityCost_*2.0*thisInf;
1035    double lastObj = progress->lastObjective(1);
1036    lastObj += infeasibilityCost_*2.0*lastInf;
1037    double lastObj3 = progress->lastObjective(3);
1038    lastObj3 += infeasibilityCost_*2.0*lastInf3;
1039    if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-7
1040        &&firstFree_<0) {
1041      if (handler_->logLevel()==63)
1042        printf("lastobj %g this %g force %d ",lastObj,thisObj,forceFactorization_);
1043      int maxFactor = factorization_->maximumPivots();
1044      if (maxFactor>10) {
1045        if (forceFactorization_<0)
1046          forceFactorization_= maxFactor;
1047        forceFactorization_ = CoinMax(1,(forceFactorization_>>2));
1048        if (handler_->logLevel()==63)
1049          printf("Reducing factorization frequency to %d\n",forceFactorization_);
1050      }
1051    } else if (lastObj3<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj3))-1.0e-7
1052        &&firstFree_<0) {
1053      if (handler_->logLevel()==63)
1054        printf("lastobj3 %g this3 %g `force %d ",lastObj3,thisObj,forceFactorization_);
1055      int maxFactor = factorization_->maximumPivots();
1056      if (maxFactor>10) {
1057        if (forceFactorization_<0)
1058          forceFactorization_= maxFactor;
1059        forceFactorization_ = CoinMax(1,(forceFactorization_*2)/3);
1060        if (handler_->logLevel()==63)
1061        printf("Reducing factorization frequency to %d\n",forceFactorization_);
1062      }
1063    } else if(lastInf<testValue||trueInfeasibility==1.123456e10) {
1064      if (infeasibilityCost_<1.0e14) {
1065        infeasibilityCost_ *= 1.5;
1066        // reset looping criterion
1067        progress->reset();
1068        if (handler_->logLevel()==63)
1069          printf("increasing weight to %g\n",infeasibilityCost_);
1070        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1071      }
1072    }
1073  }
1074  // we may wish to say it is optimal even if infeasible
1075  bool alwaysOptimal = (specialOptions_&1)!=0;
1076  // give code benefit of doubt
1077  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1078      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1079    // say optimal (with these bounds etc)
1080    numberDualInfeasibilities_ = 0;
1081    sumDualInfeasibilities_ = 0.0;
1082    numberPrimalInfeasibilities_ = 0;
1083    sumPrimalInfeasibilities_ = 0.0;
1084    // But check if in sprint
1085    if (originalModel) {
1086      // Carry on and re-do
1087      numberDualInfeasibilities_ = -776;
1088    }
1089    // But if real primal infeasibilities nonzero carry on
1090    if (nonLinearCost_->numberInfeasibilities()) {
1091      // most likely to happen if infeasible
1092      double relaxedToleranceP=primalTolerance_;
1093      // we can't really trust infeasibilities if there is primal error
1094      double error = CoinMin(1.0e-2,largestPrimalError_);
1095      // allow tolerance at least slightly bigger than standard
1096      relaxedToleranceP = relaxedToleranceP +  error;
1097      int ninfeas = nonLinearCost_->numberInfeasibilities();
1098      double sum = nonLinearCost_->sumInfeasibilities();
1099      double average = sum/ ((double) ninfeas);
1100#ifdef COIN_DEVELOP
1101      if (handler_->logLevel()>0)
1102        printf("nonLinearCost says infeasible %d summing to %g\n",
1103               ninfeas,sum);
1104#endif
1105      if (average>relaxedToleranceP) {
1106        sumOfRelaxedPrimalInfeasibilities_ = sum;
1107        numberPrimalInfeasibilities_ = ninfeas;
1108        sumPrimalInfeasibilities_ = sum;
1109#ifdef COIN_DEVELOP
1110        bool unflagged = 
1111#endif
1112        unflag();
1113#ifdef COIN_DEVELOP
1114        if (unflagged&&handler_->logLevel()>0)
1115          printf(" - but flagged variables\n");
1116#endif
1117      }
1118    }
1119  } 
1120  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
1121  if ((dualFeasible()||problemStatus_==-4)&&!ifValuesPass) {
1122    // see if extra helps
1123    if (nonLinearCost_->numberInfeasibilities()&&
1124         (nonLinearCost_->sumInfeasibilities()>1.0e-3||sumOfRelaxedPrimalInfeasibilities_)
1125        &&!alwaysOptimal) {
1126      //may need infeasiblity cost changed
1127      // we can see if we can construct a ray
1128      // make up a new objective
1129      double saveWeight = infeasibilityCost_;
1130      // save nonlinear cost as we are going to switch off costs
1131      ClpNonLinearCost * nonLinear = nonLinearCost_;
1132      // do twice to make sure Primal solution has settled
1133      // put non-basics to bounds in case tolerance moved
1134      // put back original costs
1135      createRim(4);
1136      nonLinearCost_->checkInfeasibilities(0.0);
1137      gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1138
1139      infeasibilityCost_=1.0e100;
1140      // put back original costs
1141      createRim(4);
1142      nonLinearCost_->checkInfeasibilities(primalTolerance_);
1143      // may have fixed infeasibilities - double check
1144      if (nonLinearCost_->numberInfeasibilities()==0) {
1145        // carry on
1146        problemStatus_ = -1;
1147        infeasibilityCost_=saveWeight;
1148        nonLinearCost_->checkInfeasibilities(primalTolerance_);
1149      } else {
1150        nonLinearCost_=NULL;
1151        // scale
1152        int i;
1153        for (i=0;i<numberRows_+numberColumns_;i++) 
1154          cost_[i] *= 1.0e-95;
1155        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1156        nonLinearCost_=nonLinear;
1157        infeasibilityCost_=saveWeight;
1158        if ((infeasibilityCost_>=1.0e18||
1159             numberDualInfeasibilities_==0)&&perturbation_==101) {
1160          goToDual=unPerturb(); // stop any further perturbation
1161          if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
1162            goToDual=false;
1163          nonLinearCost_->checkInfeasibilities(primalTolerance_);
1164          numberDualInfeasibilities_=1; // carry on
1165          problemStatus_=-1;
1166        }
1167        if (infeasibilityCost_>=1.0e20||
1168            numberDualInfeasibilities_==0) {
1169          // we are infeasible - use as ray
1170          delete [] ray_;
1171          ray_ = new double [numberRows_];
1172          CoinMemcpyN(dual_,numberRows_,ray_);
1173          // and get feasible duals
1174          infeasibilityCost_=0.0;
1175          createRim(4);
1176          nonLinearCost_->checkInfeasibilities(primalTolerance_);
1177          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1178          // so will exit
1179          infeasibilityCost_=1.0e30;
1180          // reset infeasibilities
1181          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
1182          numberPrimalInfeasibilities_=
1183            nonLinearCost_->numberInfeasibilities();
1184        }
1185        if (infeasibilityCost_<1.0e20) {
1186          infeasibilityCost_ *= 5.0;
1187          // reset looping criterion
1188          progress->reset();
1189          changeMade_++; // say change made
1190          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
1191            <<infeasibilityCost_
1192            <<CoinMessageEol;
1193          // put back original costs and then check
1194          createRim(4);
1195          nonLinearCost_->checkInfeasibilities(0.0);
1196          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1197          problemStatus_=-1; //continue
1198          goToDual=false;
1199        } else {
1200          // say infeasible
1201          problemStatus_ = 1;
1202        }
1203      }
1204    } else {
1205      // may be optimal
1206      if (perturbation_==101) {
1207        goToDual=unPerturb(); // stop any further perturbation
1208        if (numberRows_>20000&&!numberTimesOptimal_)
1209          goToDual=0; // Better to carry on a bit longer
1210        lastCleaned=-1; // carry on
1211      }
1212      bool unflagged = (unflag()!=0);
1213      if ( lastCleaned!=numberIterations_||unflagged) {
1214        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
1215          <<primalTolerance_
1216          <<CoinMessageEol;
1217        if (numberTimesOptimal_<4) {
1218          numberTimesOptimal_++;
1219          changeMade_++; // say change made
1220          if (numberTimesOptimal_==1) {
1221            // better to have small tolerance even if slower
1222            factorization_->zeroTolerance(1.0e-15);
1223          }
1224          lastCleaned=numberIterations_;
1225          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
1226            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
1227              <<CoinMessageEol;
1228          double oldTolerance = primalTolerance_;
1229          primalTolerance_=dblParam_[ClpPrimalTolerance];
1230#if 0
1231          double * xcost = new double[numberRows_+numberColumns_];
1232          double * xlower = new double[numberRows_+numberColumns_];
1233          double * xupper = new double[numberRows_+numberColumns_];
1234          double * xdj = new double[numberRows_+numberColumns_];
1235          double * xsolution = new double[numberRows_+numberColumns_];
1236   CoinMemcpyN(cost_,(numberRows_+numberColumns_),xcost);
1237   CoinMemcpyN(lower_,(numberRows_+numberColumns_),xlower);
1238   CoinMemcpyN(upper_,(numberRows_+numberColumns_),xupper);
1239   CoinMemcpyN(dj_,(numberRows_+numberColumns_),xdj);
1240   CoinMemcpyN(solution_,(numberRows_+numberColumns_),xsolution);
1241#endif
1242          // put back original costs and then check
1243          createRim(4);
1244          nonLinearCost_->checkInfeasibilities(oldTolerance);
1245#if 0
1246          int i;
1247          for (i=0;i<numberRows_+numberColumns_;i++) {
1248            if (cost_[i]!=xcost[i])
1249              printf("** %d old cost %g new %g sol %g\n",
1250                     i,xcost[i],cost_[i],solution_[i]);
1251            if (lower_[i]!=xlower[i])
1252              printf("** %d old lower %g new %g sol %g\n",
1253                     i,xlower[i],lower_[i],solution_[i]);
1254            if (upper_[i]!=xupper[i])
1255              printf("** %d old upper %g new %g sol %g\n",
1256                     i,xupper[i],upper_[i],solution_[i]);
1257            if (dj_[i]!=xdj[i])
1258              printf("** %d old dj %g new %g sol %g\n",
1259                     i,xdj[i],dj_[i],solution_[i]);
1260            if (solution_[i]!=xsolution[i])
1261              printf("** %d old solution %g new %g sol %g\n",
1262                     i,xsolution[i],solution_[i],solution_[i]);
1263          }
1264          delete [] xcost;
1265          delete [] xupper;
1266          delete [] xlower;
1267          delete [] xdj;
1268          delete [] xsolution;
1269#endif
1270          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1271          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1272              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1273            // say optimal (with these bounds etc)
1274            numberDualInfeasibilities_ = 0;
1275            sumDualInfeasibilities_ = 0.0;
1276            numberPrimalInfeasibilities_ = 0;
1277            sumPrimalInfeasibilities_ = 0.0;
1278          }
1279          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
1280            problemStatus_=0;
1281          else
1282            problemStatus_ = -1;
1283        } else {
1284          problemStatus_=0; // optimal
1285          if (lastCleaned<numberIterations_) {
1286            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1287              <<CoinMessageEol;
1288          }
1289        }
1290      } else {
1291        if (!alwaysOptimal||!sumOfRelaxedPrimalInfeasibilities_)
1292          problemStatus_=0; // optimal
1293        else
1294          problemStatus_=1; // infeasible
1295      }
1296    }
1297  } else {
1298    // see if looks unbounded
1299    if (problemStatus_==-5) {
1300      if (nonLinearCost_->numberInfeasibilities()) {
1301        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
1302          // back off weight
1303          infeasibilityCost_ = 1.0e13;
1304          // reset looping criterion
1305          progress->reset();
1306          unPerturb(); // stop any further perturbation
1307        }
1308        //we need infeasiblity cost changed
1309        if (infeasibilityCost_<1.0e20) {
1310          infeasibilityCost_ *= 5.0;
1311          // reset looping criterion
1312          progress->reset();
1313          changeMade_++; // say change made
1314          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
1315            <<infeasibilityCost_
1316            <<CoinMessageEol;
1317          // put back original costs and then check
1318          createRim(4);
1319          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1320          problemStatus_=-1; //continue
1321        } else {
1322          // say infeasible
1323          problemStatus_ = 1;
1324          // we are infeasible - use as ray
1325          delete [] ray_;
1326          ray_ = new double [numberRows_];
1327          CoinMemcpyN(dual_,numberRows_,ray_);
1328        }
1329      } else {
1330        // say unbounded
1331        problemStatus_ = 2;
1332      } 
1333    } else {
1334      // carry on
1335      problemStatus_ = -1;
1336      if(type==3&&problemStatus_!=-5) {
1337        //bool unflagged =
1338        unflag();
1339        if (sumDualInfeasibilities_<1.0e-3||
1340            (sumDualInfeasibilities_/(double) numberDualInfeasibilities_)<1.0e-5||
1341            progress->lastIterationNumber(0)==numberIterations_) {
1342          if (!numberPrimalInfeasibilities_) {
1343            if (numberTimesOptimal_<4) {
1344              numberTimesOptimal_++;
1345              changeMade_++; // say change made
1346            } else {
1347              problemStatus_=0;
1348              secondaryStatus_=5;
1349            }
1350          }
1351        }
1352      }
1353    }
1354  }
1355  if (problemStatus_==0) {
1356    double objVal = nonLinearCost_->feasibleCost();
1357    double tol = 1.0e-10*CoinMax(fabs(objVal),fabs(objectiveValue_))+1.0e-8;
1358    if (fabs(objVal-objectiveValue_)>tol) {
1359#ifdef COIN_DEVELOP
1360      if (handler_->logLevel()>0)
1361        printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
1362               objVal,objectiveValue_);
1363#endif
1364      objectiveValue_ = objVal;
1365    }
1366  }
1367  // save extra stuff
1368  matrix_->generalExpanded(this,5,dummy);
1369  if (type==0||type==1) {
1370    if (type!=1||!saveStatus_) {
1371      // create save arrays
1372      delete [] saveStatus_;
1373      delete [] savedSolution_;
1374      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1375      savedSolution_ = new double [numberRows_+numberColumns_];
1376    }
1377    // save arrays
1378    CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
1379    CoinMemcpyN(rowActivityWork_,
1380           numberRows_,savedSolution_+numberColumns_);
1381    CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
1382  }
1383  // see if in Cbc etc
1384  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
1385  bool disaster=false;
1386  if (disasterArea_&&inCbcOrOther&&disasterArea_->check()) {
1387    disasterArea_->saveInfo();
1388    disaster=true;
1389  }
1390  if (disaster)
1391    problemStatus_=3;
1392  if (problemStatus_<0&&!changeMade_) {
1393    problemStatus_=4; // unknown
1394  }
1395  lastGoodIteration_ = numberIterations_;
1396  if (numberIterations_>lastBadIteration_+100)
1397    moreSpecialOptions_ &= ~16; // clear check accuracy flag
1398  if (goToDual) 
1399    problemStatus_=10; // try dual
1400  // make sure first free monotonic
1401  if (firstFree_>=0&&saveFirstFree>=0) {
1402    firstFree_=saveFirstFree;
1403    nextSuperBasic(1,NULL);
1404  }
1405  if (doFactorization) {
1406    // restore weights (if saved) - also recompute infeasibility list
1407    if (tentativeStatus>-3) 
1408      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
1409    else
1410      primalColumnPivot_->saveWeights(this,3);
1411    if (saveThreshold) {
1412      // use default at present
1413      factorization_->sparseThreshold(0);
1414      factorization_->goSparse();
1415    }
1416  }
1417  // Allow matrices to be sorted etc
1418  int fake=-999; // signal sort
1419  matrix_->correctSequence(this,fake,fake);
1420}
1421/*
1422   Row array has pivot column
1423   This chooses pivot row.
1424   For speed, we may need to go to a bucket approach when many
1425   variables go through bounds
1426   On exit rhsArray will have changes in costs of basic variables
1427*/
1428void 
1429ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
1430                            CoinIndexedVector * rhsArray,
1431                            CoinIndexedVector * spareArray,
1432                            CoinIndexedVector * spareArray2,
1433                            int valuesPass)
1434{
1435  double saveDj = dualIn_;
1436  if (valuesPass&&objective_->type()<2) {
1437    dualIn_ = cost_[sequenceIn_];
1438
1439    double * work=rowArray->denseVector();
1440    int number=rowArray->getNumElements();
1441    int * which=rowArray->getIndices();
1442
1443    int iIndex;
1444    for (iIndex=0;iIndex<number;iIndex++) {
1445     
1446      int iRow = which[iIndex];
1447      double alpha = work[iIndex];
1448      int iPivot=pivotVariable_[iRow];
1449      dualIn_ -= alpha*cost_[iPivot];
1450    }
1451    // determine direction here
1452    if (dualIn_<-dualTolerance_) {
1453      directionIn_=1;
1454    } else if (dualIn_>dualTolerance_) {
1455      directionIn_=-1;
1456    } else {
1457      // towards nearest bound
1458      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
1459        directionIn_=-1;
1460        dualIn_=dualTolerance_;
1461      } else {
1462        directionIn_=1;
1463        dualIn_=-dualTolerance_;
1464      }
1465    }
1466  }
1467
1468  // sequence stays as row number until end
1469  pivotRow_=-1;
1470  int numberRemaining=0;
1471
1472  double totalThru=0.0; // for when variables flip
1473  // Allow first few iterations to take tiny
1474  double acceptablePivot=1.0e-1*acceptablePivot_;
1475  if (numberIterations_>100)
1476    acceptablePivot=acceptablePivot_;
1477  if (factorization_->pivots()>10)
1478    acceptablePivot=1.0e+3*acceptablePivot_; // if we have iterated be more strict
1479  else if (factorization_->pivots()>5)
1480    acceptablePivot=1.0e+2*acceptablePivot_; // if we have iterated be slightly more strict
1481  else if (factorization_->pivots())
1482    acceptablePivot=acceptablePivot_; // relax
1483  double bestEverPivot=acceptablePivot;
1484  int lastPivotRow = -1;
1485  double lastPivot=0.0;
1486  double lastTheta=1.0e50;
1487
1488  // use spareArrays to put ones looked at in
1489  // First one is list of candidates
1490  // We could compress if we really know we won't need any more
1491  // Second array has current set of pivot candidates
1492  // with a backup list saved in double * part of indexed vector
1493
1494  // pivot elements
1495  double * spare;
1496  // indices
1497  int * index;
1498  spareArray->clear();
1499  spare = spareArray->denseVector();
1500  index = spareArray->getIndices();
1501
1502  // we also need somewhere for effective rhs
1503  double * rhs=rhsArray->denseVector();
1504  // and we can use indices to point to alpha
1505  // that way we can store fabs(alpha)
1506  int * indexPoint = rhsArray->getIndices();
1507  //int numberFlip=0; // Those which may change if flips
1508
1509  /*
1510    First we get a list of possible pivots.  We can also see if the
1511    problem looks unbounded.
1512
1513    At first we increase theta and see what happens.  We start
1514    theta at a reasonable guess.  If in right area then we do bit by bit.
1515    We save possible pivot candidates
1516
1517   */
1518
1519  // do first pass to get possibles
1520  // We can also see if unbounded
1521
1522  double * work=rowArray->denseVector();
1523  int number=rowArray->getNumElements();
1524  int * which=rowArray->getIndices();
1525
1526  // we need to swap sign if coming in from ub
1527  double way = directionIn_;
1528  double maximumMovement;
1529  if (way>0.0) 
1530    maximumMovement = CoinMin(1.0e30,upperIn_-valueIn_);
1531  else
1532    maximumMovement = CoinMin(1.0e30,valueIn_-lowerIn_);
1533
1534  double averageTheta = nonLinearCost_->averageTheta();
1535  double tentativeTheta = CoinMin(10.0*averageTheta,maximumMovement);
1536  double upperTheta = maximumMovement;
1537  if (tentativeTheta>0.5*maximumMovement)
1538    tentativeTheta=maximumMovement;
1539  bool thetaAtMaximum=tentativeTheta==maximumMovement;
1540  // In case tiny bounds increase
1541  if (maximumMovement<1.0)
1542    tentativeTheta *= 1.1;
1543  double dualCheck = fabs(dualIn_);
1544  // but make a bit more pessimistic
1545  dualCheck=CoinMax(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1546
1547  int iIndex;
1548  int pivotOne=-1;
1549  //#define CLP_DEBUG
1550#ifdef CLP_DEBUG
1551  if (numberIterations_==-3839||numberIterations_==-3840) {
1552    double dj=cost_[sequenceIn_];
1553    printf("cost in on %d is %g, dual in %g\n",sequenceIn_,dj,dualIn_);
1554    for (iIndex=0;iIndex<number;iIndex++) {
1555
1556      int iRow = which[iIndex];
1557      double alpha = work[iIndex];
1558      int iPivot=pivotVariable_[iRow];
1559      dj -= alpha*cost_[iPivot];
1560      printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1561             iRow,iPivot,lower_[iPivot],solution_[iPivot],upper_[iPivot],
1562             alpha, solution_[iPivot]-1.0e9*alpha,cost_[iPivot],dj);
1563    }
1564  }
1565#endif
1566  while (true) {
1567    pivotOne=-1;
1568    totalThru=0.0;
1569    // We also re-compute reduced cost
1570    numberRemaining=0;
1571    dualIn_ = cost_[sequenceIn_];
1572#ifndef NDEBUG
1573    double tolerance = primalTolerance_*1.002;
1574#endif
1575    for (iIndex=0;iIndex<number;iIndex++) {
1576
1577      int iRow = which[iIndex];
1578      double alpha = work[iIndex];
1579      int iPivot=pivotVariable_[iRow];
1580      if (cost_[iPivot])
1581        dualIn_ -= alpha*cost_[iPivot];
1582      alpha *= way;
1583      double oldValue = solution_[iPivot];
1584      // get where in bound sequence
1585      // note that after this alpha is actually fabs(alpha)
1586      bool possible;
1587      // do computation same way as later on in primal
1588      if (alpha>0.0) {
1589        // basic variable going towards lower bound
1590        double bound = lower_[iPivot];
1591        // must be exactly same as when used
1592        double change = tentativeTheta*alpha;
1593        possible = (oldValue-change)<=bound+primalTolerance_;
1594        oldValue -= bound;
1595      } else {
1596        // basic variable going towards upper bound
1597        double bound = upper_[iPivot];
1598        // must be exactly same as when used
1599        double change = tentativeTheta*alpha;
1600        possible = (oldValue-change)>=bound-primalTolerance_;
1601        oldValue = bound-oldValue;
1602        alpha = - alpha;
1603      }
1604      double value;
1605      assert (oldValue>=-tolerance);
1606      if (possible) {
1607        value=oldValue-upperTheta*alpha;
1608        if (value<-primalTolerance_&&alpha>=acceptablePivot) {
1609          upperTheta = (oldValue+primalTolerance_)/alpha;
1610          pivotOne=numberRemaining;
1611        }
1612        // add to list
1613        spare[numberRemaining]=alpha;
1614        rhs[numberRemaining]=oldValue;
1615        indexPoint[numberRemaining]=iIndex;
1616        index[numberRemaining++]=iRow;
1617        totalThru += alpha;
1618        setActive(iRow);
1619        //} else if (value<primalTolerance_*1.002) {
1620        // May change if is a flip
1621        //indexRhs[numberFlip++]=iRow;
1622      }
1623    }
1624    if (upperTheta<maximumMovement&&totalThru*infeasibilityCost_>=1.0001*dualCheck) {
1625      // Can pivot here
1626      break;
1627    } else if (!thetaAtMaximum) {
1628      //printf("Going round with average theta of %g\n",averageTheta);
1629      tentativeTheta=maximumMovement;
1630      thetaAtMaximum=true; // seems to be odd compiler error
1631    } else {
1632      break;
1633    }
1634  }
1635  totalThru=0.0;
1636
1637  theta_=maximumMovement;
1638
1639  bool goBackOne = false;
1640  if (objective_->type()>1) 
1641    dualIn_=saveDj;
1642
1643  //printf("%d remain out of %d\n",numberRemaining,number);
1644  int iTry=0;
1645#define MAXTRY 1000
1646  if (numberRemaining&&upperTheta<maximumMovement) {
1647    // First check if previously chosen one will work
1648    if (pivotOne>=0&&0) {
1649      double thruCost = infeasibilityCost_*spare[pivotOne];
1650      if (thruCost>=0.99*fabs(dualIn_))
1651        printf("Could pivot on %d as change %g dj %g\n",
1652               index[pivotOne],thruCost,dualIn_);
1653      double alpha = spare[pivotOne];
1654      double oldValue = rhs[pivotOne];
1655      theta_ = oldValue/alpha;
1656      pivotRow_=pivotOne;
1657      // Stop loop
1658      iTry=MAXTRY;
1659    }
1660
1661    // first get ratio with tolerance
1662    for ( ;iTry<MAXTRY;iTry++) {
1663     
1664      upperTheta=maximumMovement;
1665      int iBest=-1;
1666      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1667
1668        double alpha = spare[iIndex];
1669        double oldValue = rhs[iIndex];
1670        double value = oldValue-upperTheta*alpha;
1671
1672        if (value<-primalTolerance_ && alpha>=acceptablePivot) {
1673          upperTheta = (oldValue+primalTolerance_)/alpha;
1674          iBest=iIndex; // just in case weird numbers
1675        }
1676      }
1677     
1678      // now look at best in this lot
1679      // But also see how infeasible small pivots will make
1680      double sumInfeasibilities=0.0;
1681      double bestPivot=acceptablePivot;
1682      pivotRow_=-1;
1683      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1684
1685        int iRow = index[iIndex];
1686        double alpha = spare[iIndex];
1687        double oldValue = rhs[iIndex];
1688        double value = oldValue-upperTheta*alpha;
1689
1690        if (value<=0||iBest==iIndex) {
1691          // how much would it cost to go thru and modify bound
1692          double trueAlpha=way*work[indexPoint[iIndex]];
1693          totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],trueAlpha,rhs[iIndex]);
1694          setActive(iRow);
1695          if (alpha>bestPivot) {
1696            bestPivot=alpha;
1697            theta_ = oldValue/bestPivot;
1698            pivotRow_=iIndex;
1699          } else if (alpha<acceptablePivot) {
1700            if (value<-primalTolerance_)
1701              sumInfeasibilities += -value-primalTolerance_;
1702          }
1703        }
1704      }
1705      if (bestPivot<0.1*bestEverPivot&&
1706          bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
1707        // back to previous one
1708        goBackOne = true;
1709        break;
1710      } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1711        if (lastPivot>acceptablePivot) {
1712          // back to previous one
1713          goBackOne = true;
1714        } else {
1715          // can only get here if all pivots so far too small
1716        }
1717        break;
1718      } else if (totalThru>=dualCheck) {
1719        if (sumInfeasibilities>primalTolerance_&&!nonLinearCost_->numberInfeasibilities()) {
1720          // Looks a bad choice
1721          if (lastPivot>acceptablePivot) {
1722            goBackOne=true;
1723          } else {
1724            // say no good
1725            dualIn_=0.0;
1726          }
1727        } 
1728        break; // no point trying another loop
1729      } else {
1730        lastPivotRow=pivotRow_;
1731        lastTheta = theta_;
1732        if (bestPivot>bestEverPivot)
1733          bestEverPivot=bestPivot;
1734      }   
1735    }
1736    // can get here without pivotRow_ set but with lastPivotRow
1737    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1738      // back to previous one
1739      pivotRow_=lastPivotRow;
1740      theta_ = lastTheta;
1741    }
1742  } else if (pivotRow_<0&&maximumMovement>1.0e20) {
1743    // looks unbounded
1744    valueOut_=COIN_DBL_MAX; // say odd
1745    if (nonLinearCost_->numberInfeasibilities()) {
1746      // but infeasible??
1747      // move variable but don't pivot
1748      tentativeTheta=1.0e50;
1749      for (iIndex=0;iIndex<number;iIndex++) {
1750        int iRow = which[iIndex];
1751        double alpha = work[iIndex];
1752        int iPivot=pivotVariable_[iRow];
1753        alpha *= way;
1754        double oldValue = solution_[iPivot];
1755        // get where in bound sequence
1756        // note that after this alpha is actually fabs(alpha)
1757        if (alpha>0.0) {
1758          // basic variable going towards lower bound
1759          double bound = lower_[iPivot];
1760          oldValue -= bound;
1761        } else {
1762          // basic variable going towards upper bound
1763          double bound = upper_[iPivot];
1764          oldValue = bound-oldValue;
1765          alpha = - alpha;
1766        }
1767        if (oldValue-tentativeTheta*alpha<0.0) {
1768          tentativeTheta = oldValue/alpha;
1769        }
1770      }
1771      // If free in then see if we can get to 0.0
1772      if (lowerIn_<-1.0e20&&upperIn_>1.0e20) {
1773        if (dualIn_*valueIn_>0.0) {
1774          if (fabs(valueIn_)<1.0e-2&&(tentativeTheta<fabs(valueIn_)||tentativeTheta>1.0e20)) {
1775            tentativeTheta = fabs(valueIn_);
1776          }
1777        }
1778      }
1779      if (tentativeTheta<1.0e10)
1780        valueOut_=valueIn_+way*tentativeTheta;
1781    }
1782  }
1783  //if (iTry>50)
1784  //printf("** %d tries\n",iTry);
1785  if (pivotRow_>=0) {
1786    int position=pivotRow_; // position in list
1787    pivotRow_=index[position];
1788    alpha_=work[indexPoint[position]];
1789    // translate to sequence
1790    sequenceOut_ = pivotVariable_[pivotRow_];
1791    valueOut_ = solution(sequenceOut_);
1792    lowerOut_=lower_[sequenceOut_];
1793    upperOut_=upper_[sequenceOut_];
1794#define MINIMUMTHETA 1.0e-12
1795    // Movement should be minimum for anti-degeneracy - unless
1796    // fixed variable out
1797    double minimumTheta;
1798    if (upperOut_>lowerOut_)
1799      minimumTheta=MINIMUMTHETA;
1800    else
1801      minimumTheta=0.0;
1802    // But can't go infeasible
1803    double distance;
1804    if (alpha_*way>0.0) 
1805      distance=valueOut_-lowerOut_;
1806    else
1807      distance=upperOut_-valueOut_;
1808    if (distance-minimumTheta*fabs(alpha_)<-primalTolerance_)
1809      minimumTheta = CoinMax(0.0,(distance+0.5*primalTolerance_)/fabs(alpha_));
1810    // will we need to increase tolerance
1811    //#define CLP_DEBUG
1812    double largestInfeasibility = primalTolerance_;
1813    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
1814      theta_=minimumTheta;
1815      for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
1816        largestInfeasibility = CoinMax(largestInfeasibility,
1817                                    -(rhs[iIndex]-spare[iIndex]*theta_));
1818      }
1819//#define CLP_DEBUG
1820#ifdef CLP_DEBUG
1821      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
1822        printf("Primal tolerance increased from %g to %g\n",
1823               primalTolerance_,largestInfeasibility);
1824#endif
1825//#undef CLP_DEBUG
1826      primalTolerance_ = CoinMax(primalTolerance_,largestInfeasibility);
1827    }
1828    // Need to look at all in some cases
1829    if (theta_>tentativeTheta) {
1830      for (iIndex=0;iIndex<number;iIndex++) 
1831        setActive(which[iIndex]);
1832    }
1833    if (way<0.0) 
1834      theta_ = - theta_;
1835    double newValue = valueOut_ - theta_*alpha_;
1836    // If  4 bit set - Force outgoing variables to exact bound (primal)
1837    if (alpha_*way<0.0) {
1838      directionOut_=-1;      // to upper bound
1839      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1840        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1841      } else {
1842        upperOut_ = newValue;
1843      }
1844    } else {
1845      directionOut_=1;      // to lower bound
1846      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1847        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1848      } else {
1849        lowerOut_ = newValue;
1850      }
1851    }
1852    dualOut_ = reducedCost(sequenceOut_);
1853  } else if (maximumMovement<1.0e20) {
1854    // flip
1855    pivotRow_ = -2; // so we can tell its a flip
1856    sequenceOut_ = sequenceIn_;
1857    valueOut_ = valueIn_;
1858    dualOut_ = dualIn_;
1859    lowerOut_ = lowerIn_;
1860    upperOut_ = upperIn_;
1861    alpha_ = 0.0;
1862    if (way<0.0) {
1863      directionOut_=1;      // to lower bound
1864      theta_ = lowerOut_ - valueOut_;
1865    } else {
1866      directionOut_=-1;      // to upper bound
1867      theta_ = upperOut_ - valueOut_;
1868    }
1869  }
1870
1871  double theta1 = CoinMax(theta_,1.0e-12);
1872  double theta2 = numberIterations_*nonLinearCost_->averageTheta();
1873  // Set average theta
1874  nonLinearCost_->setAverageTheta((theta1+theta2)/((double) (numberIterations_+1)));
1875  //if (numberIterations_%1000==0)
1876  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
1877
1878  // clear arrays
1879
1880  CoinZeroN(spare,numberRemaining);
1881
1882  // put back original bounds etc
1883  CoinMemcpyN(index,numberRemaining,rhsArray->getIndices());
1884  rhsArray->setNumElements(numberRemaining);
1885  rhsArray->setPacked();
1886  nonLinearCost_->goBackAll(rhsArray);
1887  rhsArray->clear();
1888
1889}
1890/*
1891   Chooses primal pivot column
1892   updateArray has cost updates (also use pivotRow_ from last iteration)
1893   Would be faster with separate region to scan
1894   and will have this (with square of infeasibility) when steepest
1895   For easy problems we can just choose one of the first columns we look at
1896*/
1897void 
1898ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
1899                               CoinIndexedVector * spareRow1,
1900                               CoinIndexedVector * spareRow2,
1901                               CoinIndexedVector * spareColumn1,
1902                               CoinIndexedVector * spareColumn2)
1903{
1904  sequenceIn_ = primalColumnPivot_->pivotColumn(updates,spareRow1,
1905                                               spareRow2,spareColumn1,
1906                                               spareColumn2);
1907  if (sequenceIn_>=0) {
1908    valueIn_=solution_[sequenceIn_];
1909    dualIn_=dj_[sequenceIn_];
1910    if (nonLinearCost_->lookBothWays()) {
1911      // double check
1912      ClpSimplex::Status status = getStatus(sequenceIn_);
1913     
1914      switch(status) {
1915      case ClpSimplex::atUpperBound:
1916        if (dualIn_<0.0) {
1917          // move to other side
1918          printf("For %d U (%g, %g, %g) dj changed from %g",
1919                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1920                 upper_[sequenceIn_],dualIn_);
1921          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
1922          printf(" to %g\n",dualIn_);
1923          nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
1924          setStatus(sequenceIn_,ClpSimplex::atLowerBound);
1925        }
1926        break;
1927      case ClpSimplex::atLowerBound:
1928        if (dualIn_>0.0) {
1929          // move to other side
1930          printf("For %d L (%g, %g, %g) dj changed from %g",
1931                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1932                 upper_[sequenceIn_],dualIn_);
1933          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
1934          printf(" to %g\n",dualIn_);
1935          nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
1936          setStatus(sequenceIn_,ClpSimplex::atUpperBound);
1937        }
1938        break;
1939      default:
1940        break;
1941      }
1942    }
1943    lowerIn_=lower_[sequenceIn_];
1944    upperIn_=upper_[sequenceIn_];
1945    if (dualIn_>0.0)
1946      directionIn_ = -1;
1947    else 
1948      directionIn_ = 1;
1949  } else {
1950    sequenceIn_ = -1;
1951  }
1952}
1953/* The primals are updated by the given array.
1954   Returns number of infeasibilities.
1955   After rowArray will have list of cost changes
1956*/
1957int 
1958ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
1959                                        double theta,
1960                                        double & objectiveChange,
1961                                        int valuesPass)
1962{
1963  // Cost on pivot row may change - may need to change dualIn
1964  double oldCost=0.0;
1965  if (pivotRow_>=0)
1966    oldCost = cost_[sequenceOut_];
1967  //rowArray->scanAndPack();
1968  double * work=rowArray->denseVector();
1969  int number=rowArray->getNumElements();
1970  int * which=rowArray->getIndices();
1971
1972  int newNumber = 0;
1973  int pivotPosition = -1;
1974  nonLinearCost_->setChangeInCost(0.0);
1975  //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
1976  //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
1977  // allow for case where bound+tolerance == bound
1978  //double tolerance = 0.999999*primalTolerance_;
1979  double relaxedTolerance = 1.001*primalTolerance_;
1980  int iIndex;
1981  if (!valuesPass) {
1982    for (iIndex=0;iIndex<number;iIndex++) {
1983     
1984      int iRow = which[iIndex];
1985      double alpha = work[iIndex];
1986      work[iIndex]=0.0;
1987      int iPivot=pivotVariable_[iRow];
1988      double change = theta*alpha;
1989      double value = solution_[iPivot] - change;
1990      solution_[iPivot]=value;
1991#ifndef NDEBUG
1992      // check if not active then okay
1993      if (!active(iRow)&&(specialOptions_&4)==0&&pivotRow_!=-1) {
1994        // But make sure one going out is feasible
1995        if (change>0.0) {
1996          // going down
1997          if (value<=lower_[iPivot]+primalTolerance_) {
1998            if (iPivot==sequenceOut_&&value>lower_[iPivot]-relaxedTolerance)
1999              value=lower_[iPivot];
2000            double difference = nonLinearCost_->setOne(iPivot,value);
2001            assert (!difference||fabs(change)>1.0e9);
2002          }
2003        } else {
2004          // going up
2005          if (value>=upper_[iPivot]-primalTolerance_) {
2006            if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
2007              value=upper_[iPivot];
2008            double difference = nonLinearCost_->setOne(iPivot,value);
2009            assert (!difference||fabs(change)>1.0e9);
2010          }
2011        }
2012      }
2013#endif   
2014      if (active(iRow)||theta_<0.0) {
2015        clearActive(iRow);
2016        // But make sure one going out is feasible
2017        if (change>0.0) {
2018          // going down
2019          if (value<=lower_[iPivot]+primalTolerance_) {
2020            if (iPivot==sequenceOut_&&value>=lower_[iPivot]-relaxedTolerance)
2021              value=lower_[iPivot];
2022            double difference = nonLinearCost_->setOne(iPivot,value);
2023            if (difference) {
2024              if (iRow==pivotRow_)
2025                pivotPosition=newNumber;
2026              work[newNumber] = difference;
2027              //change reduced cost on this
2028              dj_[iPivot] = -difference;
2029              which[newNumber++]=iRow;
2030            }
2031          }
2032        } else {
2033          // going up
2034          if (value>=upper_[iPivot]-primalTolerance_) {
2035            if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
2036              value=upper_[iPivot];
2037            double difference = nonLinearCost_->setOne(iPivot,value);
2038            if (difference) {
2039              if (iRow==pivotRow_)
2040                pivotPosition=newNumber;
2041              work[newNumber] = difference;
2042              //change reduced cost on this
2043              dj_[iPivot] = -difference;
2044              which[newNumber++]=iRow;
2045            }
2046          }
2047        }
2048      }
2049    }
2050  } else {
2051    // values pass so look at all
2052    for (iIndex=0;iIndex<number;iIndex++) {
2053     
2054      int iRow = which[iIndex];
2055      double alpha = work[iIndex];
2056      work[iIndex]=0.0;
2057      int iPivot=pivotVariable_[iRow];
2058      double change = theta*alpha;
2059      double value = solution_[iPivot] - change;
2060      solution_[iPivot]=value;
2061      clearActive(iRow);
2062      // But make sure one going out is feasible
2063      if (change>0.0) {
2064        // going down
2065        if (value<=lower_[iPivot]+primalTolerance_) {
2066          if (iPivot==sequenceOut_&&value>lower_[iPivot]-relaxedTolerance)
2067            value=lower_[iPivot];
2068          double difference = nonLinearCost_->setOne(iPivot,value);
2069          if (difference) {
2070            if (iRow==pivotRow_)
2071              pivotPosition=newNumber;
2072            work[newNumber] = difference;
2073            //change reduced cost on this
2074            dj_[iPivot] = -difference;
2075            which[newNumber++]=iRow;
2076          }
2077        }
2078      } else {
2079        // going up
2080        if (value>=upper_[iPivot]-primalTolerance_) {
2081          if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
2082            value=upper_[iPivot];
2083          double difference = nonLinearCost_->setOne(iPivot,value);
2084          if (difference) {
2085            if (iRow==pivotRow_)
2086              pivotPosition=newNumber;
2087            work[newNumber] = difference;
2088            //change reduced cost on this
2089            dj_[iPivot] = -difference;
2090            which[newNumber++]=iRow;
2091          }
2092        }
2093      }
2094    }
2095  }
2096  objectiveChange += nonLinearCost_->changeInCost();
2097  rowArray->setPacked();
2098#if 0
2099  rowArray->setNumElements(newNumber);
2100  rowArray->expand();
2101  if (pivotRow_>=0) {
2102    dualIn_ += (oldCost-cost_[sequenceOut_]);
2103    // update change vector to include pivot
2104    rowArray->add(pivotRow_,-dualIn_);
2105    // and convert to packed
2106    rowArray->scanAndPack();
2107  } else {
2108    // and convert to packed
2109    rowArray->scanAndPack();
2110  }
2111#else
2112  if (pivotRow_>=0) {
2113    double dualIn = dualIn_+(oldCost-cost_[sequenceOut_]);
2114    // update change vector to include pivot
2115    if (pivotPosition>=0) {
2116      work[pivotPosition] -= dualIn;
2117    } else {
2118      work[newNumber]=-dualIn;
2119      which[newNumber++]=pivotRow_;
2120    }
2121  }
2122  rowArray->setNumElements(newNumber);
2123#endif
2124  return 0;
2125}
2126// Perturbs problem
2127void 
2128ClpSimplexPrimal::perturb(int type)
2129{
2130  if (perturbation_>100)
2131    return; //perturbed already
2132  if (perturbation_==100)
2133    perturbation_=50; // treat as normal
2134  int savePerturbation = perturbation_;
2135  int i;
2136  if (!numberIterations_)
2137    cleanStatus(); // make sure status okay
2138  // Make sure feasible bounds
2139  if (nonLinearCost_)
2140    nonLinearCost_->feasibleBounds();
2141  // look at element range
2142  double smallestNegative;
2143  double largestNegative;
2144  double smallestPositive;
2145  double largestPositive;
2146  matrix_->rangeOfElements(smallestNegative, largestNegative,
2147                           smallestPositive, largestPositive);
2148  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2149  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2150  double elementRatio = largestPositive/smallestPositive;
2151  if (!numberIterations_&&perturbation_==50) {
2152    // See if we need to perturb
2153    double * sort = new double[numberRows_];
2154    for (i=0;i<numberRows_;i++) {
2155      double lo = fabs(rowLower_[i]);
2156      double up = fabs(rowUpper_[i]);
2157      double value=0.0;
2158      if (lo&&lo<1.0e20) {
2159        if (up&&up<1.0e20)
2160          value = 0.5*(lo+up);
2161        else
2162          value=lo;
2163      } else {
2164        if (up&&up<1.0e20)
2165          value = up;
2166      }
2167      sort[i]=value;
2168    }
2169    std::sort(sort,sort+numberRows_);
2170    int number=1;
2171    double last = sort[0];
2172    for (i=1;i<numberRows_;i++) {
2173      if (last!=sort[i])
2174        number++;
2175      last=sort[i];
2176    }
2177    delete [] sort;
2178    //printf("ratio number diff rhs %g, element ratio %g\n",((double)number)/((double) numberRows_),
2179    //                                                                elementRatio);
2180    if (number*3>numberRows_||elementRatio>1.0e12) {
2181      perturbation_=100;
2182      return; // good enough
2183    }
2184  }
2185  // primal perturbation
2186  double perturbation=1.0e-20;
2187  int numberNonZero=0;
2188  // maximum fraction of rhs/bounds to perturb
2189  double maximumFraction = 1.0e-5;
2190  if (perturbation_>=50) {
2191    perturbation = 1.0e-4;
2192    for (i=0;i<numberColumns_+numberRows_;i++) {
2193      if (upper_[i]>lower_[i]+primalTolerance_) {
2194        double lowerValue, upperValue;
2195        if (lower_[i]>-1.0e20)
2196          lowerValue = fabs(lower_[i]);
2197        else
2198          lowerValue=0.0;
2199        if (upper_[i]<1.0e20)
2200          upperValue = fabs(upper_[i]);
2201        else
2202          upperValue=0.0;
2203        double value = CoinMax(fabs(lowerValue),fabs(upperValue));
2204        value = CoinMin(value,upper_[i]-lower_[i]);
2205#if 1
2206        if (value) {
2207          perturbation += value;
2208          numberNonZero++;
2209        }
2210#else
2211        perturbation = CoinMax(perturbation,value);
2212#endif
2213      }
2214    }
2215    if (numberNonZero) 
2216      perturbation /= (double) numberNonZero;
2217    else
2218      perturbation = 1.0e-1;
2219  } else if (perturbation_<100) {
2220    perturbation = pow(10.0,perturbation_);
2221    // user is in charge
2222    maximumFraction = 1.0;
2223  }
2224  double largestZero=0.0;
2225  double largest=0.0;
2226  double largestPerCent=0.0;
2227  bool printOut=(handler_->logLevel()==63);
2228  printOut=false; //off
2229  // Check if all slack
2230  int number=0;
2231  int iSequence;
2232  for (iSequence=0;iSequence<numberRows_;iSequence++) {
2233    if (getRowStatus(iSequence)==basic) 
2234      number++;
2235  }
2236  if (rhsScale_>100.0) {
2237    // tone down perturbation
2238    maximumFraction *= 0.1;
2239  }
2240  if (number!=numberRows_)
2241    type=1;
2242  // modify bounds
2243  // Change so at least 1.0e-5 and no more than 0.1
2244  // For now just no more than 0.1
2245  // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
2246  if (type==1) {
2247    double tolerance = 100.0*primalTolerance_;
2248    //double multiplier = perturbation*maximumFraction;
2249    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2250      if (getStatus(iSequence)==basic) {
2251        double lowerValue = lower_[iSequence];
2252        double upperValue = upper_[iSequence];
2253        if (upperValue>lowerValue+tolerance) {
2254          double solutionValue = solution_[iSequence];
2255          double difference = upperValue-lowerValue;
2256          difference = CoinMin(difference,perturbation);
2257          difference = CoinMin(difference,fabs(solutionValue)+1.0);
2258          double value = maximumFraction*(difference+1.0);
2259          value = CoinMin(value,0.1);
2260          value *= randomNumberGenerator_.randomDouble();
2261          if (solutionValue-lowerValue<=primalTolerance_) {
2262            lower_[iSequence] -= value;
2263          } else if (upperValue-solutionValue<=primalTolerance_) {
2264            upper_[iSequence] += value;
2265          } else {
2266#if 0
2267            if (iSequence>=numberColumns_) {
2268              // may not be at bound - but still perturb (unless free)
2269              if (upperValue>1.0e30&&lowerValue<-1.0e30)
2270                value=0.0;
2271              else
2272                value = - value; // as -1.0 in matrix
2273            } else {
2274              value = 0.0;
2275            }
2276#else
2277            value=0.0;
2278#endif
2279          }
2280          if (value) {
2281            if (printOut)
2282              printf("col %d lower from %g to %g, upper from %g to %g\n",
2283                     iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
2284            if (solutionValue) {
2285              largest = CoinMax(largest,value);
2286              if (value>(fabs(solutionValue)+1.0)*largestPerCent)
2287                largestPerCent=value/(fabs(solutionValue)+1.0);
2288            } else {
2289              largestZero = CoinMax(largestZero,value);
2290            } 
2291          }
2292        }
2293      }
2294    }
2295  } else {
2296    double tolerance = 100.0*primalTolerance_;
2297    for (i=0;i<numberColumns_;i++) {
2298      double lowerValue=lower_[i], upperValue=upper_[i];
2299      if (upperValue>lowerValue+primalTolerance_) {
2300        double value = perturbation*maximumFraction;
2301        value = CoinMin(value,0.1);
2302        value *= randomNumberGenerator_.randomDouble();
2303        if (savePerturbation!=50) {
2304          if (fabs(value)<=primalTolerance_)
2305            value=0.0;
2306          if (lowerValue>-1.0e20&&lowerValue)
2307            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2308          if (upperValue<1.0e20&&upperValue)
2309            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2310        } else if (value) {
2311          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2312          // get in range
2313          if (valueL<=tolerance) {
2314            valueL *= 10.0;
2315            while (valueL<=tolerance) 
2316              valueL *= 10.0;
2317          } else if (valueL>1.0) {
2318            valueL *= 0.1;
2319            while (valueL>1.0) 
2320              valueL *= 0.1;
2321          }
2322          if (lowerValue>-1.0e20&&lowerValue)
2323            lowerValue -= valueL;
2324          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2325          // get in range
2326          if (valueU<=tolerance) {
2327            valueU *= 10.0;
2328            while (valueU<=tolerance) 
2329              valueU *= 10.0;
2330          } else if (valueU>1.0) {
2331            valueU *= 0.1;
2332            while (valueU>1.0) 
2333              valueU *= 0.1;
2334          }
2335          if (upperValue<1.0e20&&upperValue)
2336            upperValue += valueU;
2337        }
2338        if (lowerValue!=lower_[i]) {
2339          double difference = fabs(lowerValue-lower_[i]);
2340          largest = CoinMax(largest,difference);
2341          if (difference>fabs(lower_[i])*largestPerCent)
2342            largestPerCent=fabs(difference/lower_[i]);
2343        } 
2344        if (upperValue!=upper_[i]) {
2345          double difference = fabs(upperValue-upper_[i]);
2346          largest = CoinMax(largest,difference);
2347          if (difference>fabs(upper_[i])*largestPerCent)
2348            largestPerCent=fabs(difference/upper_[i]);
2349        } 
2350        if (printOut)
2351          printf("col %d lower from %g to %g, upper from %g to %g\n",
2352                 i,lower_[i],lowerValue,upper_[i],upperValue);
2353      }
2354      lower_[i]=lowerValue;
2355      upper_[i]=upperValue;
2356    }
2357    for (;i<numberColumns_+numberRows_;i++) {
2358      double lowerValue=lower_[i], upperValue=upper_[i];
2359      double value = perturbation*maximumFraction;
2360      value = CoinMin(value,0.1);
2361      value *= randomNumberGenerator_.randomDouble();
2362      if (upperValue>lowerValue+tolerance) {
2363        if (savePerturbation!=50) {
2364          if (fabs(value)<=primalTolerance_)
2365            value=0.0;
2366          if (lowerValue>-1.0e20&&lowerValue)
2367            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2368          if (upperValue<1.0e20&&upperValue)
2369            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2370        } else if (value) {
2371          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2372          // get in range
2373          if (valueL<=tolerance) {
2374            valueL *= 10.0;
2375            while (valueL<=tolerance) 
2376              valueL *= 10.0;
2377          } else if (valueL>1.0) {
2378            valueL *= 0.1;
2379            while (valueL>1.0) 
2380              valueL *= 0.1;
2381          }
2382          if (lowerValue>-1.0e20&&lowerValue)
2383            lowerValue -= valueL;
2384          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2385          // get in range
2386          if (valueU<=tolerance) {
2387            valueU *= 10.0;
2388            while (valueU<=tolerance) 
2389              valueU *= 10.0;
2390          } else if (valueU>1.0) {
2391            valueU *= 0.1;
2392            while (valueU>1.0) 
2393              valueU *= 0.1;
2394          }
2395          if (upperValue<1.0e20&&upperValue)
2396            upperValue += valueU;
2397        }
2398      } else if (upperValue>0.0) {
2399        upperValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2400        lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2401      } else if (upperValue<0.0) {
2402        upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2403        lowerValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2404      } else {
2405      }
2406      if (lowerValue!=lower_[i]) {
2407        double difference = fabs(lowerValue-lower_[i]);
2408        largest = CoinMax(largest,difference);
2409        if (difference>fabs(lower_[i])*largestPerCent)
2410          largestPerCent=fabs(difference/lower_[i]);
2411      } 
2412      if (upperValue!=upper_[i]) {
2413        double difference = fabs(upperValue-upper_[i]);
2414        largest = CoinMax(largest,difference);
2415        if (difference>fabs(upper_[i])*largestPerCent)
2416          largestPerCent=fabs(difference/upper_[i]);
2417      } 
2418      if (printOut)
2419        printf("row %d lower from %g to %g, upper from %g to %g\n",
2420               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
2421      lower_[i]=lowerValue;
2422      upper_[i]=upperValue;
2423    }
2424  }
2425  // Clean up
2426  for (i=0;i<numberColumns_+numberRows_;i++) {
2427    switch(getStatus(i)) {
2428     
2429    case basic:
2430      break;
2431    case atUpperBound:
2432      solution_[i]=upper_[i];
2433      break;
2434    case isFixed:
2435    case atLowerBound:
2436      solution_[i]=lower_[i];
2437      break;
2438    case isFree:
2439      break;
2440    case superBasic:
2441      break;
2442    }
2443  }
2444  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2445    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
2446    <<CoinMessageEol;
2447  // redo nonlinear costs
2448  // say perturbed
2449  perturbation_=101;
2450}
2451// un perturb
2452bool
2453ClpSimplexPrimal::unPerturb()
2454{
2455  if (perturbation_!=101)
2456    return false;
2457  // put back original bounds and costs
2458  createRim(1+4);
2459  sanityCheck();
2460  // unflag
2461  unflag();
2462  // get a valid nonlinear cost function
2463  delete nonLinearCost_;
2464  nonLinearCost_= new ClpNonLinearCost(this);
2465  perturbation_ = 102; // stop any further perturbation
2466  // move non basic variables to new bounds
2467  nonLinearCost_->checkInfeasibilities(0.0);
2468#if 1
2469  // Try using dual
2470  return true;
2471#else
2472  gutsOfSolution(NULL,NULL,ifValuesPass!=0);
2473  return false;
2474#endif
2475 
2476}
2477// Unflag all variables and return number unflagged
2478int 
2479ClpSimplexPrimal::unflag()
2480{
2481  int i;
2482  int number = numberRows_+numberColumns_;
2483  int numberFlagged=0;
2484  // we can't really trust infeasibilities if there is dual error
2485  // allow tolerance bigger than standard to check on duals
2486  double relaxedToleranceD=dualTolerance_ + CoinMin(1.0e-2,10.0*largestDualError_);
2487  for (i=0;i<number;i++) {
2488    if (flagged(i)) {
2489      clearFlagged(i);
2490      // only say if reasonable dj
2491      if (fabs(dj_[i])>relaxedToleranceD)
2492        numberFlagged++;
2493    }
2494  }
2495  numberFlagged += matrix_->generalExpanded(this,8,i);
2496  if (handler_->logLevel()>2&&numberFlagged&&objective_->type()>1)
2497    printf("%d unflagged\n",numberFlagged);
2498  return numberFlagged;
2499}
2500// Do not change infeasibility cost and always say optimal
2501void 
2502ClpSimplexPrimal::alwaysOptimal(bool onOff)
2503{
2504  if (onOff)
2505    specialOptions_ |= 1;
2506  else
2507    specialOptions_ &= ~1;
2508}
2509bool 
2510ClpSimplexPrimal::alwaysOptimal() const
2511{
2512  return (specialOptions_&1)!=0;
2513}
2514// Flatten outgoing variables i.e. - always to exact bound
2515void 
2516ClpSimplexPrimal::exactOutgoing(bool onOff)
2517{
2518  if (onOff)
2519    specialOptions_ |= 4;
2520  else
2521    specialOptions_ &= ~4;
2522}
2523bool 
2524ClpSimplexPrimal::exactOutgoing() const
2525{
2526  return (specialOptions_&4)!=0;
2527}
2528/*
2529  Reasons to come out (normal mode/user mode):
2530  -1 normal
2531  -2 factorize now - good iteration/ NA
2532  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
2533  -4 inaccuracy - refactorize - no iteration/ NA
2534  -5 something flagged - go round again/ pivot not possible
2535  +2 looks unbounded
2536  +3 max iterations (iteration done)
2537*/
2538int
2539ClpSimplexPrimal::pivotResult(int ifValuesPass)
2540{
2541
2542  bool roundAgain=true;
2543  int returnCode=-1;
2544
2545  // loop round if user setting and doing refactorization
2546  while (roundAgain) {
2547    roundAgain=false;
2548    returnCode=-1;
2549    pivotRow_=-1;
2550    sequenceOut_=-1;
2551    rowArray_[1]->clear();
2552#if 0
2553    {
2554      int seq[]={612,643};
2555      int k;
2556      for (k=0;k<sizeof(seq)/sizeof(int);k++) {
2557        int iSeq=seq[k];
2558        if (getColumnStatus(iSeq)!=basic) {
2559          double djval;
2560          double * work;
2561          int number;
2562          int * which;
2563         
2564          int iIndex;
2565          unpack(rowArray_[1],iSeq);
2566          factorization_->updateColumn(rowArray_[2],rowArray_[1]);
2567          djval = cost_[iSeq];
2568          work=rowArray_[1]->denseVector();
2569          number=rowArray_[1]->getNumElements();
2570          which=rowArray_[1]->getIndices();
2571         
2572          for (iIndex=0;iIndex<number;iIndex++) {
2573           
2574            int iRow = which[iIndex];
2575            double alpha = work[iRow];
2576            int iPivot=pivotVariable_[iRow];
2577            djval -= alpha*cost_[iPivot];
2578          }
2579          double comp = 1.0e-8 + 1.0e-7*(CoinMax(fabs(dj_[iSeq]),fabs(djval)));
2580          if (fabs(djval-dj_[iSeq])>comp)
2581            printf("Bad dj %g for %d - true is %g\n",
2582                   dj_[iSeq],iSeq,djval);
2583          assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
2584          rowArray_[1]->clear();
2585        }
2586      }
2587    }
2588#endif
2589
2590    // we found a pivot column
2591    // update the incoming column
2592    unpackPacked(rowArray_[1]);
2593    // save reduced cost
2594    double saveDj = dualIn_;
2595    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2596    // Get extra rows
2597    matrix_->extendUpdated(this,rowArray_[1],0);
2598    // do ratio test and re-compute dj
2599    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2600              ifValuesPass);
2601    if (ifValuesPass) {
2602      saveDj=dualIn_;
2603      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
2604      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2605        if(fabs(dualIn_)<1.0e2*dualTolerance_&&objective_->type()<2) {
2606          // try other way
2607          directionIn_=-directionIn_;
2608          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2609                    0);
2610        }
2611        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2612          if (solveType_==1) {
2613            // reject it
2614            char x = isColumn(sequenceIn_) ? 'C' :'R';
2615            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2616              <<x<<sequenceWithin(sequenceIn_)
2617              <<CoinMessageEol;
2618            setFlagged(sequenceIn_);
2619            progress_.clearBadTimes();
2620            lastBadIteration_ = numberIterations_; // say be more cautious
2621            clearAll();
2622            pivotRow_=-1;
2623          }
2624          returnCode=-5;
2625          break;
2626        }
2627      }
2628    }
2629    // need to clear toIndex_ in gub
2630    // ? when can I clear stuff
2631    // Clean up any gub stuff
2632    matrix_->extendUpdated(this,rowArray_[1],1);
2633    double checkValue=1.0e-2;
2634    if (largestDualError_>1.0e-5)
2635      checkValue=1.0e-1;
2636    if (!ifValuesPass&&solveType_==1&&(saveDj*dualIn_<1.0e-20||
2637        fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj))||
2638                        fabs(dualIn_)<dualTolerance_)) {
2639      char x = isColumn(sequenceIn_) ? 'C' :'R';
2640      handler_->message(CLP_PRIMAL_DJ,messages_)
2641        <<x<<sequenceWithin(sequenceIn_)
2642        <<saveDj<<dualIn_
2643        <<CoinMessageEol;
2644      if(lastGoodIteration_ != numberIterations_) {
2645        clearAll();
2646        pivotRow_=-1; // say no weights update
2647        returnCode=-4;
2648        if(lastGoodIteration_+1 == numberIterations_) {
2649          // not looking wonderful - try cleaning bounds
2650          // put non-basics to bounds in case tolerance moved
2651          nonLinearCost_->checkInfeasibilities(0.0);
2652        }
2653        sequenceOut_=-1;
2654        break;
2655      } else {
2656        // take on more relaxed criterion
2657        if (saveDj*dualIn_<1.0e-20||
2658            fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))||
2659            fabs(dualIn_)<dualTolerance_) {
2660          // need to reject something
2661          char x = isColumn(sequenceIn_) ? 'C' :'R';
2662          handler_->message(CLP_SIMPLEX_FLAG,messages_)
2663            <<x<<sequenceWithin(sequenceIn_)
2664            <<CoinMessageEol;
2665          setFlagged(sequenceIn_);
2666          progress_.clearBadTimes();
2667          lastBadIteration_ = numberIterations_; // say be more cautious
2668          clearAll();
2669          pivotRow_=-1;
2670          returnCode=-5;
2671          sequenceOut_=-1;
2672          break;
2673        }
2674      }
2675    }
2676    if (pivotRow_>=0) {
2677      if (solveType_==2) {
2678        // **** Coding for user interface
2679        // do ray
2680        primalRay(rowArray_[1]);
2681        // update duals
2682        // as packed need to find pivot row
2683        //assert (rowArray_[1]->packedMode());
2684        //int i;
2685
2686        //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2687        CoinAssert (fabs(alpha_)>1.0e-8);
2688        double multiplier = dualIn_/alpha_;
2689        rowArray_[0]->insert(pivotRow_,multiplier);
2690        factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
2691        // put row of tableau in rowArray[0] and columnArray[0]
2692        matrix_->transposeTimes(this,-1.0,
2693                                rowArray_[0],columnArray_[1],columnArray_[0]);
2694        // update column djs
2695        int i;
2696        int * index = columnArray_[0]->getIndices();
2697        int number = columnArray_[0]->getNumElements();
2698        double * element = columnArray_[0]->denseVector();
2699        for (i=0;i<number;i++) {
2700          int ii = index[i];
2701          dj_[ii] += element[ii];
2702          reducedCost_[ii] = dj_[ii];
2703          element[ii]=0.0;
2704        }
2705        columnArray_[0]->setNumElements(0);
2706        // and row djs
2707        index = rowArray_[0]->getIndices();
2708        number = rowArray_[0]->getNumElements();
2709        element = rowArray_[0]->denseVector();
2710        for (i=0;i<number;i++) {
2711          int ii = index[i];
2712          dj_[ii+numberColumns_] += element[ii];
2713          dual_[ii] = dj_[ii+numberColumns_];
2714          element[ii]=0.0;
2715        }
2716        rowArray_[0]->setNumElements(0);
2717        // check incoming
2718        CoinAssert (fabs(dj_[sequenceIn_])<1.0e-1);
2719      }
2720      // if stable replace in basis
2721      // If gub or odd then alpha and pivotRow may change
2722      int updateType=0;
2723      int updateStatus = matrix_->generalExpanded(this,3,updateType);
2724      if (updateType>=0)
2725        updateStatus = factorization_->replaceColumn(this,
2726                                                     rowArray_[2],
2727                                                     rowArray_[1],
2728                                                     pivotRow_,
2729                                                     alpha_,
2730                                                     (moreSpecialOptions_&16)!=0);
2731
2732      // if no pivots, bad update but reasonable alpha - take and invert
2733      if (updateStatus==2&&
2734          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2735        updateStatus=4;
2736      if (updateStatus==1||updateStatus==4) {
2737        // slight error
2738        if (factorization_->pivots()>5||updateStatus==4) {
2739          returnCode=-3;
2740        }
2741      } else if (updateStatus==2) {
2742        // major error
2743        // better to have small tolerance even if slower
2744        factorization_->zeroTolerance(1.0e-15);
2745        int maxFactor = factorization_->maximumPivots();
2746        if (maxFactor>10) {
2747          if (forceFactorization_<0)
2748            forceFactorization_= maxFactor;
2749          forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2750        } 
2751        // later we may need to unwind more e.g. fake bounds
2752        if(lastGoodIteration_ != numberIterations_) {
2753          clearAll();
2754          pivotRow_=-1;
2755          if (solveType_==1) {
2756            returnCode=-4;
2757            break;
2758          } else {
2759            // user in charge - re-factorize
2760            int lastCleaned=0;
2761            ClpSimplexProgress dummyProgress;
2762            if (saveStatus_)
2763              statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2764            else
2765              statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2766            roundAgain=true;
2767            continue;
2768          }
2769        } else {
2770          // need to reject something
2771          if (solveType_==1) {
2772            char x = isColumn(sequenceIn_) ? 'C' :'R';
2773            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2774              <<x<<sequenceWithin(sequenceIn_)
2775              <<CoinMessageEol;
2776            setFlagged(sequenceIn_);
2777            progress_.clearBadTimes();
2778          }
2779          lastBadIteration_ = numberIterations_; // say be more cautious
2780          clearAll();
2781          pivotRow_=-1;
2782          sequenceOut_=-1;
2783          returnCode = -5;
2784          break;
2785
2786        }
2787      } else if (updateStatus==3) {
2788        // out of memory
2789        // increase space if not many iterations
2790        if (factorization_->pivots()<
2791            0.5*factorization_->maximumPivots()&&
2792            factorization_->pivots()<200)
2793          factorization_->areaFactor(
2794                                     factorization_->areaFactor() * 1.1);
2795        returnCode =-2; // factorize now
2796      } else if (updateStatus==5) {
2797        problemStatus_=-2; // factorize now
2798      }
2799      // here do part of steepest - ready for next iteration
2800      if (!ifValuesPass)
2801        primalColumnPivot_->updateWeights(rowArray_[1]);
2802    } else {
2803      if (pivotRow_==-1) {
2804        // no outgoing row is valid
2805        if (valueOut_!=COIN_DBL_MAX) {
2806          double objectiveChange=0.0;
2807          theta_=valueOut_-valueIn_;
2808          updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2809          solution_[sequenceIn_] += theta_;
2810        }
2811        rowArray_[0]->clear();
2812        if (!factorization_->pivots()&&acceptablePivot_<=1.0e-8) {
2813          returnCode = 2; //say looks unbounded
2814          // do ray
2815          primalRay(rowArray_[1]);
2816        } else if (solveType_==2) {
2817          // refactorize
2818          int lastCleaned=0;
2819          ClpSimplexProgress dummyProgress;
2820          if (saveStatus_)
2821            statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2822          else
2823            statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2824          roundAgain=true;
2825          continue;
2826        } else {
2827          acceptablePivot_=1.0e-8;
2828          returnCode = 4; //say looks unbounded but has iterated
2829        }
2830        break;
2831      } else {
2832        // flipping from bound to bound
2833      }
2834    }
2835
2836    double oldCost = 0.0;
2837    if (sequenceOut_>=0)
2838      oldCost=cost_[sequenceOut_];
2839    // update primal solution
2840   
2841    double objectiveChange=0.0;
2842    // after this rowArray_[1] is not empty - used to update djs
2843    // If pivot row >= numberRows then may be gub
2844    int savePivot = pivotRow_;
2845    if (pivotRow_>=numberRows_)
2846      pivotRow_=-1;
2847    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2848    pivotRow_=savePivot;
2849   
2850    double oldValue = valueIn_;
2851    if (directionIn_==-1) {
2852      // as if from upper bound
2853      if (sequenceIn_!=sequenceOut_) {
2854        // variable becoming basic
2855        valueIn_ -= fabs(theta_);
2856      } else {
2857        valueIn_=lowerIn_;
2858      }
2859    } else {
2860      // as if from lower bound
2861      if (sequenceIn_!=sequenceOut_) {
2862        // variable becoming basic
2863        valueIn_ += fabs(theta_);
2864      } else {
2865        valueIn_=upperIn_;
2866      }
2867    }
2868    objectiveChange += dualIn_*(valueIn_-oldValue);
2869    // outgoing
2870    if (sequenceIn_!=sequenceOut_) {
2871      if (directionOut_>0) {
2872        valueOut_ = lowerOut_;
2873      } else {
2874        valueOut_ = upperOut_;
2875      }
2876      if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2877        valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2878      else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2879        valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2880      // may not be exactly at bound and bounds may have changed
2881      // Make sure outgoing looks feasible
2882      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2883      // May have got inaccurate
2884      //if (oldCost!=cost_[sequenceOut_])
2885      //printf("costchange on %d from %g to %g\n",sequenceOut_,
2886      //       oldCost,cost_[sequenceOut_]);
2887      if (solveType_!=2)
2888        dj_[sequenceOut_]=cost_[sequenceOut_]-oldCost; // normally updated next iteration
2889      solution_[sequenceOut_]=valueOut_;
2890    }
2891    // change cost and bounds on incoming if primal
2892    nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2893    int whatNext=housekeeping(objectiveChange);
2894    //nonLinearCost_->validate();
2895#if CLP_DEBUG >1
2896    {
2897      double sum;
2898      int ninf= matrix_->checkFeasible(this,sum);
2899      if (ninf)
2900        printf("infeas %d\n",ninf);
2901    }
2902#endif
2903    if (whatNext==1) {
2904        returnCode =-2; // refactorize
2905    } else if (whatNext==2) {
2906      // maximum iterations or equivalent
2907      returnCode=3;
2908    } else if(numberIterations_ == lastGoodIteration_
2909              + 2 * factorization_->maximumPivots()) {
2910      // done a lot of flips - be safe
2911      returnCode =-2; // refactorize
2912    }
2913    // Check event
2914    {
2915      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2916      if (status>=0) {
2917        problemStatus_=5;
2918        secondaryStatus_=ClpEventHandler::endOfIteration;
2919        returnCode=3;
2920      }
2921    }
2922  }
2923  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
2924    // refactorize here
2925    int lastCleaned=0;
2926    ClpSimplexProgress dummyProgress;
2927    if (saveStatus_)
2928      statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2929    else
2930      statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2931    if (problemStatus_==5) {
2932      printf("Singular basis\n");
2933      problemStatus_=-1;
2934      returnCode=5;
2935    }
2936  }
2937#ifdef CLP_DEBUG
2938  {
2939    int i;
2940    // not [1] as may have information
2941    for (i=0;i<4;i++) {
2942      if (i!=1)
2943        rowArray_[i]->checkClear();
2944    }   
2945    for (i=0;i<2;i++) {
2946      columnArray_[i]->checkClear();
2947    }   
2948  }     
2949#endif
2950  return returnCode;
2951}
2952// Create primal ray
2953void 
2954ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
2955{
2956  delete [] ray_;
2957  ray_ = new double [numberColumns_];
2958  CoinZeroN(ray_,numberColumns_);
2959  int number=rowArray->getNumElements();
2960  int * index = rowArray->getIndices();
2961  double * array = rowArray->denseVector();
2962  double way=-directionIn_;
2963  int i;
2964  double zeroTolerance=1.0e-12;
2965  if (sequenceIn_<numberColumns_)
2966    ray_[sequenceIn_]=directionIn_;
2967  if (!rowArray->packedMode()) {
2968    for (i=0;i<number;i++) {
2969      int iRow=index[i];
2970      int iPivot=pivotVariable_[iRow];
2971      double arrayValue = array[iRow];
2972      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2973        ray_[iPivot] = way* arrayValue;
2974    }
2975  } else {
2976    for (i=0;i<number;i++) {
2977      int iRow=index[i];
2978      int iPivot=pivotVariable_[iRow];
2979      double arrayValue = array[i];
2980      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2981        ray_[iPivot] = way* arrayValue;
2982    }
2983  }
2984}
2985/* Get next superbasic -1 if none,
2986   Normal type is 1
2987   If type is 3 then initializes sorted list
2988   if 2 uses list.
2989*/
2990int 
2991ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
2992{
2993  if (firstFree_>=0&&superBasicType) {
2994    int returnValue=-1;
2995    bool finished=false;
2996    while (!finished) {
2997      returnValue=firstFree_;
2998      int iColumn=firstFree_+1;
2999      if (superBasicType>1) {
3000        if (superBasicType>2) {
3001          // Initialize list
3002          // Wild guess that lower bound more natural than upper
3003          int number=0;
3004          double * work=columnArray->denseVector();
3005          int * which=columnArray->getIndices();
3006          for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
3007            if (!flagged(iColumn)) {
3008              if (getStatus(iColumn)==superBasic) {
3009                if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
3010                  solution_[iColumn]=lower_[iColumn];
3011                  setStatus(iColumn,atLowerBound);
3012                } else if (fabs(solution_[iColumn]-upper_[iColumn])
3013                           <=primalTolerance_) {
3014                  solution_[iColumn]=upper_[iColumn];
3015                  setStatus(iColumn,atUpperBound);
3016                } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
3017                  setStatus(iColumn,isFree);
3018                  break;
3019                } else if (!flagged(iColumn)) {
3020                  // put ones near bounds at end after sorting
3021                  work[number]= - CoinMin(0.1*(solution_[iColumn]-lower_[iColumn]),
3022                                      upper_[iColumn]-solution_[iColumn]);
3023                  which[number++] = iColumn;
3024                }
3025              }
3026            }
3027          }
3028          CoinSort_2(work,work+number,which);
3029          columnArray->setNumElements(number);
3030          CoinZeroN(work,number);
3031        }
3032        int * which=columnArray->getIndices();
3033        int number = columnArray->getNumElements();
3034        if (!number) {
3035          // finished
3036          iColumn = numberRows_+numberColumns_;
3037          returnValue=-1;
3038        } else {
3039          number--;
3040          returnValue=which[number];
3041          iColumn=returnValue;
3042          columnArray->setNumElements(number);
3043        }     
3044      } else {
3045        for (;iColumn<numberRows_+numberColumns_;iColumn++) {
3046          if (!flagged(iColumn)) {
3047            if (getStatus(iColumn)==superBasic) {
3048              if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
3049                solution_[iColumn]=lower_[iColumn];
3050                setStatus(iColumn,atLowerBound);
3051              } else if (fabs(solution_[iColumn]-upper_[iColumn])
3052                         <=primalTolerance_) {
3053                solution_[iColumn]=upper_[iColumn];
3054                setStatus(iColumn,atUpperBound);
3055              } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
3056                setStatus(iColumn,isFree);
3057                break;
3058              } else {
3059                break;
3060              }
3061            }
3062          }
3063        }
3064      }
3065      firstFree_ = iColumn;
3066      finished=true;
3067      if (firstFree_==numberRows_+numberColumns_)
3068        firstFree_=-1;
3069      if (returnValue>=0&&getStatus(returnValue)!=superBasic&&getStatus(returnValue)!=isFree)
3070        finished=false; // somehow picked up odd one
3071    }
3072    return returnValue;
3073  } else {
3074    return -1;
3075  }
3076}
3077void
3078ClpSimplexPrimal::clearAll()
3079{
3080  // Clean up any gub stuff
3081  matrix_->extendUpdated(this,rowArray_[1],1);
3082  int number=rowArray_[1]->getNumElements();
3083  int * which=rowArray_[1]->getIndices();
3084 
3085  int iIndex;
3086  for (iIndex=0;iIndex<number;iIndex++) {
3087   
3088    int iRow = which[iIndex];
3089    clearActive(iRow);
3090  }
3091  rowArray_[1]->clear();
3092  // make sure any gub sets are clean
3093  matrix_->generalExpanded(this,11,sequenceIn_);
3094 
3095}
3096// Sort of lexicographic resolve
3097int
3098ClpSimplexPrimal::lexSolve()
3099{
3100  algorithm_ = +1;
3101  //specialOptions_ |= 4;
3102
3103  // save data
3104  ClpDataSave data = saveData();
3105  matrix_->refresh(this); // make sure matrix okay
3106 
3107  // Save so can see if doing after dual
3108  int initialStatus=problemStatus_;
3109  int initialIterations = numberIterations_;
3110  int initialNegDjs=-1;
3111  // initialize - maybe values pass and algorithm_ is +1
3112  int ifValuesPass=0;
3113#if 0
3114  // if so - put in any superbasic costed slacks
3115  if (ifValuesPass&&specialOptions_<0x01000000) {
3116    // Get column copy
3117    const CoinPackedMatrix * columnCopy = matrix();
3118    const int * row = columnCopy->getIndices();
3119    const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3120    const int * columnLength = columnCopy->getVectorLengths();
3121    //const double * element = columnCopy->getElements();
3122    int n=0;
3123    for (int iColumn = 0;iColumn<numberColumns_;iColumn++) {
3124      if (columnLength[iColumn]==1) {
3125        Status status = getColumnStatus(iColumn);
3126        if (status!=basic&&status!=isFree) {
3127          double value = columnActivity_[iColumn];
3128          if (fabs(value-columnLower_[iColumn])>primalTolerance_&&
3129              fabs(value-columnUpper_[iColumn])>primalTolerance_) {
3130            int iRow = row[columnStart[iColumn]];
3131            if (getRowStatus(iRow)==basic) {
3132              setRowStatus(iRow,superBasic);
3133              setColumnStatus(iColumn,basic);
3134              n++;
3135            }
3136          }
3137        }   
3138      }
3139    }
3140    printf("%d costed slacks put in basis\n",n);
3141  }
3142#endif
3143  double * originalCost = NULL;
3144  double * originalLower = NULL;
3145  double * originalUpper = NULL;
3146  if (!startup(0,0)) {
3147   
3148    // Set average theta
3149    nonLinearCost_->setAverageTheta(1.0e3);
3150    int lastCleaned=0; // last time objective or bounds cleaned up
3151   
3152    // Say no pivot has occurred (for steepest edge and updates)
3153    pivotRow_=-2;
3154   
3155    // This says whether to restore things etc
3156    int factorType=0;
3157    if (problemStatus_<0&&perturbation_<100) {
3158      perturb(0);
3159      // Can't get here if values pass
3160      assert (!ifValuesPass);
3161      gutsOfSolution(NULL,NULL);
3162      if (handler_->logLevel()>2) {
3163        handler_->message(CLP_SIMPLEX_STATUS,messages_)
3164          <<numberIterations_<<objectiveValue();
3165        handler_->printing(sumPrimalInfeasibilities_>0.0)
3166          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
3167        handler_->printing(sumDualInfeasibilities_>0.0)
3168          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
3169        handler_->printing(numberDualInfeasibilitiesWithoutFree_
3170                           <numberDualInfeasibilities_)
3171                             <<numberDualInfeasibilitiesWithoutFree_;
3172        handler_->message()<<CoinMessageEol;
3173      }
3174    }
3175    ClpSimplex * saveModel=NULL;
3176    int stopSprint=-1;
3177    int sprintPass=0;
3178    int reasonableSprintIteration=0;
3179    int lastSprintIteration=0;
3180    double lastObjectiveValue=COIN_DBL_MAX;
3181    // Start check for cycles
3182    progress_.fillFromModel(this);
3183    progress_.startCheck();
3184    /*
3185      Status of problem:
3186      0 - optimal
3187      1 - infeasible
3188      2 - unbounded
3189      -1 - iterating
3190      -2 - factorization wanted
3191      -3 - redo checking without factorization
3192      -4 - looks infeasible
3193      -5 - looks unbounded
3194    */
3195    originalCost = CoinCopyOfArray(cost_,numberColumns_+numberRows_);
3196    originalLower = CoinCopyOfArray(lower_,numberColumns_+numberRows_);
3197    originalUpper = CoinCopyOfArray(upper_,numberColumns_+numberRows_);
3198    while (problemStatus_<0) {
3199      int iRow,iColumn;
3200      // clear
3201      for (iRow=0;iRow<4;iRow++) {
3202        rowArray_[iRow]->clear();
3203      }   
3204     
3205      for (iColumn=0;iColumn<2;iColumn++) {
3206        columnArray_[iColumn]->clear();
3207      }   
3208     
3209      // give matrix (and model costs and bounds a chance to be
3210      // refreshed (normally null)
3211      matrix_->refresh(this);
3212      // If getting nowhere - why not give it a kick
3213#if 1
3214      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)&&(specialOptions_&4)==0
3215          &&initialStatus!=10) {
3216        perturb(1);
3217        matrix_->rhsOffset(this,true,false);
3218      }
3219#endif
3220      // If we have done no iterations - special
3221      if (lastGoodIteration_==numberIterations_&&factorType)
3222        factorType=3;
3223      if (saveModel) {
3224        // Doing sprint
3225        if (sequenceIn_<0||numberIterations_>=stopSprint) {
3226          problemStatus_=-1;
3227          originalModel(saveModel);
3228          saveModel=NULL;
3229          if (sequenceIn_<0&&numberIterations_<reasonableSprintIteration&&
3230              sprintPass>100)
3231            primalColumnPivot_->switchOffSprint();
3232          //lastSprintIteration=numberIterations_;
3233          printf("End small model\n");
3234        }
3235      }
3236         
3237      // may factorize, checks if problem finished
3238      statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
3239      if (initialStatus==10) {
3240        // cleanup phase
3241        if(initialIterations != numberIterations_) {
3242          if (numberDualInfeasibilities_>10000&&numberDualInfeasibilities_>10*initialNegDjs) {
3243            // getting worse - try perturbing
3244            if (perturbation_<101&&(specialOptions_&4)==0) {
3245              perturb(1);
3246              matrix_->rhsOffset(this,true,false);
3247              statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
3248            }
3249          }
3250        } else {
3251          // save number of negative djs
3252          if (!numberPrimalInfeasibilities_)
3253            initialNegDjs=numberDualInfeasibilities_;
3254          // make sure weight won't be changed
3255          if (infeasibilityCost_==1.0e10)
3256            infeasibilityCost_=1.000001e10;
3257        }
3258      }
3259      // See if sprint says redo because of problems
3260      if (numberDualInfeasibilities_==-776) {
3261        // Need new set of variables
3262        problemStatus_=-1;
3263        originalModel(saveModel);
3264        saveModel=NULL;
3265        //lastSprintIteration=numberIterations_;
3266        printf("End small model after\n");
3267        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
3268      } 
3269      int numberSprintIterations=0;
3270      int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
3271      if (problemStatus_==777) {
3272        // problems so do one pass with normal
3273        problemStatus_=-1;
3274        originalModel(saveModel);
3275        saveModel=NULL;
3276        // Skip factorization
3277        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3278        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
3279      } else if (problemStatus_<0&&!saveModel&&numberSprintColumns&&firstFree_<0) {
3280        int numberSort=0;
3281        int numberFixed=0;
3282        int numberBasic=0;
3283        reasonableSprintIteration = numberIterations_ + 100;
3284        int * whichColumns = new int[numberColumns_];
3285        double * weight = new double[numberColumns_];
3286        int numberNegative=0;
3287        double sumNegative = 0.0;
3288        // now massage weight so all basic in plus good djs
3289        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3290          double dj = dj_[iColumn];
3291          switch(getColumnStatus(iColumn)) {
3292           
3293          case basic:
3294            dj = -1.0e50;
3295            numberBasic++;
3296            break;
3297          case atUpperBound:
3298            dj = -dj;
3299            break;
3300          case isFixed:
3301            dj=1.0e50;
3302            numberFixed++;
3303            break;
3304          case atLowerBound:
3305            dj = dj;
3306            break;
3307          case isFree:
3308            dj = -100.0*fabs(dj);
3309              break;
3310          case superBasic:
3311            dj = -100.0*fabs(dj);
3312            break;
3313          }
3314          if (dj<-dualTolerance_&&dj>-1.0e50) {
3315            numberNegative++;
3316            sumNegative -= dj;
3317          }
3318          weight[iColumn]=dj;
3319          whichColumns[iColumn] = iColumn;
3320        }
3321        handler_->message(CLP_SPRINT,messages_)
3322          <<sprintPass<<numberIterations_-lastSprintIteration<<objectiveValue()<<sumNegative
3323          <<numberNegative
3324          <<CoinMessageEol;
3325        sprintPass++;
3326        lastSprintIteration=numberIterations_;
3327        if (objectiveValue()*optimizationDirection_>lastObjectiveValue-1.0e-7&&sprintPass>5) {
3328          // switch off
3329          printf("Switching off sprint\n");
3330          primalColumnPivot_->switchOffSprint();
3331        } else {
3332          lastObjectiveValue = objectiveValue()*optimizationDirection_;
3333          // sort
3334          CoinSort_2(weight,weight+numberColumns_,whichColumns);
3335          numberSort = CoinMin(numberColumns_-numberFixed,numberBasic+numberSprintColumns);
3336          // Sort to make consistent ?
3337          std::sort(whichColumns,whichColumns+numberSort);
3338          saveModel = new ClpSimplex(this,numberSort,whichColumns);
3339          delete [] whichColumns;
3340          delete [] weight;
3341          // Skip factorization
3342          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3343          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
3344          stopSprint = numberIterations_+numberSprintIterations;
3345          printf("Sprint with %d columns for %d iterations\n",
3346                 numberSprintColumns,numberSprintIterations);
3347        }
3348      }
3349     
3350      // Say good factorization
3351      factorType=1;
3352     
3353      // Say no pivot has occurred (for steepest edge and updates)
3354      pivotRow_=-2;
3355
3356      // exit if victory declared
3357      if (problemStatus_>=0) {
3358        if (originalCost) {
3359          // find number nonbasic with zero reduced costs
3360          int numberDegen=0;
3361          int numberTotal = numberColumns_; //+numberRows_;
3362          for (int i=0;i<numberTotal;i++) {
3363            cost_[i]=0.0;
3364            if (getStatus(i)==atLowerBound) {
3365              if (dj_[i]<=dualTolerance_) {
3366                cost_[i]=numberTotal-i+randomNumberGenerator_.randomDouble()*0.5;
3367                numberDegen++;
3368              } else {
3369                // fix
3370                cost_[i]=1.0e10;//upper_[i]=lower_[i];
3371              }
3372            } else if (getStatus(i)==atUpperBound) {
3373              if (dj_[i]>=-dualTolerance_) {
3374                cost_[i]=(numberTotal-i)+randomNumberGenerator_.randomDouble()*0.5;
3375                numberDegen++;
3376              } else {
3377                // fix
3378                cost_[i]=-1.0e10;//lower_[i]=upper_[i];
3379              } 
3380            } else if (getStatus(i)==basic) {
3381              cost_[i] = (numberTotal-i)+randomNumberGenerator_.randomDouble()*0.5;
3382            }
3383          }
3384          problemStatus_=-1;
3385          lastObjectiveValue=COIN_DBL_MAX;
3386          // Start check for cycles
3387          progress_.fillFromModel(this);
3388          progress_.startCheck();
3389          printf("%d degenerate after %d iterations\n",numberDegen,
3390                 numberIterations_);
3391          if (!numberDegen) {
3392            CoinMemcpyN(originalCost,numberTotal,cost_);
3393            delete [] originalCost;
3394            originalCost=NULL;
3395            CoinMemcpyN(originalLower,numberTotal,lower_);
3396            delete [] originalLower;
3397            CoinMemcpyN(originalUpper,numberTotal,upper_);
3398            delete [] originalUpper;
3399          }
3400          delete nonLinearCost_;
3401          nonLinearCost_ = new ClpNonLinearCost(this);
3402          progress_.endOddState();
3403          continue;
3404        } else { 
3405          printf("exiting after %d iterations\n",numberIterations_);
3406          break;
3407        }
3408      }
3409     
3410      // test for maximum iterations
3411      if (hitMaximumIterations()||(ifValuesPass==2&&firstFree_<0)) {
3412        problemStatus_=3;
3413        break;
3414      }
3415
3416      if (firstFree_<0) {
3417        if (ifValuesPass) {
3418          // end of values pass
3419          ifValuesPass=0;
3420          int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
3421          if (status>=0) {
3422            problemStatus_=5;
3423            secondaryStatus_=ClpEventHandler::endOfValuesPass;
3424            break;
3425          }
3426        }
3427      }
3428      // Check event
3429      {
3430        int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3431        if (status>=0) {
3432          problemStatus_=5;
3433          secondaryStatus_=ClpEventHandler::endOfFactorization;
3434          break;
3435        }
3436      }
3437      // Iterate
3438      whileIterating(ifValuesPass ? 1 : 0);
3439    }
3440  }
3441  assert (!originalCost);
3442  // if infeasible get real values
3443  //printf("XXXXY final cost %g\n",infeasibilityCost_);
3444  progress_.initialWeight_=0.0;
3445  if (problemStatus_==1&&secondaryStatus_!=6) {
3446    infeasibilityCost_=0.0;
3447    createRim(1+4);
3448    nonLinearCost_->checkInfeasibilities(0.0);
3449    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
3450    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
3451    // and get good feasible duals
3452    computeDuals(NULL);
3453  }
3454  // clean up
3455  unflag();
3456  finish(0);
3457  restoreData(data);
3458  return problemStatus_;
3459}
Note: See TracBrowser for help on using the repository browser.