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

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

changes for cholesky including code from Anshul Gupta

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