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

Last change on this file since 871 was 871, checked in by forrest, 15 years ago

possible newline error

  • 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        difference = CoinMin(difference,perturbation);
2055        difference = CoinMin(difference,fabs(solutionValue)+1.0);
2056        double value = maximumFraction*(difference+1.0);
2057        value = CoinMin(value,0.1);
2058        value *= CoinDrand48();
2059        if (solutionValue-lowerValue<=primalTolerance_) {
2060          lower_[iSequence] -= value;
2061        } else if (upperValue-solutionValue<=primalTolerance_) {
2062          upper_[iSequence] += value;
2063        } else {
2064#if 0
2065          if (iSequence>=numberColumns_) {
2066            // may not be at bound - but still perturb (unless free)
2067            if (upperValue>1.0e30&&lowerValue<-1.0e30)
2068              value=0.0;
2069            else
2070              value = - value; // as -1.0 in matrix
2071          } else {
2072            value = 0.0;
2073          }
2074#else
2075          value=0.0;
2076#endif
2077        }
2078        if (value) {
2079          if (printOut)
2080            printf("col %d lower from %g to %g, upper from %g to %g\n",
2081                   iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
2082          if (solutionValue) {
2083            largest = CoinMax(largest,value);
2084            if (value>(fabs(solutionValue)+1.0)*largestPerCent)
2085              largestPerCent=value/(fabs(solutionValue)+1.0);
2086          } else {
2087            largestZero = CoinMax(largestZero,value);
2088          } 
2089        }
2090      }
2091    }
2092  } else {
2093    double tolerance = 100.0*primalTolerance_;
2094    for (i=0;i<numberColumns_;i++) {
2095      double lowerValue=lower_[i], upperValue=upper_[i];
2096      if (upperValue>lowerValue+primalTolerance_) {
2097        double value = perturbation*maximumFraction;
2098        value = CoinMin(value,0.1);
2099        value *= CoinDrand48();
2100        if (savePerturbation!=50) {
2101          if (fabs(value)<=primalTolerance_)
2102            value=0.0;
2103          if (lowerValue>-1.0e20&&lowerValue)
2104            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2105          if (upperValue<1.0e20&&upperValue)
2106            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2107        } else if (value) {
2108          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2109          // get in range
2110          if (valueL<=tolerance) {
2111            valueL *= 10.0;
2112            while (valueL<=tolerance) 
2113              valueL *= 10.0;
2114          } else if (valueL>1.0) {
2115            valueL *= 0.1;
2116            while (valueL>1.0) 
2117              valueL *= 0.1;
2118          }
2119          if (lowerValue>-1.0e20&&lowerValue)
2120            lowerValue -= valueL;
2121          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2122          // get in range
2123          if (valueU<=tolerance) {
2124            valueU *= 10.0;
2125            while (valueU<=tolerance) 
2126              valueU *= 10.0;
2127          } else if (valueU>1.0) {
2128            valueU *= 0.1;
2129            while (valueU>1.0) 
2130              valueU *= 0.1;
2131          }
2132          if (upperValue<1.0e20&&upperValue)
2133            upperValue += valueU;
2134        }
2135        if (lowerValue!=lower_[i]) {
2136          double difference = fabs(lowerValue-lower_[i]);
2137          largest = CoinMax(largest,difference);
2138          if (difference>fabs(lower_[i])*largestPerCent)
2139            largestPerCent=fabs(difference/lower_[i]);
2140        } 
2141        if (upperValue!=upper_[i]) {
2142          double difference = fabs(upperValue-upper_[i]);
2143          largest = CoinMax(largest,difference);
2144          if (difference>fabs(upper_[i])*largestPerCent)
2145            largestPerCent=fabs(difference/upper_[i]);
2146        } 
2147        if (printOut)
2148          printf("col %d lower from %g to %g, upper from %g to %g\n",
2149                 i,lower_[i],lowerValue,upper_[i],upperValue);
2150      }
2151      lower_[i]=lowerValue;
2152      upper_[i]=upperValue;
2153    }
2154    for (;i<numberColumns_+numberRows_;i++) {
2155      double lowerValue=lower_[i], upperValue=upper_[i];
2156      double value = perturbation*maximumFraction;
2157      value = CoinMin(value,0.1);
2158      value *= CoinDrand48();
2159      if (upperValue>lowerValue+tolerance) {
2160        if (savePerturbation!=50) {
2161          if (fabs(value)<=primalTolerance_)
2162            value=0.0;
2163          if (lowerValue>-1.0e20&&lowerValue)
2164            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2165          if (upperValue<1.0e20&&upperValue)
2166            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2167        } else if (value) {
2168          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2169          // get in range
2170          if (valueL<=tolerance) {
2171            valueL *= 10.0;
2172            while (valueL<=tolerance) 
2173              valueL *= 10.0;
2174          } else if (valueL>1.0) {
2175            valueL *= 0.1;
2176            while (valueL>1.0) 
2177              valueL *= 0.1;
2178          }
2179          if (lowerValue>-1.0e20&&lowerValue)
2180            lowerValue -= valueL;
2181          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2182          // get in range
2183          if (valueU<=tolerance) {
2184            valueU *= 10.0;
2185            while (valueU<=tolerance) 
2186              valueU *= 10.0;
2187          } else if (valueU>1.0) {
2188            valueU *= 0.1;
2189            while (valueU>1.0) 
2190              valueU *= 0.1;
2191          }
2192          if (upperValue<1.0e20&&upperValue)
2193            upperValue += valueU;
2194        }
2195      } else if (upperValue>0.0) {
2196        upperValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2197        lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2198      } else if (upperValue<0.0) {
2199        upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2200        lowerValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2201      } else {
2202      }
2203      if (lowerValue!=lower_[i]) {
2204        double difference = fabs(lowerValue-lower_[i]);
2205        largest = CoinMax(largest,difference);
2206        if (difference>fabs(lower_[i])*largestPerCent)
2207          largestPerCent=fabs(difference/lower_[i]);
2208      } 
2209      if (upperValue!=upper_[i]) {
2210        double difference = fabs(upperValue-upper_[i]);
2211        largest = CoinMax(largest,difference);
2212        if (difference>fabs(upper_[i])*largestPerCent)
2213          largestPerCent=fabs(difference/upper_[i]);
2214      } 
2215      if (printOut)
2216        printf("row %d lower from %g to %g, upper from %g to %g\n",
2217               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
2218      lower_[i]=lowerValue;
2219      upper_[i]=upperValue;
2220    }
2221  }
2222  // Clean up
2223  for (i=0;i<numberColumns_+numberRows_;i++) {
2224    switch(getStatus(i)) {
2225     
2226    case basic:
2227      break;
2228    case atUpperBound:
2229      solution_[i]=upper_[i];
2230      break;
2231    case isFixed:
2232    case atLowerBound:
2233      solution_[i]=lower_[i];
2234      break;
2235    case isFree:
2236      break;
2237    case superBasic:
2238      break;
2239    }
2240  }
2241  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2242    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
2243    <<CoinMessageEol;
2244  // redo nonlinear costs
2245  // say perturbed
2246  perturbation_=101;
2247}
2248// un perturb
2249bool
2250ClpSimplexPrimal::unPerturb()
2251{
2252  if (perturbation_!=101)
2253    return false;
2254  // put back original bounds and costs
2255  createRim(1+4);
2256  sanityCheck();
2257  // unflag
2258  unflag();
2259  // get a valid nonlinear cost function
2260  delete nonLinearCost_;
2261  nonLinearCost_= new ClpNonLinearCost(this);
2262  perturbation_ = 102; // stop any further perturbation
2263  // move non basic variables to new bounds
2264  nonLinearCost_->checkInfeasibilities(0.0);
2265#if 1
2266  // Try using dual
2267  return true;
2268#else
2269  gutsOfSolution(NULL,NULL,ifValuesPass!=0);
2270  return false;
2271#endif
2272 
2273}
2274// Unflag all variables and return number unflagged
2275int 
2276ClpSimplexPrimal::unflag()
2277{
2278  int i;
2279  int number = numberRows_+numberColumns_;
2280  int numberFlagged=0;
2281  // we can't really trust infeasibilities if there is dual error
2282  // allow tolerance bigger than standard to check on duals
2283  double relaxedToleranceD=dualTolerance_ + CoinMin(1.0e-2,10.0*largestDualError_);
2284  for (i=0;i<number;i++) {
2285    if (flagged(i)) {
2286      clearFlagged(i);
2287      // only say if reasonable dj
2288      if (fabs(dj_[i])>relaxedToleranceD)
2289        numberFlagged++;
2290    }
2291  }
2292  numberFlagged += matrix_->generalExpanded(this,8,i);
2293  if (handler_->logLevel()>2&&numberFlagged&&objective_->type()>1)
2294    printf("%d unflagged\n",numberFlagged);
2295  return numberFlagged;
2296}
2297// Do not change infeasibility cost and always say optimal
2298void 
2299ClpSimplexPrimal::alwaysOptimal(bool onOff)
2300{
2301  if (onOff)
2302    specialOptions_ |= 1;
2303  else
2304    specialOptions_ &= ~1;
2305}
2306bool 
2307ClpSimplexPrimal::alwaysOptimal() const
2308{
2309  return (specialOptions_&1)!=0;
2310}
2311// Flatten outgoing variables i.e. - always to exact bound
2312void 
2313ClpSimplexPrimal::exactOutgoing(bool onOff)
2314{
2315  if (onOff)
2316    specialOptions_ |= 4;
2317  else
2318    specialOptions_ &= ~4;
2319}
2320bool 
2321ClpSimplexPrimal::exactOutgoing() const
2322{
2323  return (specialOptions_&4)!=0;
2324}
2325/*
2326  Reasons to come out (normal mode/user mode):
2327  -1 normal
2328  -2 factorize now - good iteration/ NA
2329  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
2330  -4 inaccuracy - refactorize - no iteration/ NA
2331  -5 something flagged - go round again/ pivot not possible
2332  +2 looks unbounded
2333  +3 max iterations (iteration done)
2334*/
2335int
2336ClpSimplexPrimal::pivotResult(int ifValuesPass)
2337{
2338
2339  bool roundAgain=true;
2340  int returnCode=-1;
2341
2342  // loop round if user setting and doing refactorization
2343  while (roundAgain) {
2344    roundAgain=false;
2345    returnCode=-1;
2346    pivotRow_=-1;
2347    sequenceOut_=-1;
2348    rowArray_[1]->clear();
2349#if 0
2350    {
2351      int seq[]={612,643};
2352      int k;
2353      for (k=0;k<sizeof(seq)/sizeof(int);k++) {
2354        int iSeq=seq[k];
2355        if (getColumnStatus(iSeq)!=basic) {
2356          double djval;
2357          double * work;
2358          int number;
2359          int * which;
2360         
2361          int iIndex;
2362          unpack(rowArray_[1],iSeq);
2363          factorization_->updateColumn(rowArray_[2],rowArray_[1]);
2364          djval = cost_[iSeq];
2365          work=rowArray_[1]->denseVector();
2366          number=rowArray_[1]->getNumElements();
2367          which=rowArray_[1]->getIndices();
2368         
2369          for (iIndex=0;iIndex<number;iIndex++) {
2370           
2371            int iRow = which[iIndex];
2372            double alpha = work[iRow];
2373            int iPivot=pivotVariable_[iRow];
2374            djval -= alpha*cost_[iPivot];
2375          }
2376          double comp = 1.0e-8 + 1.0e-7*(CoinMax(fabs(dj_[iSeq]),fabs(djval)));
2377          if (fabs(djval-dj_[iSeq])>comp)
2378            printf("Bad dj %g for %d - true is %g\n",
2379                   dj_[iSeq],iSeq,djval);
2380          assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
2381          rowArray_[1]->clear();
2382        }
2383      }
2384    }
2385#endif
2386       
2387    // we found a pivot column
2388    // update the incoming column
2389    unpackPacked(rowArray_[1]);
2390    // save reduced cost
2391    double saveDj = dualIn_;
2392    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2393    // Get extra rows
2394    matrix_->extendUpdated(this,rowArray_[1],0);
2395    // do ratio test and re-compute dj
2396    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2397              ifValuesPass);
2398    if (ifValuesPass) {
2399      saveDj=dualIn_;
2400      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
2401      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2402        if(fabs(dualIn_)<1.0e2*dualTolerance_&&objective_->type()<2) {
2403          // try other way
2404          directionIn_=-directionIn_;
2405          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2406                    0);
2407        }
2408        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2409          if (solveType_==1) {
2410            // reject it
2411            char x = isColumn(sequenceIn_) ? 'C' :'R';
2412            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2413              <<x<<sequenceWithin(sequenceIn_)
2414              <<CoinMessageEol;
2415            setFlagged(sequenceIn_);
2416            progress_->clearBadTimes();
2417            lastBadIteration_ = numberIterations_; // say be more cautious
2418            clearAll();
2419            pivotRow_=-1;
2420          }
2421          returnCode=-5;
2422          break;
2423        }
2424      }
2425    }
2426    // need to clear toIndex_ in gub
2427    // ? when can I clear stuff
2428    // Clean up any gub stuff
2429    matrix_->extendUpdated(this,rowArray_[1],1);
2430    double checkValue=1.0e-2;
2431    if (largestDualError_>1.0e-5)
2432      checkValue=1.0e-1;
2433    if (!ifValuesPass&&solveType_==1&&(saveDj*dualIn_<1.0e-20||
2434        fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj))||
2435                        fabs(dualIn_)<dualTolerance_)) {
2436      char x = isColumn(sequenceIn_) ? 'C' :'R';
2437      handler_->message(CLP_PRIMAL_DJ,messages_)
2438        <<x<<sequenceWithin(sequenceIn_)
2439        <<saveDj<<dualIn_
2440        <<CoinMessageEol;
2441      if(lastGoodIteration_ != numberIterations_) {
2442        clearAll();
2443        pivotRow_=-1; // say no weights update
2444        returnCode=-4;
2445        if(lastGoodIteration_+1 == numberIterations_) {
2446          // not looking wonderful - try cleaning bounds
2447          // put non-basics to bounds in case tolerance moved
2448          nonLinearCost_->checkInfeasibilities(0.0);
2449        }
2450        sequenceOut_=-1;
2451        break;
2452      } else {
2453        // take on more relaxed criterion
2454        if (saveDj*dualIn_<1.0e-20||
2455            fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))||
2456            fabs(dualIn_)<dualTolerance_) {
2457          // need to reject something
2458          char x = isColumn(sequenceIn_) ? 'C' :'R';
2459          handler_->message(CLP_SIMPLEX_FLAG,messages_)
2460            <<x<<sequenceWithin(sequenceIn_)
2461            <<CoinMessageEol;
2462          setFlagged(sequenceIn_);
2463          progress_->clearBadTimes();
2464          lastBadIteration_ = numberIterations_; // say be more cautious
2465          clearAll();
2466          pivotRow_=-1;
2467          returnCode=-5;
2468          sequenceOut_=-1;
2469          break;
2470        }
2471      }
2472    }
2473    if (pivotRow_>=0) {
2474      if (solveType_==2) {
2475        // **** Coding for user interface
2476        // do ray
2477        primalRay(rowArray_[1]);
2478        // update duals
2479        // as packed need to find pivot row
2480        //assert (rowArray_[1]->packedMode());
2481        //int i;
2482       
2483        //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2484        CoinAssert (fabs(alpha_)>1.0e-8);
2485        double multiplier = dualIn_/alpha_;
2486        rowArray_[0]->insert(pivotRow_,multiplier);
2487        factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
2488        // put row of tableau in rowArray[0] and columnArray[0]
2489        matrix_->transposeTimes(this,-1.0,
2490                                rowArray_[0],columnArray_[1],columnArray_[0]);
2491        // update column djs
2492        int i;
2493        int * index = columnArray_[0]->getIndices();
2494        int number = columnArray_[0]->getNumElements();
2495        double * element = columnArray_[0]->denseVector();
2496        for (i=0;i<number;i++) {
2497          int ii = index[i];
2498          dj_[ii] += element[ii];
2499          reducedCost_[ii] = dj_[ii];
2500          element[ii]=0.0;
2501        }
2502        columnArray_[0]->setNumElements(0);
2503        // and row djs
2504        index = rowArray_[0]->getIndices();
2505        number = rowArray_[0]->getNumElements();
2506        element = rowArray_[0]->denseVector();
2507        for (i=0;i<number;i++) {
2508          int ii = index[i];
2509          dj_[ii+numberColumns_] += element[ii];
2510          dual_[ii] = dj_[ii+numberColumns_];
2511          element[ii]=0.0;
2512        }
2513        rowArray_[0]->setNumElements(0);
2514        // check incoming
2515        CoinAssert (fabs(dj_[sequenceIn_])<1.0e-1);
2516      }
2517      // if stable replace in basis
2518      // If gub or odd then alpha and pivotRow may change
2519      int updateType=0;
2520      int updateStatus = matrix_->generalExpanded(this,3,updateType);
2521      if (updateType>=0)
2522        updateStatus = factorization_->replaceColumn(this,
2523                                                     rowArray_[2],
2524                                                     rowArray_[1],
2525                                                     pivotRow_,
2526                                                     alpha_);
2527
2528      // if no pivots, bad update but reasonable alpha - take and invert
2529      if (updateStatus==2&&
2530          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2531        updateStatus=4;
2532      if (updateStatus==1||updateStatus==4) {
2533        // slight error
2534        if (factorization_->pivots()>5||updateStatus==4) {
2535          returnCode=-3;
2536        }
2537      } else if (updateStatus==2) {
2538        // major error
2539        // better to have small tolerance even if slower
2540        factorization_->zeroTolerance(1.0e-15);
2541        int maxFactor = factorization_->maximumPivots();
2542        if (maxFactor>10) {
2543          if (forceFactorization_<0)
2544            forceFactorization_= maxFactor;
2545          forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2546        } 
2547        // later we may need to unwind more e.g. fake bounds
2548        if(lastGoodIteration_ != numberIterations_) {
2549          clearAll();
2550          pivotRow_=-1;
2551          if (solveType_==1) {
2552            returnCode=-4;
2553            break;
2554          } else {
2555            // user in charge - re-factorize
2556            int lastCleaned;
2557            ClpSimplexProgress dummyProgress;
2558            if (saveStatus_)
2559              statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2560            else
2561              statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2562            roundAgain=true;
2563            continue;
2564          }
2565        } else {
2566          // need to reject something
2567          if (solveType_==1) {
2568            char x = isColumn(sequenceIn_) ? 'C' :'R';
2569            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2570              <<x<<sequenceWithin(sequenceIn_)
2571              <<CoinMessageEol;
2572            setFlagged(sequenceIn_);
2573            progress_->clearBadTimes();
2574          }
2575          lastBadIteration_ = numberIterations_; // say be more cautious
2576          clearAll();
2577          pivotRow_=-1;
2578          sequenceOut_=-1;
2579          returnCode = -5;
2580          break;
2581
2582        }
2583      } else if (updateStatus==3) {
2584        // out of memory
2585        // increase space if not many iterations
2586        if (factorization_->pivots()<
2587            0.5*factorization_->maximumPivots()&&
2588            factorization_->pivots()<200)
2589          factorization_->areaFactor(
2590                                     factorization_->areaFactor() * 1.1);
2591        returnCode =-2; // factorize now
2592      } else if (updateStatus==5) {
2593        problemStatus_=-2; // factorize now
2594      }
2595      // here do part of steepest - ready for next iteration
2596      if (!ifValuesPass)
2597        primalColumnPivot_->updateWeights(rowArray_[1]);
2598    } else {
2599      if (pivotRow_==-1) {
2600        // no outgoing row is valid
2601        if (valueOut_!=COIN_DBL_MAX) {
2602          double objectiveChange=0.0;
2603          theta_=valueOut_-valueIn_;
2604          updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2605          solution_[sequenceIn_] += theta_;
2606        }
2607        rowArray_[0]->clear();
2608        if (!factorization_->pivots()&&acceptablePivot_<=1.0e-8) {
2609          returnCode = 2; //say looks unbounded
2610          // do ray
2611          primalRay(rowArray_[1]);
2612        } else if (solveType_==2) {
2613          // refactorize
2614          int lastCleaned;
2615          ClpSimplexProgress dummyProgress;
2616          if (saveStatus_)
2617            statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2618          else
2619            statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2620          roundAgain=true;
2621          continue;
2622        } else {
2623          acceptablePivot_=1.0e-8;
2624          returnCode = 4; //say looks unbounded but has iterated
2625        }
2626        break;
2627      } else {
2628        // flipping from bound to bound
2629      }
2630    }
2631
2632    double oldCost = 0.0;
2633    if (sequenceOut_>=0)
2634      oldCost=cost_[sequenceOut_];
2635    // update primal solution
2636   
2637    double objectiveChange=0.0;
2638    // after this rowArray_[1] is not empty - used to update djs
2639    // If pivot row >= numberRows then may be gub
2640    int savePivot = pivotRow_;
2641    if (pivotRow_>=numberRows_)
2642      pivotRow_=-1;
2643    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2644    pivotRow_=savePivot;
2645   
2646    double oldValue = valueIn_;
2647    if (directionIn_==-1) {
2648      // as if from upper bound
2649      if (sequenceIn_!=sequenceOut_) {
2650        // variable becoming basic
2651        valueIn_ -= fabs(theta_);
2652      } else {
2653        valueIn_=lowerIn_;
2654      }
2655    } else {
2656      // as if from lower bound
2657      if (sequenceIn_!=sequenceOut_) {
2658        // variable becoming basic
2659        valueIn_ += fabs(theta_);
2660      } else {
2661        valueIn_=upperIn_;
2662      }
2663    }
2664    objectiveChange += dualIn_*(valueIn_-oldValue);
2665    // outgoing
2666    if (sequenceIn_!=sequenceOut_) {
2667      if (directionOut_>0) {
2668        valueOut_ = lowerOut_;
2669      } else {
2670        valueOut_ = upperOut_;
2671      }
2672      if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2673        valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2674      else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2675        valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2676      // may not be exactly at bound and bounds may have changed
2677      // Make sure outgoing looks feasible
2678      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2679      // May have got inaccurate
2680      //if (oldCost!=cost_[sequenceOut_])
2681      //printf("costchange on %d from %g to %g\n",sequenceOut_,
2682      //       oldCost,cost_[sequenceOut_]);
2683      if (solveType_!=2)
2684        dj_[sequenceOut_]=cost_[sequenceOut_]-oldCost; // normally updated next iteration
2685      solution_[sequenceOut_]=valueOut_;
2686    }
2687    // change cost and bounds on incoming if primal
2688    nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2689    int whatNext=housekeeping(objectiveChange);
2690    //nonLinearCost_->validate();
2691#if CLP_DEBUG >1
2692    {
2693      double sum;
2694      int ninf= matrix_->checkFeasible(this,sum);
2695      if (ninf)
2696        printf("infeas %d\n",ninf);
2697    }
2698#endif
2699    if (whatNext==1) {
2700        returnCode =-2; // refactorize
2701    } else if (whatNext==2) {
2702      // maximum iterations or equivalent
2703      returnCode=3;
2704    } else if(numberIterations_ == lastGoodIteration_
2705              + 2 * factorization_->maximumPivots()) {
2706      // done a lot of flips - be safe
2707      returnCode =-2; // refactorize
2708    }
2709    // Check event
2710    {
2711      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2712      if (status>=0) {
2713        problemStatus_=5;
2714        secondaryStatus_=ClpEventHandler::endOfIteration;
2715        returnCode=4;
2716      }
2717    }
2718  }
2719  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
2720    // refactorize here
2721    int lastCleaned;
2722    ClpSimplexProgress dummyProgress;
2723    if (saveStatus_)
2724      statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2725    else
2726      statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2727    if (problemStatus_==5) {
2728      printf("Singular basis\n");
2729      problemStatus_=-1;
2730      returnCode=5;
2731    }
2732  }
2733#ifdef CLP_DEBUG
2734  {
2735    int i;
2736    // not [1] as may have information
2737    for (i=0;i<4;i++) {
2738      if (i!=1)
2739        rowArray_[i]->checkClear();
2740    }   
2741    for (i=0;i<2;i++) {
2742      columnArray_[i]->checkClear();
2743    }   
2744  }     
2745#endif
2746  return returnCode;
2747}
2748// Create primal ray
2749void 
2750ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
2751{
2752  delete [] ray_;
2753  ray_ = new double [numberColumns_];
2754  CoinZeroN(ray_,numberColumns_);
2755  int number=rowArray->getNumElements();
2756  int * index = rowArray->getIndices();
2757  double * array = rowArray->denseVector();
2758  double way=-directionIn_;
2759  int i;
2760  double zeroTolerance=1.0e-12;
2761  if (sequenceIn_<numberColumns_)
2762    ray_[sequenceIn_]=directionIn_;
2763  if (!rowArray->packedMode()) {
2764    for (i=0;i<number;i++) {
2765      int iRow=index[i];
2766      int iPivot=pivotVariable_[iRow];
2767      double arrayValue = array[iRow];
2768      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2769        ray_[iPivot] = way* arrayValue;
2770    }
2771  } else {
2772    for (i=0;i<number;i++) {
2773      int iRow=index[i];
2774      int iPivot=pivotVariable_[iRow];
2775      double arrayValue = array[i];
2776      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2777        ray_[iPivot] = way* arrayValue;
2778    }
2779  }
2780}
2781/* Get next superbasic -1 if none,
2782   Normal type is 1
2783   If type is 3 then initializes sorted list
2784   if 2 uses list.
2785*/
2786int 
2787ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
2788{
2789  if (firstFree_>=0&&superBasicType) {
2790    int returnValue=-1;
2791    bool finished=false;
2792    while (!finished) {
2793      returnValue=firstFree_;
2794      int iColumn=firstFree_+1;
2795      if (superBasicType>1) {
2796        if (superBasicType>2) {
2797          // Initialize list
2798          // Wild guess that lower bound more natural than upper
2799          int number=0;
2800          double * work=columnArray->denseVector();
2801          int * which=columnArray->getIndices();
2802          for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2803            if (!flagged(iColumn)) {
2804              if (getStatus(iColumn)==superBasic) {
2805                if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2806                  solution_[iColumn]=lower_[iColumn];
2807                  setStatus(iColumn,atLowerBound);
2808                } else if (fabs(solution_[iColumn]-upper_[iColumn])
2809                           <=primalTolerance_) {
2810                  solution_[iColumn]=upper_[iColumn];
2811                  setStatus(iColumn,atUpperBound);
2812                } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2813                  setStatus(iColumn,isFree);
2814                  break;
2815                } else if (!flagged(iColumn)) {
2816                  // put ones near bounds at end after sorting
2817                  work[number]= - CoinMin(0.1*(solution_[iColumn]-lower_[iColumn]),
2818                                      upper_[iColumn]-solution_[iColumn]);
2819                  which[number++] = iColumn;
2820                }
2821              }
2822            }
2823          }
2824          CoinSort_2(work,work+number,which);
2825          columnArray->setNumElements(number);
2826          CoinZeroN(work,number);
2827        }
2828        int * which=columnArray->getIndices();
2829        int number = columnArray->getNumElements();
2830        if (!number) {
2831          // finished
2832          iColumn = numberRows_+numberColumns_;
2833          returnValue=-1;
2834        } else {
2835          number--;
2836          returnValue=which[number];
2837          iColumn=returnValue;
2838          columnArray->setNumElements(number);
2839        }     
2840      } else {
2841        for (;iColumn<numberRows_+numberColumns_;iColumn++) {
2842          if (!flagged(iColumn)) {
2843            if (getStatus(iColumn)==superBasic) {
2844              if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2845                solution_[iColumn]=lower_[iColumn];
2846                setStatus(iColumn,atLowerBound);
2847              } else if (fabs(solution_[iColumn]-upper_[iColumn])
2848                         <=primalTolerance_) {
2849                solution_[iColumn]=upper_[iColumn];
2850                setStatus(iColumn,atUpperBound);
2851              } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2852                setStatus(iColumn,isFree);
2853                break;
2854              } else {
2855                break;
2856              }
2857            }
2858          }
2859        }
2860      }
2861      firstFree_ = iColumn;
2862      finished=true;
2863      if (firstFree_==numberRows_+numberColumns_)
2864        firstFree_=-1;
2865      if (returnValue>=0&&getStatus(returnValue)!=superBasic)
2866        finished=false; // somehow picked up odd one
2867    }
2868    return returnValue;
2869  } else {
2870    return -1;
2871  }
2872}
2873void
2874ClpSimplexPrimal::clearAll()
2875{
2876  // Clean up any gub stuff
2877  matrix_->extendUpdated(this,rowArray_[1],1);
2878  int number=rowArray_[1]->getNumElements();
2879  int * which=rowArray_[1]->getIndices();
2880 
2881  int iIndex;
2882  for (iIndex=0;iIndex<number;iIndex++) {
2883   
2884    int iRow = which[iIndex];
2885    clearActive(iRow);
2886  }
2887  rowArray_[1]->clear();
2888  // make sure any gub sets are clean
2889  matrix_->generalExpanded(this,11,sequenceIn_);
2890 
2891}
Note: See TracBrowser for help on using the repository browser.