# source:trunk/Clp/src/ClpSimplexPrimal.cpp@2030

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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