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

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

minor changes for stability

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