source: trunk/Clp/src/AbcSimplexPrimal.cpp @ 1879

Last change on this file since 1879 was 1879, checked in by forrest, 7 years ago

fix a few problems with Aboca

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