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

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

minor changes to implement Aboca

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