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

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