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

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

for threaded random numbers

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