source: branches/devel/Clp/src/ClpSimplexPrimal.cpp @ 989

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

this may be a mistake

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