# source:stable/1.15/Clp/src/AbcSimplexPrimal.cpp@2006

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

changes for parallel and idiot

• Property svn:keywords set to `Id`
File size: 138.8 KB
Line
1/* \$Id: AbcSimplexPrimal.cpp 2006 2013-12-12 15:40:41Z forrest \$ */
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
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) {
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
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
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        }
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);
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
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;
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) {
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) {
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();
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_++;
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();
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_++;
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        }
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        }
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
2634                               CoinPartitionedVector * spareRow2,
2635                               CoinPartitionedVector * spareColumn1)
2636{
2637  for (int i=0;i<4;i++)
2638    multipleSequenceIn_[i]=-1;
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)
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  //double multiplier = perturbation*maximumFraction;
3018  for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3019    if (getInternalStatus(iSequence) == basic) {
3020      double lowerValue = abcLower_[iSequence];
3021      double upperValue = abcUpper_[iSequence];
3022      if (upperValue > lowerValue + tolerance) {
3023        double solutionValue = abcSolution_[iSequence];
3024        double difference = upperValue - lowerValue;
3025        difference = CoinMin(difference, perturbation);
3026        difference = CoinMin(difference, fabs(solutionValue) + 1.0);
3027        double value = maximumFraction * (difference + bias);
3028        value = CoinMin(value, 0.1);
3029        value = CoinMax(value,primalTolerance_);
3030        double perturbationValue=overallMultiplier*randomNumberGenerator_.randomDouble();
3031        value *= perturbationValue;
3032        if (solutionValue - lowerValue <= primalTolerance_) {
3033          abcLower_[iSequence] -= value;
3034          // change correct saved value
3035          if (abcNonLinearCost_->getCurrentStatus(iSequence)==CLP_FEASIBLE)
3036            abcPerturbation_[iSequence]=abcLower_[iSequence];
3037          else
3038            abcPerturbation_[iSequence+numberTotal_]=abcLower_[iSequence];
3039        } else if (upperValue - solutionValue <= primalTolerance_) {
3040          abcUpper_[iSequence] += value;
3041          // change correct saved value
3042          if (abcNonLinearCost_->getCurrentStatus(iSequence)==CLP_FEASIBLE)
3043            abcPerturbation_[iSequence+numberTotal_]=abcUpper_[iSequence];
3044          else
3045            abcPerturbation_[iSequence]=abcUpper_[iSequence];
3046        } else {
3047#if 0
3048          if (iSequence >= numberColumns_) {
3049            // may not be at bound - but still perturb (unless free)
3050            if (upperValue > 1.0e30 && lowerValue < -1.0e30)
3051              value = 0.0;
3052            else
3053              value = - value; // as -1.0 in matrix
3054          } else {
3055            value = 0.0;
3056          }
3057#else
3058          value = 0.0;
3059#endif
3060        }
3061        if (value) {
3062          numberChanged++;
3063          if (printOut)
3064            printf("col %d lower from %g to %g, upper from %g to %g\n",
3065                   iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue);
3066          if (solutionValue) {
3067            largest = CoinMax(largest, value);
3068            if (value > (fabs(solutionValue) + 1.0)*largestPerCent)
3069              largestPerCent = value / (fabs(solutionValue) + 1.0);
3070          } else {
3071            largestZero = CoinMax(largestZero, value);
3072          }
3073        }
3074      }
3075    }
3076  }
3077  if (!numberChanged) {
3078    // do non basic columns?
3079    for (iSequence = 0; iSequence < maximumAbcNumberRows_ + numberColumns_; iSequence++) {
3080      if (getInternalStatus(iSequence) != basic) {
3081        double lowerValue = abcLower_[iSequence];
3082        double upperValue = abcUpper_[iSequence];
3083        if (upperValue > lowerValue + tolerance) {
3084          double solutionValue = abcSolution_[iSequence];
3085          double difference = upperValue - lowerValue;
3086          difference = CoinMin(difference, perturbation);
3087          difference = CoinMin(difference, fabs(solutionValue) + 1.0);
3088          double value = maximumFraction * (difference + bias);
3089          value = CoinMin(value, 0.1);
3090          value = CoinMax(value,primalTolerance_);
3091          double perturbationValue=overallMultiplier*randomNumberGenerator_.randomDouble();
3092          value *= perturbationValue;
3093          if (solutionValue - lowerValue <= primalTolerance_) {
3094            abcLower_[iSequence] -= value;
3095            // change correct saved value
3096            if (abcNonLinearCost_->getCurrentStatus(iSequence)==CLP_FEASIBLE)
3097              abcPerturbation_[iSequence]=abcLower_[iSequence];
3098            else
3099              abcPerturbation_[iSequence+numberTotal_]=abcLower_[iSequence];
3100          } else if (upperValue - solutionValue <= primalTolerance_) {
3101            abcUpper_[iSequence] += value;
3102            // change correct saved value
3103            if (abcNonLinearCost_->getCurrentStatus(iSequence)==CLP_FEASIBLE)
3104              abcPerturbation_[iSequence+numberTotal_]=abcUpper_[iSequence];
3105            else
3106              abcPerturbation_[iSequence]=abcUpper_[iSequence];
3107          } else {
3108            value = 0.0;
3109          }
3110          if (value) {
3111            if (printOut)
3112              printf("col %d lower from %g to %g, upper from %g to %g\n",
3113                     iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue);
3114            if (solutionValue) {
3115              largest = CoinMax(largest, value);
3116              if (value > (fabs(solutionValue) + 1.0)*largestPerCent)
3117                largestPerCent = value / (fabs(solutionValue) + 1.0);
3118            } else {
3119              largestZero = CoinMax(largestZero, value);
3120            }
3121          }
3122        }
3123      }
3124    }
3125  }
3126  // Clean up
3127  for (i = 0; i < numberColumns_ + numberRows_; i++) {
3128    switch(getInternalStatus(i)) {
3129
3130    case basic:
3131      break;
3132    case atUpperBound:
3133      abcSolution_[i] = abcUpper_[i];
3134      break;
3135    case isFixed:
3136    case atLowerBound:
3137      abcSolution_[i] = abcLower_[i];
3138      break;
3139    case isFree:
3140      break;
3141    case superBasic:
3142      break;
3143    }
3144  }
3145  handler_->message(CLP_SIMPLEX_PERTURB, messages_)
3146    << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
3147    << CoinMessageEol;
3148  // redo nonlinear costs
3149  //delete abcNonLinearCost_;abcNonLinearCost_=new AbcNonLinearCost(this);//abort();// something elseabcNonLinearCost_->refresh();
3150  moveToBasic();
3151  if (!numberChanged) {
3152    // we changed nonbasic
3153    gutsOfPrimalSolution(3);
3154    // Make sure feasible bounds
3155    if (abcNonLinearCost_) {
3156      //abcNonLinearCost_->refresh();
3157      abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
3158      //abcNonLinearCost_->feasibleBounds();
3159    }
3160    abcPrimalColumnPivot_->saveWeights(this, 3);
3161  }
3162  // say perturbed
3163  perturbation_ = 101;
3164}
3165// un perturb
3166bool
3167AbcSimplexPrimal::unPerturb()
3168{
3169  if (perturbation_ != 101)
3170    return false;
3171  // put back original bounds and costs
3172  copyFromSaved();
3173  // copy bounds to perturbation
3174  CoinAbcMemcpy(abcPerturbation_,abcLower_,numberTotal_);
3175  CoinAbcMemcpy(abcPerturbation_+numberTotal_,abcUpper_,numberTotal_);
3176  //sanityCheck();
3177  // unflag
3178  unflag();
3179  abcProgress_.clearTimesFlagged();
3180  // get a valid nonlinear cost function
3181  delete abcNonLinearCost_;
3182  abcNonLinearCost_ = new AbcNonLinearCost(this);
3183  perturbation_ = 102; // stop any further perturbation
3184  // move non basic variables to new bounds
3185  abcNonLinearCost_->checkInfeasibilities(0.0);
3186  gutsOfSolution(3);
3187  abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
3188  // Try using dual
3189  return abcNonLinearCost_->sumInfeasibilities()!=0.0;
3190
3191}
3192// Unflag all variables and return number unflagged
3193int
3194AbcSimplexPrimal::unflag()
3195{
3196  int i;
3197  int number = numberRows_ + numberColumns_;
3198  int numberFlagged = 0;
3199  // we can't really trust infeasibilities if there is dual error
3200  // allow tolerance bigger than standard to check on duals
3201  double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
3202  for (i = 0; i < number; i++) {
3203    if (flagged(i)) {
3204      clearFlagged(i);
3205      // only say if reasonable dj
3206      if (fabs(abcDj_[i]) > relaxedToleranceD)
3207        numberFlagged++;
3208    }
3209  }
3210#if ABC_NORMAL_DEBUG>0
3211  if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
3212    printf("%d unflagged\n", numberFlagged);
3213#endif
3214  return numberFlagged;
3215}
3216// Do not change infeasibility cost and always say optimal
3217void
3218AbcSimplexPrimal::alwaysOptimal(bool onOff)
3219{
3220  if (onOff)
3221    specialOptions_ |= 1;
3222  else
3223    specialOptions_ &= ~1;
3224}
3225bool
3226AbcSimplexPrimal::alwaysOptimal() const
3227{
3228  return (specialOptions_ & 1) != 0;
3229}
3230// Flatten outgoing variables i.e. - always to exact bound
3231void
3232AbcSimplexPrimal::exactOutgoing(bool onOff)
3233{
3234  if (onOff)
3235    specialOptions_ |= 4;
3236  else
3237    specialOptions_ &= ~4;
3238}
3239bool
3240AbcSimplexPrimal::exactOutgoing() const
3241{
3242  return (specialOptions_ & 4) != 0;
3243}
3244/*
3245  Reasons to come out (normal mode/user mode):
3246  -1 normal
3247  -2 factorize now - good iteration/ NA
3248  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
3249  -4 inaccuracy - refactorize - no iteration/ NA
3250  -5 something flagged - go round again/ pivot not possible
3251  +2 looks unbounded
3252  +3 max iterations (iteration done)
3253*/
3254int
3255AbcSimplexPrimal::pivotResult(int ifValuesPass)
3256{
3257
3258  bool roundAgain = true;
3259  int returnCode = -1;
3260
3261  // loop round if user setting and doing refactorization
3262  while (roundAgain) {
3263    roundAgain = false;
3264    returnCode = -1;
3265    pivotRow_ = -1;
3266    sequenceOut_ = -1;
3267    usefulArray_[arrayForFtran_].clear();
3268
3269    // we found a pivot column
3270    // update the incoming column
3271    unpack(usefulArray_[arrayForFtran_]);
3272    // save reduced cost
3273    double saveDj = dualIn_;
3274    abcFactorization_->updateColumnFT(usefulArray_[arrayForFtran_]);
3275    // do ratio test and re-compute dj
3276#ifdef CLP_USER_DRIVEN
3277    if ((moreSpecialOptions_ & 512) == 0) {
3278#endif
3279      primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_],
3280                &usefulArray_[arrayForTableauRow_],
3281                ifValuesPass);
3282#ifdef CLP_USER_DRIVEN
3283      // user can tell which use it is
3284      int status = eventHandler_->event(ClpEventHandler::pivotRow);
3285      if (status >= 0) {
3286        problemStatus_ = 5;
3287        secondaryStatus_ = ClpEventHandler::pivotRow;
3288        break;
3289      }
3290    } else {
3291      int status = eventHandler_->event(ClpEventHandler::pivotRow);
3292      if (status >= 0) {
3293        problemStatus_ = 5;
3294        secondaryStatus_ = ClpEventHandler::pivotRow;
3295        break;
3296      }
3297    }
3298#endif
3299    if (ifValuesPass) {
3300      saveDj = dualIn_;
3301      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
3302      if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3303        if(fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
3304          // try other way
3305          directionIn_ = -directionIn_;
3306          primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_],
3307                    &usefulArray_[arrayForTableauRow_],
3308                    0);
3309        }
3310        if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3311          // reject it
3312          char x = isColumn(sequenceIn_) ? 'C' : 'R';
3313          handler_->message(CLP_SIMPLEX_FLAG, messages_)
3314            << x << sequenceWithin(sequenceIn_)
3315            << CoinMessageEol;
3316          setFlagged(sequenceIn_);
3317          abcProgress_.incrementTimesFlagged();
3319          lastBadIteration_ = numberIterations_; // say be more cautious
3320          clearAll();
3321          pivotRow_ = -1;
3322          returnCode = -5;
3323          break;
3324        }
3325      }
3326    }
3327    double checkValue = 1.0e-2;
3328    if (largestDualError_ > 1.0e-5)
3329      checkValue = 1.0e-1;
3330    double test2 = dualTolerance_;
3331    double test1 = 1.0e-20;
3332#if 0 //def FEB_TRY
3333    if (abcFactorization_->pivots() < 1) {
3334      test1 = -1.0e-4;
3335      if ((saveDj < 0.0 && dualIn_ < -1.0e-5 * dualTolerance_) ||
3336          (saveDj > 0.0 && dualIn_ > 1.0e-5 * dualTolerance_))
3337        test2 = 0.0; // allow through
3338    }
3339#endif
3340    if (!ifValuesPass && (saveDj * dualIn_ < test1 ||
3341                                             fabs(saveDj - dualIn_) > checkValue*(1.0 + fabs(saveDj)) ||
3342                          fabs(dualIn_) < test2)) {
3343      if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) >
3344            1.0e5)) {
3345        char x = isColumn(sequenceIn_) ? 'C' : 'R';
3346        handler_->message(CLP_PRIMAL_DJ, messages_)
3347          << x << sequenceWithin(sequenceIn_)
3348          << saveDj << dualIn_
3349          << CoinMessageEol;
3350        if(lastGoodIteration_ != numberIterations_) {
3351          clearAll();
3352          pivotRow_ = -1; // say no weights update
3353          returnCode = -4;
3354          if(lastGoodIteration_ + 1 == numberIterations_) {
3355            // not looking wonderful - try cleaning bounds
3356            // put non-basics to bounds in case tolerance moved
3357            abcNonLinearCost_->checkInfeasibilities(0.0);
3358          }
3359          sequenceOut_ = -1;
3360          sequenceIn_=-1;
3361          if (abcFactorization_->pivots()<10&&abcFactorization_->pivotTolerance()<0.25)
3362            abcFactorization_->saferTolerances(1.0,-1.03);
3363          break;
3364        } else {
3365          // take on more relaxed criterion
3366          if (saveDj * dualIn_ < test1 ||
3367                                 fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) ||
3368              fabs(dualIn_) < test2) {
3369            // need to reject something
3370            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3371            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3372              << x << sequenceWithin(sequenceIn_)
3373              << CoinMessageEol;
3374            setFlagged(sequenceIn_);
3375            abcProgress_.incrementTimesFlagged();
3376#if 1 //def FEB_TRY
3377            // Make safer?
3378            double tolerance=abcFactorization_->pivotTolerance();
3379            abcFactorization_->saferTolerances (1.0, -1.03);
3380#endif
3382            lastBadIteration_ = numberIterations_; // say be more cautious
3383            clearAll();
3384            pivotRow_ = -1;
3385            returnCode = -5;
3386            if(tolerance<abcFactorization_->pivotTolerance())
3387              returnCode=-4;
3388            sequenceOut_ = -1;
3389            sequenceIn_=-1;
3390            break;
3391          }
3392        }
3393      } else {
3394        //printf("%d %g %g\n",numberIterations_,saveDj,dualIn_);
3395      }
3396    }
3397    if (pivotRow_ >= 0) {
3398#ifdef CLP_USER_DRIVEN1
3399      // Got good pivot - may need to unflag stuff
3400      userChoiceWasGood(this);
3401#endif
3402      // if stable replace in basis
3403      // check update
3404      //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
3405      //usefulArray_[arrayForReplaceColumn_].print();
3406      ftAlpha_=abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
3408      abcFactorization_->replaceColumnPart3(this,
3409                                            &usefulArray_[arrayForReplaceColumn_],
3410                                            &usefulArray_[arrayForFtran_],
3411                                            pivotRow_,
3412                                            ftAlpha_?ftAlpha_:alpha_);
3413
3414      // if no pivots, bad update but reasonable alpha - take and invert
3415      if (updateStatus == 2 &&
3416          lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
3419        // slight error
3420        if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
3421          returnCode = -3;
3422        }
3423      } else if (updateStatus == 2) {
3424        // major error
3425        // better to have small tolerance even if slower
3426        abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3427        int maxFactor = abcFactorization_->maximumPivots();
3428        if (maxFactor > 10) {
3429          if (forceFactorization_ < 0)
3430            forceFactorization_ = maxFactor;
3431          forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3432        }
3433        // later we may need to unwind more e.g. fake bounds
3434        if(lastGoodIteration_ != numberIterations_) {
3435          clearAll();
3436          pivotRow_ = -1;
3437          sequenceIn_ = -1;
3438          sequenceOut_ = -1;
3439          returnCode = -4;
3440          break;
3441        } else {
3442          // need to reject something
3443          char x = isColumn(sequenceIn_) ? 'C' : 'R';
3444          handler_->message(CLP_SIMPLEX_FLAG, messages_)
3445            << x << sequenceWithin(sequenceIn_)
3446            << CoinMessageEol;
3447          setFlagged(sequenceIn_);
3448          abcProgress_.incrementTimesFlagged();
3450          lastBadIteration_ = numberIterations_; // say be more cautious
3451          clearAll();
3452          pivotRow_ = -1;
3453          sequenceIn_ = -1;
3454          sequenceOut_ = -1;
3455          returnCode = -5;
3456          break;
3457
3458        }
3459      } else if (updateStatus == 3) {
3460        // out of memory
3461        // increase space if not many iterations
3462        if (abcFactorization_->pivots() <
3463            0.5 * abcFactorization_->maximumPivots() &&
3464            abcFactorization_->pivots() < 200)
3465          abcFactorization_->areaFactor(
3466                                        abcFactorization_->areaFactor() * 1.1);
3467        returnCode = -2; // factorize now
3468      } else if (updateStatus == 5) {
3469        problemStatus_ = -2; // factorize now
3470      }
3471      // here do part of steepest - ready for next iteration
3472      if (!ifValuesPass)
3473        abcPrimalColumnPivot_->updateWeights(&usefulArray_[arrayForFtran_]);
3474    } else {
3475      if (pivotRow_ == -1) {
3476        // no outgoing row is valid
3477        if (valueOut_ != COIN_DBL_MAX) {
3478          double objectiveChange = 0.0;
3479          theta_ = valueOut_ - valueIn_;
3480          updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass);
3481          abcSolution_[sequenceIn_] += theta_;
3482        }
3483#ifdef CLP_USER_DRIVEN1
3484        /* Note if valueOut_ < COIN_DBL_MAX and
3485           theta_ reasonable then this may be a valid sub flip */
3486        if(!userChoiceValid2(this)) {
3487          if (abcFactorization_->pivots()<5) {
3488            // flag variable
3489            char x = isColumn(sequenceIn_) ? 'C' : 'R';
3490            handler_->message(CLP_SIMPLEX_FLAG, messages_)
3491              << x << sequenceWithin(sequenceIn_)
3492              << CoinMessageEol;
3493            setFlagged(sequenceIn_);
3494            abcProgress_.incrementTimesFlagged();
3496            roundAgain = true;
3497            continue;
3498          } else {
3499            // try refactorizing first
3500            returnCode = 4; //say looks odd but has iterated
3501            break;
3502          }
3503        }
3504#endif
3505        if (!abcFactorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3506          returnCode = 2; //say looks unbounded
3507          // do ray
3508          primalRay(&usefulArray_[arrayForFtran_]);
3509        } else {
3510          acceptablePivot_ = 1.0e-8;
3511          returnCode = 4; //say looks unbounded but has iterated
3512        }
3513        break;
3514      } else {
3515        // flipping from bound to bound
3516      }
3517    }
3518
3519    double oldCost = 0.0;
3520    if (sequenceOut_ >= 0)
3521      oldCost = abcCost_[sequenceOut_];
3522    // update primal solution
3523
3524    double objectiveChange = 0.0;
3525    // after this usefulArray_[arrayForFtran_] is not empty - used to update djs
3526#ifdef CLP_USER_DRIVEN
3527    if (theta_<0.0) {
3528      if (theta_>=-1.0e-12)
3529        theta_=0.0;
3530      //else
3531      //printf("negative theta %g\n",theta_);
3532    }
3533#endif
3534    updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass);
3535
3536    double oldValue = valueIn_;
3537    if (directionIn_ == -1) {
3538      // as if from upper bound
3539      if (sequenceIn_ != sequenceOut_) {
3540        // variable becoming basic
3541        valueIn_ -= fabs(theta_);
3542      } else {
3543        valueIn_ = lowerIn_;
3544      }
3545    } else {
3546      // as if from lower bound
3547      if (sequenceIn_ != sequenceOut_) {
3548        // variable becoming basic
3549        valueIn_ += fabs(theta_);
3550      } else {
3551        valueIn_ = upperIn_;
3552      }
3553    }
3554    objectiveChange += dualIn_ * (valueIn_ - oldValue);
3555    // outgoing
3556    if (sequenceIn_ != sequenceOut_) {
3557      if (directionOut_ > 0) {
3558        valueOut_ = lowerOut_;
3559      } else {
3560        valueOut_ = upperOut_;
3561      }
3562      if(valueOut_ < abcLower_[sequenceOut_] - primalTolerance_)
3563        valueOut_ = abcLower_[sequenceOut_] - 0.9 * primalTolerance_;
3564      else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_)
3565        valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_;
3566      // may not be exactly at bound and bounds may have changed
3567      // Make sure outgoing looks feasible
3568      directionOut_ = abcNonLinearCost_->setOneOutgoing(pivotRow_, valueOut_);
3569      // May have got inaccurate
3570      //if (oldCost!=abcCost_[sequenceOut_])
3571      //printf("costchange on %d from %g to %g\n",sequenceOut_,
3572      //       oldCost,abcCost_[sequenceOut_]);
3573      abcDj_[sequenceOut_] = abcCost_[sequenceOut_] - oldCost; // normally updated next iteration
3574      abcSolution_[sequenceOut_] = valueOut_;
3575    }
3576    // change cost and bounds on incoming if primal
3577    abcNonLinearCost_->setOne(sequenceIn_, valueIn_);
3578    int whatNext = housekeeping(/*objectiveChange*/);
3579    swapPrimalStuff();
3580    if (whatNext == 1) {
3581      returnCode = -2; // refactorize
3582    } else if (whatNext == 2) {
3583      // maximum iterations or equivalent
3584      returnCode = 3;
3585    } else if(numberIterations_ == lastGoodIteration_
3586              + 2 * abcFactorization_->maximumPivots()) {
3587      // done a lot of flips - be safe
3588      returnCode = -2; // refactorize
3589    }
3590    // Check event
3591    {
3592      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3593      if (status >= 0) {
3594        problemStatus_ = 5;
3595        secondaryStatus_ = ClpEventHandler::endOfIteration;
3596        returnCode = 3;
3597      }
3598    }
3599  }
3600  return returnCode;
3601}
3602/*
3603  Reasons to come out (normal mode/user mode):
3604  -1 normal
3605  -2 factorize now - good iteration/ NA
3606  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
3607  -4 inaccuracy - refactorize - no iteration/ NA
3608  -5 something flagged - go round again/ pivot not possible
3609  +2 looks unbounded
3610  +3 max iterations (iteration done)
3611*/
3612int
3613AbcSimplexPrimal::pivotResult4(int ifValuesPass)
3614{
3615
3616  int returnCode = -1;
3617  int numberMinor=0;
3618  int numberDone=0;
3619  CoinIndexedVector * vector[16];
3620  pivotStruct stuff[4];
3621  vector[0]=&usefulArray_[arrayForFtran_];
3622  vector[1]=&usefulArray_[arrayForFlipBounds_];
3623  vector[2]=&usefulArray_[arrayForTableauRow_];
3624  vector[3]=&usefulArray_[arrayForBtran_];
3625  /*
3626    For pricing we need to get a vector with difference in costs
3627    later we could modify code so only costBasic changed
3628    and store in first [4*x+1] array to go out the indices of changed
3629    plus at most four for pivots
3630  */
3631  //double * saveCosts=CoinCopyOfArray(costBasic_,numberRows_);
3632  //double saveCostsIn[4];
3633  for (int iMinor=0;iMinor<4;iMinor++) {
3634    int sequenceIn=multipleSequenceIn_[iMinor];
3635    if (sequenceIn<0)
3636      break;
3637    stuff[iMinor].valuesPass_=ifValuesPass;
3638    stuff[iMinor].lowerIn_=abcLower_[sequenceIn];
3639    stuff[iMinor].upperIn_=abcUpper_[sequenceIn];
3640    stuff[iMinor].valueIn_=abcSolution_[sequenceIn];
3641    stuff[iMinor].sequenceIn_=sequenceIn;
3642    //saveCostsIn[iMinor]=abcCost_[sequenceIn];
3643    numberMinor++;
3644    if (iMinor) {
3645      vector[4*iMinor]=rowArray_[2*iMinor-2];
3646      vector[4*iMinor+1]=rowArray_[2*iMinor-1];
3647      vector[4*iMinor+2]=columnArray_[2*iMinor-2];
3648      vector[4*iMinor+3]=columnArray_[2*iMinor-1];
3649    }
3650    for (int i=0;i<4;i++)
3651      vector[4*iMinor+i]->checkClear();
3652    unpack(*vector[4*iMinor],sequenceIn);
3653  }
3654  int numberLeft=numberMinor;
3655  // parallel (with cpu)
3656  for (int iMinor=1;iMinor<numberMinor;iMinor++) {
3657    // update the incoming columns
3658    cilk_spawn abcFactorization_->updateColumnFT(*vector[4*iMinor],*vector[4*iMinor+3],iMinor);
3659  }
3660  abcFactorization_->updateColumnFT(*vector[0],*vector[+3],0);
3661  cilk_sync;
3662  for (int iMinor=0;iMinor<numberMinor;iMinor++) {
3663    // find best (or first if values pass)
3664    int numberDo=1;
3665    int jMinor=0;
3666    while (jMinor<numberDo) {
3667      int sequenceIn=stuff[jMinor].sequenceIn_;
3668      double dj=abcDj_[sequenceIn];
3671        stuff[jMinor].dualIn_=dj;
3672        stuff[jMinor].saveDualIn_=dj;
3673        if (dj<0.0)
3674          stuff[jMinor].directionIn_=1;
3675        else
3676          stuff[jMinor].directionIn_=-1;
3677        jMinor++;
3678      } else {
3679        numberDo--;
3680        numberLeft--;
3681        // throw away
3682        for (int i=0;i<4;i++) {
3683          vector[4*jMinor+i]->clear();
3684          vector[4*jMinor+i]=vector[4*numberDo+i];
3685          vector[4*numberDo+i]=NULL;
3686        }
3687        stuff[jMinor]=stuff[numberDo];
3688      }
3689    }
3690    for (int jMinor=0;jMinor<numberDo;jMinor++) {
3691      // do ratio test and re-compute dj
3692      primalRow(vector[4*jMinor],vector[4*jMinor+1],vector[4*jMinor+2],stuff[jMinor]);
3693    }
3694    // choose best
3695    int iBest=-1;
3696    double bestMovement=-COIN_DBL_MAX;
3697    for (int jMinor=0;jMinor<numberDo;jMinor++) {
3698      double movement=stuff[jMinor].theta_*fabs(stuff[jMinor].dualIn_);
3699      if (movement>bestMovement) {
3700        bestMovement=movement;
3701        iBest=jMinor;
3702      }
3703    }
3704#if 0 //ndef MULTIPLE_PRICE
3705    if (maximumIterations()!=100000)
3706      iBest=0;
3707#endif
3708    if (iBest>=0) {
3709      dualIn_=stuff[iBest].dualIn_;
3710      dualOut_=stuff[iBest].dualOut_;
3711      lowerOut_=stuff[iBest].lowerOut_;
3712      upperOut_=stuff[iBest].upperOut_;
3713      valueOut_=stuff[iBest].valueOut_;
3714      sequenceOut_=stuff[iBest].sequenceOut_;
3715      sequenceIn_=stuff[iBest].sequenceIn_;
3716      lowerIn_=stuff[iBest].lowerIn_;
3717      upperIn_=stuff[iBest].upperIn_;
3718      valueIn_=stuff[iBest].valueIn_;
3719      directionIn_=stuff[iBest].directionIn_;
3720      directionOut_=stuff[iBest].directionOut_;
3721      pivotRow_=stuff[iBest].pivotRow_;
3722      theta_=stuff[iBest].theta_;
3723      alpha_=stuff[iBest].alpha_;
3724#ifdef MULTIPLE_PRICE
3725      for (int i=0;i<4*numberLeft;i++)
3726        vector[i]->clear();
3727      return 0;
3728#endif
3729      // maybe do this on more?
3730      double theta1 = CoinMax(theta_, 1.0e-12);
3731      double theta2 = numberIterations_ * abcNonLinearCost_->averageTheta();
3732      // Set average theta
3733      abcNonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast<double> (numberIterations_ + 1)));
3734      if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3735        if(fabs(dualIn_) < 1.0e2 * dualTolerance_) {
3736          // try other way
3737          stuff[iBest].directionIn_ = -directionIn_;
3738          stuff[iBest].valuesPass_=0;
3739          primalRow(vector[4*iBest],vector[4*iBest+1],vector[4*iBest+2],stuff[iBest]);
3740          dualIn_=stuff[iBest].dualIn_;
3741          dualOut_=stuff[iBest].dualOut_;
3742          lowerOut_=stuff[iBest].lowerOut_;
3743          upperOut_=stuff[iBest].upperOut_;
3744          valueOut_=stuff[iBest].valueOut_;
3745          sequenceOut_=stuff[iBest].sequenceOut_;
3746          sequenceIn_=stuff[iBest].sequenceIn_;
3747          lowerIn_=stuff[iBest].lowerIn_;
3748          upperIn_=stuff[iBest].upperIn_;
3749          valueIn_=stuff[iBest].valueIn_;
3750          directionIn_=stuff[iBest].directionIn_;
3751          directionOut_=stuff[iBest].directionOut_;
3752          pivotRow_=stuff[iBest].pivotRow_;
3753          theta_=stuff[iBest].theta_;
3754          alpha_=stuff[iBest].alpha_;
3755        }
3756        if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3757          // reject it
3758          char x = isColumn(sequenceIn_) ? 'C' : 'R';
3759          handler_->message(CLP_SIMPLEX_FLAG, messages_)
3760            << x << sequenceWithin(sequenceIn_)
3761            << CoinMessageEol;
3762          setFlagged(sequenceIn_);
3763          abcProgress_.incrementTimesFlagged();
3765          lastBadIteration_ = numberIterations_; // say be more cautious
3766          clearAll();
3767          pivotRow_ = -1;
3768          returnCode = -5;
3769          break;
3770        }
3771      }
3772      CoinIndexedVector * bestUpdate=NULL;
3773      if (pivotRow_ >= 0) {
3774        // Normal pivot
3775        // if stable replace in basis
3776        // check update
3777        CoinIndexedVector * tempVector[4];
3778        memcpy(tempVector,vector+4*iBest,4*sizeof(CoinIndexedVector *));
3779        returnCode = cilk_spawn doFTUpdate(tempVector);
3780        bestUpdate=tempVector[0];
3781        // after this bestUpdate is not empty - used to update djs ??
3782        cilk_spawn updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass);
3783        numberDone++;
3784        numberLeft--;
3785        // throw away
3786        for (int i=0;i<4;i++) {
3787          vector[4*iBest+i]=vector[4*numberLeft+i];
3788          vector[4*numberLeft+i]=NULL;
3789        }
3790        stuff[iBest]=stuff[numberLeft];
3791        // update pi and other vectors and FT stuff
3792        // parallel (can go 8 way?)
3793        for (int jMinor=0;jMinor<numberLeft;jMinor++) {
3794          cilk_spawn updatePartialUpdate(*vector[4*jMinor+3]);
3795        }
3796        // parallel
3797        for (int jMinor=0;jMinor<numberLeft;jMinor++) {
3798          stuff[jMinor].dualIn_=cilk_spawn updateMinorCandidate(*bestUpdate,*vector[4*jMinor],stuff[jMinor].sequenceIn_);
3799        }
3800        cilk_sync;
3801        // throw away
3802        for (int i=1;i<4;i++)
3803          tempVector[i]->clear();
3804        if (returnCode<-3) {
3805          clearAll();
3806          break;
3807        }
3808        // end Normal pivot
3809      } else {
3810        // start Flip
3811        vector[4*iBest+3]->clear();
3812        if (pivotRow_ == -1) {
3813          // no outgoing row is valid
3814          if (valueOut_ != COIN_DBL_MAX) {
3815            theta_ = valueOut_ - valueIn_;
3816            updatePrimalsInPrimal(*vector[4*iBest], theta_, ifValuesPass);
3817            abcSolution_[sequenceIn_] += theta_;
3818          }
3819          if (!abcFactorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3820            returnCode = 2; //say looks unbounded
3821            // do ray
3822            primalRay(vector[4*iBest]);
3823          } else {
3824            acceptablePivot_ = 1.0e-8;
3825            returnCode = 4; //say looks unbounded but has iterated
3826          }
3827          break;
3828        } else {
3829          // flipping from bound to bound
3830          bestUpdate=vector[4*iBest];
3831          // after this bestUpdate is not empty - used to update djs ??
3832          updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass);
3833          // throw away best BUT remember we are updating costs somehow
3834          numberDone++;
3835          numberLeft--;
3836          // throw away
3837          for (int i=0;i<4;i++) {
3838            if (i)
3839              vector[4*iBest+i]->clear();
3840            vector[4*iBest+i]=vector[4*numberLeft+i];
3841            vector[4*numberLeft+i]=NULL;
3842          }
3843          stuff[iBest]=stuff[numberLeft];
3844          // parallel
3845          for (int jMinor=0;jMinor<numberLeft;jMinor++) {
3846            stuff[jMinor].dualIn_=cilk_spawn updateMinorCandidate(*bestUpdate,*vector[4*jMinor],stuff[jMinor].sequenceIn_);
3847          }
3848          cilk_sync;
3849        }
3850        // end Flip
3851      }
3852      bestUpdate->clear(); // as only works in values pass
3853
3854      if (directionIn_ == -1) {
3855        // as if from upper bound
3856        if (sequenceIn_ != sequenceOut_) {
3857          // variable becoming basic
3858          valueIn_ -= fabs(theta_);
3859        } else {
3860          valueIn_ = lowerIn_;
3861        }
3862      } else {
3863        // as if from lower bound
3864        if (sequenceIn_ != sequenceOut_) {
3865          // variable becoming basic
3866          valueIn_ += fabs(theta_);
3867        } else {
3868          valueIn_ = upperIn_;
3869        }
3870      }
3871      // outgoing
3872      if (sequenceIn_ != sequenceOut_) {
3873        if (directionOut_ > 0) {
3874          valueOut_ = lowerOut_;
3875        } else {
3876          valueOut_ = upperOut_;
3877        }
3878        if(valueOut_ < abcLower_[sequenceOut_] - primalTolerance_)
3879          valueOut_ = abcLower_[sequenceOut_] - 0.9 * primalTolerance_;
3880        else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_)
3881          valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_;
3882        // may not be exactly at bound and bounds may have changed
3883        // Make sure outgoing looks feasible
3884        directionOut_ = abcNonLinearCost_->setOneOutgoing(pivotRow_, valueOut_);
3885        abcSolution_[sequenceOut_] = valueOut_;
3886      }
3887      // change cost and bounds on incoming if primal
3888      abcNonLinearCost_->setOne(sequenceIn_, valueIn_);
3889      int whatNext = housekeeping(/*objectiveChange*/);
3890      swapPrimalStuff();
3891      if (whatNext == 1) {
3892        returnCode = -2; // refactorize
3893      } else if (whatNext == 2) {
3894        // maximum iterations or equivalent
3895        returnCode = 3;
3896      } else if(numberIterations_ == lastGoodIteration_
3897                + 2 * abcFactorization_->maximumPivots()) {
3898        // done a lot of flips - be safe
3899        returnCode = -2; // refactorize
3900      }
3901      // Check event
3902      {
3903        int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3904        if (status >= 0) {
3905          problemStatus_ = 5;
3906          secondaryStatus_ = ClpEventHandler::endOfIteration;
3907          returnCode = 3;
3908        }
3909      }
3910    }
3911    if (returnCode!=-1)
3912      break;
3913  }
3914  for (int i=0;i<16;i++) {
3915    if (vector[i])
3916      vector[i]->clear();
3917    if (vector[i])
3918      vector[i]->checkClear();
3919  }
3920  //delete [] saveCosts;
3921  return returnCode;
3922}
3923// Do FT update as separate function for minor iterations (nonzero return code on problems)
3924int
3925AbcSimplexPrimal::doFTUpdate(CoinIndexedVector * vector[4])
3926{
3927  ftAlpha_=abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],
3928                                                vector[3],pivotRow_);
3931    abcFactorization_->replaceColumnPart3(this,
3932                                          &usefulArray_[arrayForReplaceColumn_],
3933                                          vector[0],
3934                                          vector[3],
3935                                          pivotRow_,
3936                                          ftAlpha_?ftAlpha_:alpha_);
3937  else
3938    vector[3]->clear();
3939  int returnCode=0;
3940  // if no pivots, bad update but reasonable alpha - take and invert
3941  if (updateStatus == 2 &&
3942      lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
3945    // slight error
3946    if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
3947      returnCode = -3;
3948    }
3949  } else if (updateStatus == 2) {
3950    // major error
3951    // better to have small tolerance even if slower
3952    abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3953    int maxFactor = abcFactorization_->maximumPivots();
3954    if (maxFactor > 10) {
3955      if (forceFactorization_ < 0)
3956        forceFactorization_ = maxFactor;
3957      forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3958    }
3959    // later we may need to unwind more e.g. fake bounds
3960    if(lastGoodIteration_ != numberIterations_) {
3961      pivotRow_ = -1;
3962      sequenceIn_ = -1;
3963      sequenceOut_ = -1;
3964      returnCode = -4;
3965    } else {
3966      // need to reject something
3967      char x = isColumn(sequenceIn_) ? 'C' : 'R';
3968      handler_->message(CLP_SIMPLEX_FLAG, messages_)
3969        << x << sequenceWithin(sequenceIn_)
3970        << CoinMessageEol;
3971      setFlagged(sequenceIn_);
3972      abcProgress_.incrementTimesFlagged();
3974      lastBadIteration_ = numberIterations_; // say be more cautious
3975      pivotRow_ = -1;
3976      sequenceIn_ = -1;
3977      sequenceOut_ = -1;
3978      returnCode = -5;
3979
3980    }
3981  } else if (updateStatus == 3) {
3982    // out of memory
3983    // increase space if not many iterations
3984    if (abcFactorization_->pivots() <
3985        0.5 * abcFactorization_->maximumPivots() &&
3986        abcFactorization_->pivots() < 200)
3987      abcFactorization_->areaFactor(
3988                                    abcFactorization_->areaFactor() * 1.1);
3989    returnCode = -2; // factorize now
3990  } else if (updateStatus == 5) {
3991    problemStatus_ = -2; // factorize now
3992    returnCode=-2;
3993  }
3994  return returnCode;
3995}
3996/* The primals are updated by the given array.
3997   costs are changed
3998*/
3999void
4000AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector & rowArray,
4001                                        double theta,bool valuesPass)
4002{
4003  // Cost on pivot row may change - may need to change dualIn
4004  //double oldCost = 0.0;
4005  //if (pivotRow_ >= 0)
4006  //oldCost = abcCost_[sequenceOut_];
4007  double * work = rowArray.denseVector();
4008  int number = rowArray.getNumElements();
4009  int * which = rowArray.getIndices();
4010  // allow for case where bound+tolerance == bound
4011  //double tolerance = 0.999999*primalTolerance_;
4012  double relaxedTolerance = 1.001 * primalTolerance_;
4013  int iIndex;
4014  if (!valuesPass) {
4015    for (iIndex = 0; iIndex < number; iIndex++) {
4016
4017      int iRow = which[iIndex];
4018      double alpha = work[iRow];
4019      //work[iRow] = 0.0;
4020      int iPivot = abcPivotVariable_[iRow];
4021      double change = theta * alpha;
4022      double value = abcSolution_[iPivot] - change;
4023      abcSolution_[iPivot] = value;
4024      value = solutionBasic_[iRow] - change;
4025      solutionBasic_[iRow] = value;
4026      if (active(iRow) || theta_ < 0.0) {
4027        clearActive(iRow);
4028        // But make sure one going out is feasible
4029        if (change > 0.0) {
4030          // going down
4031          if (value <= abcLower_[iPivot] + primalTolerance_) {
4032            if (iPivot == sequenceOut_ && value >= abcLower_[iPivot] - relaxedTolerance)
4033              value = abcLower_[iPivot];
4034            abcNonLinearCost_->setOneBasic(iRow, value);
4035          }
4036        } else {
4037          // going up
4038          if (value >= abcUpper_[iPivot] - primalTolerance_) {
4039            if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
4040              value = abcUpper_[iPivot];
4041            abcNonLinearCost_->setOneBasic(iRow, value);
4042          }
4043        }
4044      }
4045    }
4046  } else {
4047    // values pass so look at all
4048    for (iIndex = 0; iIndex < number; iIndex++) {
4049
4050      int iRow = which[iIndex];
4051      double alpha = work[iRow];
4052      //work[iRow] = 0.0;
4053      int iPivot = abcPivotVariable_[iRow];
4054      double change = theta * alpha;
4055      double value = abcSolution_[iPivot] - change;
4056      abcSolution_[iPivot] = value;
4057      value = solutionBasic_[iRow] - change;
4058      solutionBasic_[iRow] = value;
4059      clearActive(iRow);
4060      // But make sure one going out is feasible
4061      if (change > 0.0) {
4062        // going down
4063        if (value <= abcLower_[iPivot] + primalTolerance_) {
4064          if (iPivot == sequenceOut_ && value > abcLower_[iPivot] - relaxedTolerance)
4065            value = abcLower_[iPivot];
4066          abcNonLinearCost_->setOneBasic(iRow, value);
4067        }
4068      } else {
4069        // going up
4070        if (value >= abcUpper_[iPivot] - primalTolerance_) {
4071          if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
4072            value = abcUpper_[iPivot];
4073          abcNonLinearCost_->setOneBasic(iRow, value);
4074        }
4075      }
4076    }
4077  }
4078  //rowArray.setNumElements(0);
4079}
4080/* After rowArray will have cost changes for use next major iteration
4081 */
4082void
4083AbcSimplexPrimal::createUpdateDuals(CoinIndexedVector & rowArray,
4084                                    const double * originalCost,
4085                                    const double extraCost[4],
4086                                    double & /*objectiveChange*/,
4087                                    int /*valuesPass*/)
4088{
4089  int number=0;
4090  double * work = rowArray.denseVector();
4091  int * which = rowArray.getIndices();
4092  for (int iRow=0;iRow<numberRows_;iRow++) {
4093    if (originalCost[iRow]!=costBasic_[iRow]) {
4094      work[iRow]=costBasic_[iRow]-originalCost[iRow];
4095      which[number++]=iRow;
4096    }
4097  }
4098  rowArray.setNumElements(number);
4099}
4100// Update minor candidate vector
4101double
4102AbcSimplexPrimal::updateMinorCandidate(const CoinIndexedVector & updateBy,
4103                                       CoinIndexedVector & candidate,
4104                                       int sequenceIn)
4105{
4106  const double * COIN_RESTRICT regionUpdate = updateBy.denseVector (  );
4107  const int * COIN_RESTRICT regionIndexUpdate = updateBy.getIndices (  );
4108  int numberNonZeroUpdate = updateBy.getNumElements (  );
4109  double * COIN_RESTRICT regionCandidate = candidate.denseVector (  );
4110  int * COIN_RESTRICT regionIndexCandidate = candidate.getIndices (  );
4111  int numberNonZeroCandidate = candidate.getNumElements (  );
4112  assert (fabs(alpha_-regionUpdate[pivotRow_])<1.0e-5);
4113  double value=regionCandidate[pivotRow_];
4114  if (value) {
4115    value /= alpha_;
4116    for (int i=0;i<numberNonZeroUpdate;i++) {
4117      int iRow=regionIndexUpdate[i];
4118      double oldValue = regionCandidate[iRow];
4119      double newValue = oldValue-value*regionUpdate[iRow];
4120      if (oldValue) {
4121        if (!newValue)
4122          newValue=COIN_INDEXED_REALLY_TINY_ELEMENT;
4123        regionCandidate[iRow]=newValue;
4124      } else {
4125        assert (newValue);
4126        regionCandidate[iRow]=newValue;
4127        regionIndexCandidate[numberNonZeroCandidate++]=iRow;
4128      }
4129    }
4130    double oldValue = regionCandidate[pivotRow_];
4131    double newValue = value;
4132    if (oldValue) {
4133      if (!newValue)
4134        newValue=COIN_INDEXED_REALLY_TINY_ELEMENT;
4135      regionCandidate[pivotRow_]=newValue;
4136    } else {
4137      assert (newValue);
4138      regionCandidate[pivotRow_]=newValue;
4139      regionIndexCandidate[numberNonZeroCandidate++]=pivotRow_;
4140    }
4141    candidate.setNumElements(numberNonZeroCandidate);
4142  }
4143  double newDj=abcCost_[sequenceIn];
4144  for (int i = 0; i < numberNonZeroCandidate; i++) {
4145    int iRow = regionIndexCandidate[i];
4146    double alpha = regionCandidate[iRow];
4147    newDj -= alpha * costBasic_[iRow];
4148  }
4149  return newDj;
4150}
4151// Update partial Ftran by R update
4152void
4153AbcSimplexPrimal::updatePartialUpdate(CoinIndexedVector & partialUpdate)
4154{
4155  CoinAbcFactorization * factorization=
4156    dynamic_cast<CoinAbcFactorization *>(abcFactorization_->factorization());
4157  if (factorization)
4158    factorization->updatePartialUpdate(partialUpdate);
4159}
4160// Create primal ray
4161void
4162AbcSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
4163{
4164  delete [] ray_;
4165  ray_ = new double [numberColumns_];
4166  CoinZeroN(ray_, numberColumns_);
4167  int number = rowArray->getNumElements();
4168  int * index = rowArray->getIndices();
4169  double * array = rowArray->denseVector();
4170  double way = -directionIn_;
4171  int i;
4172  double zeroTolerance = 1.0e-12;
4173  if (sequenceIn_ < numberColumns_)
4174    ray_[sequenceIn_] = directionIn_;
4175  for (i = 0; i < number; i++) {
4176    int iRow = index[i];
4177    int iPivot = abcPivotVariable_[iRow];
4178    double arrayValue = array[iRow];
4179    if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4180      ray_[iPivot] = way * arrayValue;
4181  }
4182}
4183/* Get next superbasic -1 if none,
4184   Normal type is 1
4185   If type is 3 then initializes sorted list
4186   if 2 uses list.
4187*/
4188int
4189AbcSimplexPrimal::nextSuperBasic(int superBasicType,
4190                                 CoinIndexedVector * columnArray)
4191{
4192  int returnValue = -1;
4193  bool finished = false;
4194  while (!finished) {
4195    returnValue = firstFree_;
4196    int iColumn = firstFree_ + 1;
4197    if (superBasicType > 1) {
4198      if (superBasicType > 2) {
4199        // Initialize list
4200        // Wild guess that lower bound more natural than upper
4201        int number = 0;
4202        double * work = columnArray->denseVector();
4203        int * which = columnArray->getIndices();
4204        for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
4205          if (!flagged(iColumn)) {
4206            if (getInternalStatus(iColumn) == superBasic) {
4207              if (fabs(abcSolution_[iColumn] - abcLower_[iColumn]) <= primalTolerance_) {
4208                abcSolution_[iColumn] = abcLower_[iColumn];
4209                setInternalStatus(iColumn, atLowerBound);
4210              } else if (fabs(abcSolution_[iColumn] - abcUpper_[iColumn])
4211                         <= primalTolerance_) {
4212                abcSolution_[iColumn] = abcUpper_[iColumn];
4213                setInternalStatus(iColumn, atUpperBound);
4214              } else if (abcLower_[iColumn] < -1.0e20 && abcUpper_[iColumn] > 1.0e20) {
4215                setInternalStatus(iColumn, isFree);
4216                break;
4217              } else if (!flagged(iColumn)) {
4218                // put ones near bounds at end after sorting
4219                work[number] = - CoinMin(0.1 * (abcSolution_[iColumn] - abcLower_[iColumn]),
4220                                         abcUpper_[iColumn] - abcSolution_[iColumn]);
4221                which[number++] = iColumn;
4222              }
4223            }
4224          }
4225        }
4226        CoinSort_2(work, work + number, which);
4227        columnArray->setNumElements(number);
4228        CoinZeroN(work, number);
4229      }
4230      int * which = columnArray->getIndices();
4231      int number = columnArray->getNumElements();
4232      if (!number) {
4233        // finished
4234        iColumn = numberRows_ + numberColumns_;
4235        returnValue = -1;
4236      } else {
4237        number--;
4238        returnValue = which[number];
4239        iColumn = returnValue;
4240        columnArray->setNumElements(number);
4241      }
4242    } else {
4243      for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
4244        if (!flagged(iColumn)) {
4245          if (getInternalStatus(iColumn) == superBasic||
4246              getInternalStatus(iColumn) == isFree) {
4247            if (fabs(abcSolution_[iColumn] - abcLower_[iColumn]) <= primalTolerance_) {
4248              abcSolution_[iColumn] = abcLower_[iColumn];
4249              setInternalStatus(iColumn, atLowerBound);
4250            } else if (fabs(abcSolution_[iColumn] - abcUpper_[iColumn])
4251                       <= primalTolerance_) {
4252              abcSolution_[iColumn] = abcUpper_[iColumn];
4253              setInternalStatus(iColumn, atUpperBound);
4254            } else if (abcLower_[iColumn] < -1.0e20 && abcUpper_[iColumn] > 1.0e20) {
4255              setInternalStatus(iColumn, isFree);
4256              if (fabs(abcDj_[iColumn])>10.0*dualTolerance_)
4257                break;
4258            } else {
4259              break;
4260            }
4261          }
4262        }
4263      }
4264    }
4265    firstFree_ = iColumn;
4266    finished = true;
4267    if (firstFree_ == numberRows_ + numberColumns_)
4268      firstFree_ = -1;
4269    if (returnValue >= 0 && getInternalStatus(returnValue) != superBasic && getInternalStatus(returnValue) != isFree)
4270      finished = false; // somehow picked up odd one
4271  }
4272  return returnValue;
4273}
4274void
4275AbcSimplexPrimal::clearAll()
4276{
4277  int number = usefulArray_[arrayForFtran_].getNumElements();
4278  int * which = usefulArray_[arrayForFtran_].getIndices();
4279
4280  int iIndex;
4281  for (iIndex = 0; iIndex < number; iIndex++) {
4282
4283    int iRow = which[iIndex];
4284    clearActive(iRow);
4285  }
4286  usefulArray_[arrayForFtran_].clear();
4287  for (int i=0;i<6;i++) {
4288    if (rowArray_[i])
4289      rowArray_[i]->clear();
4290    if (columnArray_[i])
4291      columnArray_[i]->clear();
4292  }
4293}
Note: See TracBrowser for help on using the repository browser.