source: trunk/Clp/src/ClpSimplexPrimal.cpp @ 1723

Last change on this file since 1723 was 1723, checked in by forrest, 9 years ago

out some printf statements

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