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

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

trying to make faster

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