# source:trunk/Clp/src/ClpSimplexPrimal.cpp@1517

Last change on this file since 1517 was 1517, checked in by forrest, 11 years ago

fix odd behavior on crossover from idiot

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