source: branches/devel/Clp/src/ClpSimplexPrimal.cpp @ 889

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

for sprint

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