source: stable/1.4/Clp/src/ClpSimplexPrimal.cpp @ 1058

Last change on this file since 1058 was 1058, checked in by forrest, 14 years ago

don't perturb fixed variables!

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