source: trunk/ClpSimplexPrimal.cpp @ 437

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

quadratic

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