source: stable/1.6/Clp/src/ClpSimplexPrimal.cpp @ 1609

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

fix very rare bug

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