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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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