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

Last change on this file since 2074 was 2074, checked in by forrest, 5 years ago

changes to try and make more stable

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