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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 164.0 KB
Line 
1/* $Id: ClpSimplexPrimal.cpp 1665 2011-01-04 17:55:54Z lou $ */
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                         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                    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                         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                         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)
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 ", 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 ", 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          double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
1442          if (fabs(objVal - objectiveValue_) > tol) {
1443#ifdef COIN_DEVELOP
1444               if (handler_->logLevel() > 0)
1445                    printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
1446                           objVal, objectiveValue_);
1447#endif
1448               objectiveValue_ = objVal;
1449          }
1450     }
1451     // save extra stuff
1452     matrix_->generalExpanded(this, 5, dummy);
1453     if (type == 0 || type == 1) {
1454          if (type != 1 || !saveStatus_) {
1455               // create save arrays
1456               delete [] saveStatus_;
1457               delete [] savedSolution_;
1458               saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1459               savedSolution_ = new double [numberRows_+numberColumns_];
1460          }
1461          // save arrays
1462          CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
1463          CoinMemcpyN(rowActivityWork_,
1464                      numberRows_, savedSolution_ + numberColumns_);
1465          CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
1466     }
1467     // see if in Cbc etc
1468     bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
1469     bool disaster = false;
1470     if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
1471          disasterArea_->saveInfo();
1472          disaster = true;
1473     }
1474     if (disaster)
1475          problemStatus_ = 3;
1476     if (problemStatus_ < 0 && !changeMade_) {
1477          problemStatus_ = 4; // unknown
1478     }
1479     lastGoodIteration_ = numberIterations_;
1480     if (numberIterations_ > lastBadIteration_ + 100)
1481          moreSpecialOptions_ &= ~16; // clear check accuracy flag
1482     if (goToDual || (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6
1483                      && largestDualError_ > 1.0e6)) {
1484          problemStatus_ = 10; // try dual
1485          // See if second call
1486          if ((moreSpecialOptions_ & 256) != 0||nonLinearCost_->numberInfeasibilities()>1.0e2) {
1487               numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
1488               sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
1489               // say infeasible
1490               if (numberPrimalInfeasibilities_)
1491                    problemStatus_ = 1;
1492          }
1493     }
1494     // make sure first free monotonic
1495     if (firstFree_ >= 0 && saveFirstFree >= 0) {
1496          firstFree_ = (numberIterations_) ? saveFirstFree : -1;
1497          nextSuperBasic(1, NULL);
1498     }
1499     if (doFactorization) {
1500          // restore weights (if saved) - also recompute infeasibility list
1501          if (tentativeStatus > -3)
1502               primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
1503          else
1504               primalColumnPivot_->saveWeights(this, 3);
1505          if (saveThreshold) {
1506               // use default at present
1507               factorization_->sparseThreshold(0);
1508               factorization_->goSparse();
1509          }
1510     }
1511     // Allow matrices to be sorted etc
1512     int fake = -999; // signal sort
1513     matrix_->correctSequence(this, fake, fake);
1514}
1515/*
1516   Row array has pivot column
1517   This chooses pivot row.
1518   For speed, we may need to go to a bucket approach when many
1519   variables go through bounds
1520   On exit rhsArray will have changes in costs of basic variables
1521*/
1522void
1523ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
1524                            CoinIndexedVector * rhsArray,
1525                            CoinIndexedVector * spareArray,
1526                            int valuesPass)
1527{
1528     double saveDj = dualIn_;
1529     if (valuesPass && objective_->type() < 2) {
1530          dualIn_ = cost_[sequenceIn_];
1531
1532          double * work = rowArray->denseVector();
1533          int number = rowArray->getNumElements();
1534          int * which = rowArray->getIndices();
1535
1536          int iIndex;
1537          for (iIndex = 0; iIndex < number; iIndex++) {
1538
1539               int iRow = which[iIndex];
1540               double alpha = work[iIndex];
1541               int iPivot = pivotVariable_[iRow];
1542               dualIn_ -= alpha * cost_[iPivot];
1543          }
1544          // determine direction here
1545          if (dualIn_ < -dualTolerance_) {
1546               directionIn_ = 1;
1547          } else if (dualIn_ > dualTolerance_) {
1548               directionIn_ = -1;
1549          } else {
1550               // towards nearest bound
1551               if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
1552                    directionIn_ = -1;
1553                    dualIn_ = dualTolerance_;
1554               } else {
1555                    directionIn_ = 1;
1556                    dualIn_ = -dualTolerance_;
1557               }
1558          }
1559     }
1560
1561     // sequence stays as row number until end
1562     pivotRow_ = -1;
1563     int numberRemaining = 0;
1564
1565     double totalThru = 0.0; // for when variables flip
1566     // Allow first few iterations to take tiny
1567     double acceptablePivot = 1.0e-1 * acceptablePivot_;
1568     if (numberIterations_ > 100)
1569          acceptablePivot = acceptablePivot_;
1570     if (factorization_->pivots() > 10)
1571          acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1572     else if (factorization_->pivots() > 5)
1573          acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1574     else if (factorization_->pivots())
1575          acceptablePivot = acceptablePivot_; // relax
1576     double bestEverPivot = acceptablePivot;
1577     int lastPivotRow = -1;
1578     double lastPivot = 0.0;
1579     double lastTheta = 1.0e50;
1580
1581     // use spareArrays to put ones looked at in
1582     // First one is list of candidates
1583     // We could compress if we really know we won't need any more
1584     // Second array has current set of pivot candidates
1585     // with a backup list saved in double * part of indexed vector
1586
1587     // pivot elements
1588     double * spare;
1589     // indices
1590     int * index;
1591     spareArray->clear();
1592     spare = spareArray->denseVector();
1593     index = spareArray->getIndices();
1594
1595     // we also need somewhere for effective rhs
1596     double * rhs = rhsArray->denseVector();
1597     // and we can use indices to point to alpha
1598     // that way we can store fabs(alpha)
1599     int * indexPoint = rhsArray->getIndices();
1600     //int numberFlip=0; // Those which may change if flips
1601
1602     /*
1603       First we get a list of possible pivots.  We can also see if the
1604       problem looks unbounded.
1605
1606       At first we increase theta and see what happens.  We start
1607       theta at a reasonable guess.  If in right area then we do bit by bit.
1608       We save possible pivot candidates
1609
1610      */
1611
1612     // do first pass to get possibles
1613     // We can also see if unbounded
1614
1615     double * work = rowArray->denseVector();
1616     int number = rowArray->getNumElements();
1617     int * which = rowArray->getIndices();
1618
1619     // we need to swap sign if coming in from ub
1620     double way = directionIn_;
1621     double maximumMovement;
1622     if (way > 0.0)
1623          maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
1624     else
1625          maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
1626
1627     double averageTheta = nonLinearCost_->averageTheta();
1628     double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
1629     double upperTheta = maximumMovement;
1630     if (tentativeTheta > 0.5 * maximumMovement)
1631          tentativeTheta = maximumMovement;
1632     bool thetaAtMaximum = tentativeTheta == maximumMovement;
1633     // In case tiny bounds increase
1634     if (maximumMovement < 1.0)
1635          tentativeTheta *= 1.1;
1636     double dualCheck = fabs(dualIn_);
1637     // but make a bit more pessimistic
1638     dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
1639
1640     int iIndex;
1641     int pivotOne = -1;
1642     //#define CLP_DEBUG
1643#ifdef CLP_DEBUG
1644     if (numberIterations_ == -3839 || numberIterations_ == -3840) {
1645          double dj = cost_[sequenceIn_];
1646          printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
1647          for (iIndex = 0; iIndex < number; iIndex++) {
1648
1649               int iRow = which[iIndex];
1650               double alpha = work[iIndex];
1651               int iPivot = pivotVariable_[iRow];
1652               dj -= alpha * cost_[iPivot];
1653               printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1654                      iRow, iPivot, lower_[iPivot], solution_[iPivot], upper_[iPivot],
1655                      alpha, solution_[iPivot] - 1.0e9 * alpha, cost_[iPivot], dj);
1656          }
1657     }
1658#endif
1659     while (true) {
1660          pivotOne = -1;
1661          totalThru = 0.0;
1662          // We also re-compute reduced cost
1663          numberRemaining = 0;
1664          dualIn_ = cost_[sequenceIn_];
1665#ifndef NDEBUG
1666          double tolerance = primalTolerance_ * 1.002;
1667#endif
1668          for (iIndex = 0; iIndex < number; iIndex++) {
1669
1670               int iRow = which[iIndex];
1671               double alpha = work[iIndex];
1672               int iPivot = pivotVariable_[iRow];
1673               if (cost_[iPivot])
1674                    dualIn_ -= alpha * cost_[iPivot];
1675               alpha *= way;
1676               double oldValue = solution_[iPivot];
1677               // get where in bound sequence
1678               // note that after this alpha is actually fabs(alpha)
1679               bool possible;
1680               // do computation same way as later on in primal
1681               if (alpha > 0.0) {
1682                    // basic variable going towards lower bound
1683                    double bound = lower_[iPivot];
1684                    // must be exactly same as when used
1685                    double change = tentativeTheta * alpha;
1686                    possible = (oldValue - change) <= bound + primalTolerance_;
1687                    oldValue -= bound;
1688               } else {
1689                    // basic variable going towards upper bound
1690                    double bound = upper_[iPivot];
1691                    // must be exactly same as when used
1692                    double change = tentativeTheta * alpha;
1693                    possible = (oldValue - change) >= bound - primalTolerance_;
1694                    oldValue = bound - oldValue;
1695                    alpha = - alpha;
1696               }
1697               double value;
1698               assert (oldValue >= -tolerance);
1699               if (possible) {
1700                    value = oldValue - upperTheta * alpha;
1701#ifdef CLP_USER_DRIVEN1
1702                    int sequenceOut=pivotVariable_[index[iIndex]];
1703                    if(!userChoiceValid1(this,sequenceOut,oldValue,
1704                                         upperTheta,alpha,work[iIndex]*way))
1705                      value =0.0; // say can't use
1706#endif
1707                    if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1708                         upperTheta = (oldValue + primalTolerance_) / alpha;
1709                         pivotOne = numberRemaining;
1710                    }
1711                    // add to list
1712                    spare[numberRemaining] = alpha;
1713                    rhs[numberRemaining] = oldValue;
1714                    indexPoint[numberRemaining] = iIndex;
1715                    index[numberRemaining++] = iRow;
1716                    totalThru += alpha;
1717                    setActive(iRow);
1718                    //} else if (value<primalTolerance_*1.002) {
1719                    // May change if is a flip
1720                    //indexRhs[numberFlip++]=iRow;
1721               }
1722          }
1723          if (upperTheta < maximumMovement && totalThru*infeasibilityCost_ >= 1.0001 * dualCheck) {
1724               // Can pivot here
1725               break;
1726          } else if (!thetaAtMaximum) {
1727               //printf("Going round with average theta of %g\n",averageTheta);
1728               tentativeTheta = maximumMovement;
1729               thetaAtMaximum = true; // seems to be odd compiler error
1730          } else {
1731               break;
1732          }
1733     }
1734     totalThru = 0.0;
1735
1736     theta_ = maximumMovement;
1737
1738     bool goBackOne = false;
1739     if (objective_->type() > 1)
1740          dualIn_ = saveDj;
1741
1742     //printf("%d remain out of %d\n",numberRemaining,number);
1743     int iTry = 0;
1744#define MAXTRY 1000
1745     if (numberRemaining && upperTheta < maximumMovement) {
1746          // First check if previously chosen one will work
1747          if (pivotOne >= 0 && 0) {
1748               double thruCost = infeasibilityCost_ * spare[pivotOne];
1749               if (thruCost >= 0.99 * fabs(dualIn_))
1750                    printf("Could pivot on %d as change %g dj %g\n",
1751                           index[pivotOne], thruCost, dualIn_);
1752               double alpha = spare[pivotOne];
1753               double oldValue = rhs[pivotOne];
1754               theta_ = oldValue / alpha;
1755               pivotRow_ = pivotOne;
1756               // Stop loop
1757               iTry = MAXTRY;
1758          }
1759
1760          // first get ratio with tolerance
1761          for ( ; iTry < MAXTRY; iTry++) {
1762
1763               upperTheta = maximumMovement;
1764               int iBest = -1;
1765               for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
1766
1767                    double alpha = spare[iIndex];
1768                    double oldValue = rhs[iIndex];
1769                    double value = oldValue - upperTheta * alpha;
1770
1771#ifdef CLP_USER_DRIVEN1
1772                    int sequenceOut=pivotVariable_[index[iIndex]];
1773                    if(!userChoiceValid1(this,sequenceOut,oldValue,
1774                                         upperTheta,alpha, 0.0))
1775                      value =0.0; // say can't use
1776#endif
1777                    if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1778                         upperTheta = (oldValue + primalTolerance_) / alpha;
1779                         iBest = iIndex; // just in case weird numbers
1780                    }
1781               }
1782
1783               // now look at best in this lot
1784               // But also see how infeasible small pivots will make
1785               double sumInfeasibilities = 0.0;
1786               double bestPivot = acceptablePivot;
1787               pivotRow_ = -1;
1788               for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
1789
1790                    int iRow = index[iIndex];
1791                    double alpha = spare[iIndex];
1792                    double oldValue = rhs[iIndex];
1793                    double value = oldValue - upperTheta * alpha;
1794
1795                    if (value <= 0 || iBest == iIndex) {
1796                         // how much would it cost to go thru and modify bound
1797                         double trueAlpha = way * work[indexPoint[iIndex]];
1798                         totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow], trueAlpha, rhs[iIndex]);
1799                         setActive(iRow);
1800                         if (alpha > bestPivot) {
1801                              bestPivot = alpha;
1802                              theta_ = oldValue / bestPivot;
1803                              pivotRow_ = iIndex;
1804                         } else if (alpha < acceptablePivot
1805#ifdef CLP_USER_DRIVEN1
1806                      ||!userChoiceValid1(this,pivotVariable_[index[iIndex]],
1807                                          oldValue,upperTheta,alpha,0.0)
1808#endif
1809                                    ) {
1810                              if (value < -primalTolerance_)
1811                                   sumInfeasibilities += -value - primalTolerance_;
1812                         }
1813                    }
1814               }
1815               if (bestPivot < 0.1 * bestEverPivot &&
1816                         bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
1817                    // back to previous one
1818                    goBackOne = true;
1819                    break;
1820               } else if (pivotRow_ == -1 && upperTheta > largeValue_) {
1821                    if (lastPivot > acceptablePivot) {
1822                         // back to previous one
1823                         goBackOne = true;
1824                    } else {
1825                         // can only get here if all pivots so far too small
1826                    }
1827                    break;
1828               } else if (totalThru >= dualCheck) {
1829                    if (sumInfeasibilities > primalTolerance_ && !nonLinearCost_->numberInfeasibilities()) {
1830                         // Looks a bad choice
1831                         if (lastPivot > acceptablePivot) {
1832                              goBackOne = true;
1833                         } else {
1834                              // say no good
1835                              dualIn_ = 0.0;
1836                         }
1837                    }
1838                    break; // no point trying another loop
1839               } else {
1840                    lastPivotRow = pivotRow_;
1841                    lastTheta = theta_;
1842                    if (bestPivot > bestEverPivot)
1843                         bestEverPivot = bestPivot;
1844               }
1845          }
1846          // can get here without pivotRow_ set but with lastPivotRow
1847          if (goBackOne || (pivotRow_ < 0 && lastPivotRow >= 0)) {
1848               // back to previous one
1849               pivotRow_ = lastPivotRow;
1850               theta_ = lastTheta;
1851          }
1852     } else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
1853          // looks unbounded
1854          valueOut_ = COIN_DBL_MAX; // say odd
1855          if (nonLinearCost_->numberInfeasibilities()) {
1856               // but infeasible??
1857               // move variable but don't pivot
1858               tentativeTheta = 1.0e50;
1859               for (iIndex = 0; iIndex < number; iIndex++) {
1860                    int iRow = which[iIndex];
1861                    double alpha = work[iIndex];
1862                    int iPivot = pivotVariable_[iRow];
1863                    alpha *= way;
1864                    double oldValue = solution_[iPivot];
1865                    // get where in bound sequence
1866                    // note that after this alpha is actually fabs(alpha)
1867                    if (alpha > 0.0) {
1868                         // basic variable going towards lower bound
1869                         double bound = lower_[iPivot];
1870                         oldValue -= bound;
1871                    } else {
1872                         // basic variable going towards upper bound
1873                         double bound = upper_[iPivot];
1874                         oldValue = bound - oldValue;
1875                         alpha = - alpha;
1876                    }
1877                    if (oldValue - tentativeTheta * alpha < 0.0) {
1878                         tentativeTheta = oldValue / alpha;
1879                    }
1880               }
1881               // If free in then see if we can get to 0.0
1882               if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
1883                    if (dualIn_ * valueIn_ > 0.0) {
1884                         if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) || tentativeTheta > 1.0e20)) {
1885                              tentativeTheta = fabs(valueIn_);
1886                         }
1887                    }
1888               }
1889               if (tentativeTheta < 1.0e10)
1890                    valueOut_ = valueIn_ + way * tentativeTheta;
1891          }
1892     }
1893     //if (iTry>50)
1894     //printf("** %d tries\n",iTry);
1895     if (pivotRow_ >= 0) {
1896          int position = pivotRow_; // position in list
1897          pivotRow_ = index[position];
1898          alpha_ = work[indexPoint[position]];
1899          // translate to sequence
1900          sequenceOut_ = pivotVariable_[pivotRow_];
1901          valueOut_ = solution(sequenceOut_);
1902          lowerOut_ = lower_[sequenceOut_];
1903          upperOut_ = upper_[sequenceOut_];
1904#define MINIMUMTHETA 1.0e-12
1905          // Movement should be minimum for anti-degeneracy - unless
1906          // fixed variable out
1907          double minimumTheta;
1908          if (upperOut_ > lowerOut_)
1909               minimumTheta = MINIMUMTHETA;
1910          else
1911               minimumTheta = 0.0;
1912          // But can't go infeasible
1913          double distance;
1914          if (alpha_ * way > 0.0)
1915               distance = valueOut_ - lowerOut_;
1916          else
1917               distance = upperOut_ - valueOut_;
1918          if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
1919               minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
1920          // will we need to increase tolerance
1921          //#define CLP_DEBUG
1922          double largestInfeasibility = primalTolerance_;
1923          if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
1924               theta_ = minimumTheta;
1925               for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
1926                    largestInfeasibility = CoinMax(largestInfeasibility,
1927                                                   -(rhs[iIndex] - spare[iIndex] * theta_));
1928               }
1929//#define CLP_DEBUG
1930#ifdef CLP_DEBUG
1931               if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
1932                    printf("Primal tolerance increased from %g to %g\n",
1933                           primalTolerance_, largestInfeasibility);
1934#endif
1935//#undef CLP_DEBUG
1936               primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
1937          }
1938          // Need to look at all in some cases
1939          if (theta_ > tentativeTheta) {
1940               for (iIndex = 0; iIndex < number; iIndex++)
1941                    setActive(which[iIndex]);
1942          }
1943          if (way < 0.0)
1944               theta_ = - theta_;
1945          double newValue = valueOut_ - theta_ * alpha_;
1946          // If  4 bit set - Force outgoing variables to exact bound (primal)
1947          if (alpha_ * way < 0.0) {
1948               directionOut_ = -1;    // to upper bound
1949               if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
1950                    upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
1951               } else {
1952                    upperOut_ = newValue;
1953               }
1954          } else {
1955               directionOut_ = 1;    // to lower bound
1956               if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
1957                    lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
1958               } else {
1959                    lowerOut_ = newValue;
1960               }
1961          }
1962          dualOut_ = reducedCost(sequenceOut_);
1963     } else if (maximumMovement < 1.0e20) {
1964          // flip
1965          pivotRow_ = -2; // so we can tell its a flip
1966          sequenceOut_ = sequenceIn_;
1967          valueOut_ = valueIn_;
1968          dualOut_ = dualIn_;
1969          lowerOut_ = lowerIn_;
1970          upperOut_ = upperIn_;
1971          alpha_ = 0.0;
1972          if (way < 0.0) {
1973               directionOut_ = 1;    // to lower bound
1974               theta_ = lowerOut_ - valueOut_;
1975          } else {
1976               directionOut_ = -1;    // to upper bound
1977               theta_ = upperOut_ - valueOut_;
1978          }
1979     }
1980
1981     double theta1 = CoinMax(theta_, 1.0e-12);
1982     double theta2 = numberIterations_ * nonLinearCost_->averageTheta();
1983     // Set average theta
1984     nonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast<double> (numberIterations_ + 1)));
1985     //if (numberIterations_%1000==0)
1986     //printf("average theta is %g\n",nonLinearCost_->averageTheta());
1987
1988     // clear arrays
1989
1990     CoinZeroN(spare, numberRemaining);
1991
1992     // put back original bounds etc
1993     CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
1994     rhsArray->setNumElements(numberRemaining);
1995     rhsArray->setPacked();
1996     nonLinearCost_->goBackAll(rhsArray);
1997     rhsArray->clear();
1998
1999}
2000/*
2001   Chooses primal pivot column
2002   updateArray has cost updates (also use pivotRow_ from last iteration)
2003   Would be faster with separate region to scan
2004   and will have this (with square of infeasibility) when steepest
2005   For easy problems we can just choose one of the first columns we look at
2006*/
2007void
2008ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
2009                               CoinIndexedVector * spareRow1,
2010                               CoinIndexedVector * spareRow2,
2011                               CoinIndexedVector * spareColumn1,
2012                               CoinIndexedVector * spareColumn2)
2013{
2014
2015     ClpMatrixBase * saveMatrix = matrix_;
2016     double * saveRowScale = rowScale_;
2017     if (scaledMatrix_) {
2018          rowScale_ = NULL;
2019          matrix_ = scaledMatrix_;
2020     }
2021     sequenceIn_ = primalColumnPivot_->pivotColumn(updates, spareRow1,
2022                   spareRow2, spareColumn1,
2023                   spareColumn2);
2024     if (scaledMatrix_) {
2025          matrix_ = saveMatrix;
2026          rowScale_ = saveRowScale;
2027     }
2028     if (sequenceIn_ >= 0) {
2029          valueIn_ = solution_[sequenceIn_];
2030          dualIn_ = dj_[sequenceIn_];
2031          if (nonLinearCost_->lookBothWays()) {
2032               // double check
2033               ClpSimplex::Status status = getStatus(sequenceIn_);
2034
2035               switch(status) {
2036               case ClpSimplex::atUpperBound:
2037                    if (dualIn_ < 0.0) {
2038                         // move to other side
2039                         printf("For %d U (%g, %g, %g) dj changed from %g",
2040                                sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
2041                                upper_[sequenceIn_], dualIn_);
2042                         dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
2043                         printf(" to %g\n", dualIn_);
2044                         nonLinearCost_->setOne(sequenceIn_, upper_[sequenceIn_] + 2.0 * currentPrimalTolerance());
2045                         setStatus(sequenceIn_, ClpSimplex::atLowerBound);
2046                    }
2047                    break;
2048               case ClpSimplex::atLowerBound:
2049                    if (dualIn_ > 0.0) {
2050                         // move to other side
2051                         printf("For %d L (%g, %g, %g) dj changed from %g",
2052                                sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
2053                                upper_[sequenceIn_], dualIn_);
2054                         dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
2055                         printf(" to %g\n", dualIn_);
2056                         nonLinearCost_->setOne(sequenceIn_, lower_[sequenceIn_] - 2.0 * currentPrimalTolerance());
2057                         setStatus(sequenceIn_, ClpSimplex::atUpperBound);
2058                    }
2059                    break;
2060               default:
2061                    break;
2062               }
2063          }
2064          lowerIn_ = lower_[sequenceIn_];
2065          upperIn_ = upper_[sequenceIn_];
2066          if (dualIn_ > 0.0)
2067               directionIn_ = -1;
2068          else
2069               directionIn_ = 1;
2070     } else {
2071          sequenceIn_ = -1;
2072     }
2073}
2074/* The primals are updated by the given array.
2075   Returns number of infeasibilities.
2076   After rowArray will have list of cost changes
2077*/
2078int
2079ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
2080                                        double theta,
2081                                        double & objectiveChange,
2082                                        int valuesPass)
2083{
2084     // Cost on pivot row may change - may need to change dualIn
2085     double oldCost = 0.0;
2086     if (pivotRow_ >= 0)
2087          oldCost = cost_[sequenceOut_];
2088     //rowArray->scanAndPack();
2089     double * work = rowArray->denseVector();
2090     int number = rowArray->getNumElements();
2091     int * which = rowArray->getIndices();
2092
2093     int newNumber = 0;
2094     int pivotPosition = -1;
2095     nonLinearCost_->setChangeInCost(0.0);
2096     //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
2097     //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
2098     // allow for case where bound+tolerance == bound
2099     //double tolerance = 0.999999*primalTolerance_;
2100     double relaxedTolerance = 1.001 * primalTolerance_;
2101     int iIndex;
2102     if (!valuesPass) {
2103          for (iIndex = 0; iIndex < number; iIndex++) {
2104
2105               int iRow = which[iIndex];
2106               double alpha = work[iIndex];
2107               work[iIndex] = 0.0;
2108               int iPivot = pivotVariable_[iRow];
2109               double change = theta * alpha;
2110               double value = solution_[iPivot] - change;
2111               solution_[iPivot] = value;
2112#ifndef NDEBUG
2113               // check if not active then okay
2114               if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
2115                    // But make sure one going out is feasible
2116                    if (change > 0.0) {
2117                         // going down
2118                         if (value <= lower_[iPivot] + primalTolerance_) {
2119                              if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2120                                   value = lower_[iPivot];
2121                              double difference = nonLinearCost_->setOne(iPivot, value);
2122                              assert (!difference || fabs(change) > 1.0e9);
2123                         }
2124                    } else {
2125                         // going up
2126                         if (value >= upper_[iPivot] - primalTolerance_) {
2127                              if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2128                                   value = upper_[iPivot];
2129                              double difference = nonLinearCost_->setOne(iPivot, value);
2130                              assert (!difference || fabs(change) > 1.0e9);
2131                         }
2132                    }
2133               }
2134#endif
2135               if (active(iRow) || theta_ < 0.0) {
2136                    clearActive(iRow);
2137                    // But make sure one going out is feasible
2138                    if (change > 0.0) {
2139                         // going down
2140                         if (value <= lower_[iPivot] + primalTolerance_) {
2141                              if (iPivot == sequenceOut_ && value >= lower_[iPivot] - relaxedTolerance)
2142                                   value = lower_[iPivot];
2143                              double difference = nonLinearCost_->setOne(iPivot, value);
2144                              if (difference) {
2145                                   if (iRow == pivotRow_)
2146                                        pivotPosition = newNumber;
2147                                   work[newNumber] = difference;
2148                                   //change reduced cost on this
2149                                   dj_[iPivot] = -difference;
2150                                   which[newNumber++] = iRow;
2151                              }
2152                         }
2153                    } else {
2154                         // going up
2155                         if (value >= upper_[iPivot] - primalTolerance_) {
2156                              if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2157                                   value = upper_[iPivot];
2158                              double difference = nonLinearCost_->setOne(iPivot, value);
2159                              if (difference) {
2160                                   if (iRow == pivotRow_)
2161                                        pivotPosition = newNumber;
2162                                   work[newNumber] = difference;
2163                                   //change reduced cost on this
2164                                   dj_[iPivot] = -difference;
2165                                   which[newNumber++] = iRow;
2166                              }
2167                         }
2168                    }
2169               }
2170          }
2171     } else {
2172          // values pass so look at all
2173          for (iIndex = 0; iIndex < number; iIndex++) {
2174
2175               int iRow = which[iIndex];
2176               double alpha = work[iIndex];
2177               work[iIndex] = 0.0;
2178               int iPivot = pivotVariable_[iRow];
2179               double change = theta * alpha;
2180               double value = solution_[iPivot] - change;
2181               solution_[iPivot] = value;
2182               clearActive(iRow);
2183               // But make sure one going out is feasible
2184               if (change > 0.0) {
2185                    // going down
2186                    if (value <= lower_[iPivot] + primalTolerance_) {
2187                         if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2188                              value = lower_[iPivot];
2189                         double difference = nonLinearCost_->setOne(iPivot, value);
2190                         if (difference) {
2191                              if (iRow == pivotRow_)
2192                                   pivotPosition = newNumber;
2193                              work[newNumber] = difference;
2194                              //change reduced cost on this
2195                              dj_[iPivot] = -difference;
2196                              which[newNumber++] = iRow;
2197                         }
2198                    }
2199               } else {
2200                    // going up
2201                    if (value >= upper_[iPivot] - primalTolerance_) {
2202                         if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2203                              value = upper_[iPivot];
2204                         double difference = nonLinearCost_->setOne(iPivot, value);
2205                         if (difference) {
2206                              if (iRow == pivotRow_)
2207                                   pivotPosition = newNumber;
2208                              work[newNumber] = difference;
2209                              //change reduced cost on this
2210                              dj_[iPivot] = -difference;
2211                              which[newNumber++] = iRow;
2212                         }
2213                    }
2214               }
2215          }
2216     }
2217     objectiveChange += nonLinearCost_->changeInCost();
2218     rowArray->setPacked();
2219#if 0
2220     rowArray->setNumElements(newNumber);
2221     rowArray->expand();
2222     if (pivotRow_ >= 0) {
2223          dualIn_ += (oldCost - cost_[sequenceOut_]);
2224          // update change vector to include pivot
2225          rowArray->add(pivotRow_, -dualIn_);
2226          // and convert to packed
2227          rowArray->scanAndPack();
2228     } else {
2229          // and convert to packed
2230          rowArray->scanAndPack();
2231     }
2232#else
2233     if (pivotRow_ >= 0) {
2234          double dualIn = dualIn_ + (oldCost - cost_[sequenceOut_]);
2235          // update change vector to include pivot
2236          if (pivotPosition >= 0) {
2237               work[pivotPosition] -= dualIn;
2238          } else {
2239               work[newNumber] = -dualIn;
2240               which[newNumber++] = pivotRow_;
2241          }
2242     }
2243     rowArray->setNumElements(newNumber);
2244#endif
2245     return 0;
2246}
2247// Perturbs problem
2248void
2249ClpSimplexPrimal::perturb(int type)
2250{
2251     if (perturbation_ > 100)
2252          return; //perturbed already
2253     if (perturbation_ == 100)
2254          perturbation_ = 50; // treat as normal
2255     int savePerturbation = perturbation_;
2256     int i;
2257     if (!numberIterations_)
2258          cleanStatus(); // make sure status okay
2259     // Make sure feasible bounds
2260     if (nonLinearCost_)
2261          nonLinearCost_->feasibleBounds();
2262     // look at element range
2263     double smallestNegative;
2264     double largestNegative;
2265     double smallestPositive;
2266     double largestPositive;
2267     matrix_->rangeOfElements(smallestNegative, largestNegative,
2268                              smallestPositive, largestPositive);
2269     smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
2270     largestPositive = CoinMax(fabs(largestNegative), largestPositive);
2271     double elementRatio = largestPositive / smallestPositive;
2272     if (!numberIterations_ && perturbation_ == 50) {
2273          // See if we need to perturb
2274          int numberTotal = CoinMax(numberRows_, numberColumns_);
2275          double * sort = new double[numberTotal];
2276          int nFixed = 0;
2277          for (i = 0; i < numberRows_; i++) {
2278               double lo = fabs(rowLower_[i]);
2279               double up = fabs(rowUpper_[i]);
2280               double value = 0.0;
2281               if (lo && lo < 1.0e20) {
2282                    if (up && up < 1.0e20) {
2283                         value = 0.5 * (lo + up);
2284                         if (lo == up)
2285                              nFixed++;
2286                    } else {
2287                         value = lo;
2288                    }
2289               } else {
2290                    if (up && up < 1.0e20)
2291                         value = up;
2292               }
2293               sort[i] = value;
2294          }
2295          std::sort(sort, sort + numberRows_);
2296          int number = 1;
2297          double last = sort[0];
2298          for (i = 1; i < numberRows_; i++) {
2299               if (last != sort[i])
2300                    number++;
2301               last = sort[i];
2302          }
2303#ifdef KEEP_GOING_IF_FIXED
2304          //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
2305          //   numberRows_,nFixed,elementRatio);
2306#endif
2307          if (number * 4 > numberRows_ || elementRatio > 1.0e12) {
2308               perturbation_ = 100;
2309               delete [] sort;
2310               return; // good enough
2311          }
2312          number = 0;
2313#ifdef KEEP_GOING_IF_FIXED
2314          if (!integerType_) {
2315               // look at columns
2316               nFixed = 0;
2317               for (i = 0; i < numberColumns_; i++) {
2318                    double lo = fabs(columnLower_[i]);
2319                    double up = fabs(columnUpper_[i]);
2320                    double value = 0.0;
2321                    if (lo && lo < 1.0e20) {
2322                         if (up && up < 1.0e20) {
2323                              value = 0.5 * (lo + up);
2324                              if (lo == up)
2325                                   nFixed++;
2326                         } else {
2327                              value = lo;
2328                         }
2329                    } else {
2330                         if (up && up < 1.0e20)
2331                              value = up;
2332                    }
2333                    sort[i] = value;
2334               }
2335               std::sort(sort, sort + numberColumns_);
2336               number = 1;
2337               last = sort[0];
2338               for (i = 1; i < numberColumns_; i++) {
2339                    if (last != sort[i])
2340                         number++;
2341                    last = sort[i];
2342               }
2343               //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
2344               //     numberColumns_,nFixed);
2345          }
2346#endif
2347          delete [] sort;
2348          if (number * 4 > numberColumns_) {
2349               perturbation_ = 100;
2350               return; // good enough
2351          }
2352     }
2353     // primal perturbation
2354     double perturbation = 1.0e-20;
2355     double bias = 1.0;
2356     int numberNonZero = 0;
2357     // maximum fraction of rhs/bounds to perturb
2358     double maximumFraction = 1.0e-5;
2359     if (perturbation_ >= 50) {
2360          perturbation = 1.0e-4;
2361          for (i = 0; i < numberColumns_ + numberRows_; i++) {
2362               if (upper_[i] > lower_[i] + primalTolerance_) {
2363                    double lowerValue, upperValue;
2364                    if (lower_[i] > -1.0e20)
2365                         lowerValue = fabs(lower_[i]);
2366                    else
2367                         lowerValue = 0.0;
2368                    if (upper_[i] < 1.0e20)
2369                         upperValue = fabs(upper_[i]);
2370                    else
2371                         upperValue = 0.0;
2372                    double value = CoinMax(fabs(lowerValue), fabs(upperValue));
2373                    value = CoinMin(value, upper_[i] - lower_[i]);
2374#if 1
2375                    if (value) {
2376                         perturbation += value;
2377                         numberNonZero++;
2378                    }
2379#else
2380                    perturbation = CoinMax(perturbation, value);
2381#endif
2382               }
2383          }
2384          if (numberNonZero)
2385               perturbation /= static_cast<double> (numberNonZero);
2386          else
2387               perturbation = 1.0e-1;
2388          if (perturbation_ > 50 && perturbation_ < 55) {
2389               // reduce
2390               while (perturbation_ > 50) {
2391                    perturbation_--;
2392                    perturbation *= 0.25;
2393                    bias *= 0.25;
2394               }
2395          } else if (perturbation_ >= 55 && perturbation_ < 60) {
2396               // increase
2397               while (perturbation_ > 55) {
2398                    perturbation_--;
2399                    perturbation *= 4.0;
2400               }
2401               perturbation_ = 50;
2402          }
2403     } else if (perturbation_ < 100) {
2404          perturbation = pow(10.0, perturbation_);
2405          // user is in charge
2406          maximumFraction = 1.0;
2407     }
2408     double largestZero = 0.0;
2409     double largest = 0.0;
2410     double largestPerCent = 0.0;
2411     bool printOut = (handler_->logLevel() == 63);
2412     printOut = false; //off
2413     // Check if all slack
2414     int number = 0;
2415     int iSequence;
2416     for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2417          if (getRowStatus(iSequence) == basic)
2418               number++;
2419     }
2420     if (rhsScale_ > 100.0) {
2421          // tone down perturbation
2422          maximumFraction *= 0.1;
2423     }
2424     if (number != numberRows_)
2425          type = 1;
2426     // modify bounds
2427     // Change so at least 1.0e-5 and no more than 0.1
2428     // For now just no more than 0.1
2429     // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
2430     // seems much slower???#define SAVE_PERT
2431#ifdef SAVE_PERT
2432     if (2 * numberColumns_ > maximumPerturbationSize_) {
2433          delete [] perturbationArray_;
2434          maximumPerturbationSize_ = 2 * numberColumns_;
2435          perturbationArray_ = new double [maximumPerturbationSize_];
2436          for (int iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
2437               perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
2438          }
2439     }
2440#endif
2441     if (type == 1) {
2442          double tolerance = 100.0 * primalTolerance_;
2443          //double multiplier = perturbation*maximumFraction;
2444          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2445               if (getStatus(iSequence) == basic) {
2446                    double lowerValue = lower_[iSequence];
2447                    double upperValue = upper_[iSequence];
2448                    if (upperValue > lowerValue + tolerance) {
2449                         double solutionValue = solution_[iSequence];
2450                         double difference = upperValue - lowerValue;
2451                         difference = CoinMin(difference, perturbation);
2452                         difference = CoinMin(difference, fabs(solutionValue) + 1.0);
2453                         double value = maximumFraction * (difference + bias);
2454                         value = CoinMin(value, 0.1);
2455#ifndef SAVE_PERT
2456                         value *= randomNumberGenerator_.randomDouble();
2457#else
2458                         value *= perturbationArray_[2*iSequence];
2459#endif
2460                         if (solutionValue - lowerValue <= primalTolerance_) {
2461                              lower_[iSequence] -= value;
2462                         } else if (upperValue - solutionValue <= primalTolerance_) {
2463                              upper_[iSequence] += value;
2464                         } else {
2465#if 0
2466                              if (iSequence >= numberColumns_) {
2467                                   // may not be at bound - but still perturb (unless free)
2468                                   if (upperValue > 1.0e30 && lowerValue < -1.0e30)
2469                                        value = 0.0;
2470                                   else
2471                                        value = - value; // as -1.0 in matrix
2472                              } else {
2473                                   value = 0.0;
2474                              }
2475#else
2476                              value = 0.0;
2477#endif
2478                         }
2479                         if (value) {
2480                              if (printOut)
2481                                   printf("col %d lower from %g to %g, upper from %g to %g\n",
2482                                          iSequence, lower_[iSequence], lowerValue, upper_[iSequence], upperValue);
2483                              if (solutionValue) {
2484                                   largest = CoinMax(largest, value);
2485                                   if (value > (fabs(solutionValue) + 1.0)*largestPerCent)
2486                                        largestPerCent = value / (fabs(solutionValue) + 1.0);
2487                              } else {
2488                                   largestZero = CoinMax(largestZero, value);
2489                              }
2490                         }
2491                    }
2492               }
2493          }
2494     } else {
2495          double tolerance = 100.0 * primalTolerance_;
2496          for (i = 0; i < numberColumns_; i++) {
2497               double lowerValue = lower_[i], upperValue = upper_[i];
2498               if (upperValue > lowerValue + primalTolerance_) {
2499                    double value = perturbation * maximumFraction;
2500                    value = CoinMin(value, 0.1);
2501#ifndef SAVE_PERT
2502                    value *= randomNumberGenerator_.randomDouble();
2503#else
2504                    value *= perturbationArray_[2*i+1];
2505#endif
2506                    value *= randomNumberGenerator_.randomDouble();
2507                    if (savePerturbation != 50) {
2508                         if (fabs(value) <= primalTolerance_)
2509                              value = 0.0;
2510                         if (lowerValue > -1.0e20 && lowerValue)
2511                              lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2512                         if (upperValue < 1.0e20 && upperValue)
2513                              upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2514                    } else if (value) {
2515                         double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2516                         // get in range
2517                         if (valueL <= tolerance) {
2518                              valueL *= 10.0;
2519                              while (valueL <= tolerance)
2520                                   valueL *= 10.0;
2521                         } else if (valueL > 1.0) {
2522                              valueL *= 0.1;
2523                              while (valueL > 1.0)
2524                                   valueL *= 0.1;
2525                         }
2526                         if (lowerValue > -1.0e20 && lowerValue)
2527                              lowerValue -= valueL;
2528                         double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2529                         // get in range
2530                         if (valueU <= tolerance) {
2531                              valueU *= 10.0;
2532                              while (valueU <= tolerance)
2533                                   valueU *= 10.0;
2534                         } else if (valueU > 1.0) {
2535                              valueU *= 0.1;
2536                              while (valueU > 1.0)
2537                                   valueU *= 0.1;
2538                         }
2539                         if (upperValue < 1.0e20 && upperValue)
2540                              upperValue += valueU;
2541                    }
2542                    if (lowerValue != lower_[i]) {
2543                         double difference = fabs(lowerValue - lower_[i]);
2544                         largest = CoinMax(largest, difference);
2545                         if (difference > fabs(lower_[i])*largestPerCent)
2546                              largestPerCent = fabs(difference / lower_[i]);
2547                    }
2548                    if (upperValue != upper_[i]) {
2549                         double difference = fabs(upperValue - upper_[i]);
2550                         largest = CoinMax(largest, difference);
2551                         if (difference > fabs(upper_[i])*largestPerCent)
2552                              largestPerCent = fabs(difference / upper_[i]);
2553                    }
2554                    if (printOut)
2555                         printf("col %d lower from %g to %g, upper from %g to %g\n",
2556                                i, lower_[i], lowerValue, upper_[i], upperValue);
2557               }
2558               lower_[i] = lowerValue;
2559               upper_[i] = upperValue;
2560          }
2561          for (; i < numberColumns_ + numberRows_; i++) {
2562               double lowerValue = lower_[i], upperValue = upper_[i];
2563               double value = perturbation * maximumFraction;
2564               value = CoinMin(value, 0.1);
2565               value *= randomNumberGenerator_.randomDouble();
2566               if (upperValue > lowerValue + tolerance) {
2567                    if (savePerturbation != 50) {
2568                         if (fabs(value) <= primalTolerance_)
2569                              value = 0.0;
2570                         if (lowerValue > -1.0e20 && lowerValue)
2571                              lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2572                         if (upperValue < 1.0e20 && upperValue)
2573                              upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2574                    } else if (value) {
2575                         double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2576                         // get in range
2577                         if (valueL <= tolerance) {
2578                              valueL *= 10.0;
2579                              while (valueL <= tolerance)
2580                                   valueL *= 10.0;
2581                         } else if (valueL > 1.0) {
2582                              valueL *= 0.1;
2583                              while (valueL > 1.0)
2584                                   valueL *= 0.1;
2585                         }
2586                         if (lowerValue > -1.0e20 && lowerValue)
2587                              lowerValue -= valueL;
2588                         double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2589                         // get in range
2590                         if (valueU <= tolerance) {
2591                              valueU *= 10.0;
2592                              while (valueU <= tolerance)
2593                                   valueU *= 10.0;
2594                         } else if (valueU > 1.0) {
2595                              valueU *= 0.1;
2596                              while (valueU > 1.0)
2597                                   valueU *= 0.1;
2598                         }
2599                         if (upperValue < 1.0e20 && upperValue)
2600                              upperValue += valueU;
2601                    }
2602               } else if (upperValue > 0.0) {
2603                    upperValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2604                    lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2605               } else if (upperValue < 0.0) {
2606                    upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2607                    lowerValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2608               } else {
2609               }
2610               if (lowerValue != lower_[i]) {
2611                    double difference = fabs(lowerValue - lower_[i]);
2612                    largest = CoinMax(largest, difference);
2613                    if (difference > fabs(lower_[i])*largestPerCent)
2614                         largestPerCent = fabs(difference / lower_[i]);
2615               }
2616               if (upperValue != upper_[i]) {
2617                    double difference = fabs(upperValue - upper_[i]);
2618                    largest = CoinMax(largest, difference);
2619                    if (difference > fabs(upper_[i])*largestPerCent)
2620                         largestPerCent = fabs(difference / upper_[i]);
2621               }
2622               if (printOut)
2623                    printf("row %d lower from %g to %g, upper from %g to %g\n",
2624                           i - numberColumns_, lower_[i], lowerValue, upper_[i], upperValue);
2625               lower_[i] = lowerValue;
2626               upper_[i] = upperValue;
2627          }
2628     }
2629     // Clean up
2630     for (i = 0; i < numberColumns_ + numberRows_; i++) {
2631          switch(getStatus(i)) {
2632
2633          case basic:
2634               break;
2635          case atUpperBound:
2636               solution_[i] = upper_[i];
2637               break;
2638          case isFixed:
2639          case atLowerBound:
2640               solution_[i] = lower_[i];
2641               break;
2642          case isFree:
2643               break;
2644          case superBasic:
2645               break;
2646          }
2647     }
2648     handler_->message(CLP_SIMPLEX_PERTURB, messages_)
2649               << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
2650               << CoinMessageEol;
2651     // redo nonlinear costs
2652     // say perturbed
2653     perturbation_ = 101;
2654}
2655// un perturb
2656bool
2657ClpSimplexPrimal::unPerturb()
2658{
2659     if (perturbation_ != 101)
2660          return false;
2661     // put back original bounds and costs
2662     createRim(1 + 4);
2663     sanityCheck();
2664     // unflag
2665     unflag();
2666     // get a valid nonlinear cost function
2667     delete nonLinearCost_;
2668     nonLinearCost_ = new ClpNonLinearCost(this);
2669     perturbation_ = 102; // stop any further perturbation
2670     // move non basic variables to new bounds
2671     nonLinearCost_->checkInfeasibilities(0.0);
2672#if 1
2673     // Try using dual
2674     return true;
2675#else
2676     gutsOfSolution(NULL, NULL, ifValuesPass != 0);
2677     return false;
2678#endif
2679
2680}
2681// Unflag all variables and return number unflagged
2682int
2683ClpSimplexPrimal::unflag()
2684{
2685     int i;
2686     int number = numberRows_ + numberColumns_;
2687     int numberFlagged = 0;
2688     // we can't really trust infeasibilities if there is dual error
2689     // allow tolerance bigger than standard to check on duals
2690     double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
2691     for (i = 0; i < number; i++) {
2692          if (flagged(i)) {
2693               clearFlagged(i);
2694               // only say if reasonable dj
2695               if (fabs(dj_[i]) > relaxedToleranceD)
2696                    numberFlagged++;
2697          }
2698     }
2699     numberFlagged += matrix_->generalExpanded(this, 8, i);
2700     if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
2701          printf("%d unflagged\n", numberFlagged);
2702     return numberFlagged;
2703}
2704// Do not change infeasibility cost and always say optimal
2705void
2706ClpSimplexPrimal::alwaysOptimal(bool onOff)
2707{
2708     if (onOff)
2709          specialOptions_ |= 1;
2710     else
2711          specialOptions_ &= ~1;
2712}
2713bool
2714ClpSimplexPrimal::alwaysOptimal() const
2715{
2716     return (specialOptions_ & 1) != 0;
2717}
2718// Flatten outgoing variables i.e. - always to exact bound
2719void
2720ClpSimplexPrimal::exactOutgoing(bool onOff)
2721{
2722     if (onOff)
2723          specialOptions_ |= 4;
2724     else
2725          specialOptions_ &= ~4;
2726}
2727bool
2728ClpSimplexPrimal::exactOutgoing() const
2729{
2730     return (specialOptions_ & 4) != 0;
2731}
2732/*
2733  Reasons to come out (normal mode/user mode):
2734  -1 normal
2735  -2 factorize now - good iteration/ NA
2736  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
2737  -4 inaccuracy - refactorize - no iteration/ NA
2738  -5 something flagged - go round again/ pivot not possible
2739  +2 looks unbounded
2740  +3 max iterations (iteration done)
2741*/
2742int
2743ClpSimplexPrimal::pivotResult(int ifValuesPass)
2744{
2745
2746     bool roundAgain = true;
2747     int returnCode = -1;
2748
2749     // loop round if user setting and doing refactorization
2750     while (roundAgain) {
2751          roundAgain = false;
2752          returnCode = -1;
2753          pivotRow_ = -1;
2754          sequenceOut_ = -1;
2755          rowArray_[1]->clear();
2756#if 0
2757          {
2758               int seq[] = {612, 643};
2759               int k;
2760               for (k = 0; k < sizeof(seq) / sizeof(int); k++) {
2761                    int iSeq = seq[k];
2762                    if (getColumnStatus(iSeq) != basic) {
2763                         double djval;
2764                         double * work;
2765                         int number;
2766                         int * which;
2767
2768                         int iIndex;
2769                         unpack(rowArray_[1], iSeq);
2770                         factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2771                         djval = cost_[iSeq];
2772                         work = rowArray_[1]->denseVector();
2773                         number = rowArray_[1]->getNumElements();
2774                         which = rowArray_[1]->getIndices();
2775
2776                         for (iIndex = 0; iIndex < number; iIndex++) {
2777
2778                              int iRow = which[iIndex];
2779                              double alpha = work[iRow];
2780                              int iPivot = pivotVariable_[iRow];
2781                              djval -= alpha * cost_[iPivot];
2782                         }
2783                         double comp = 1.0e-8 + 1.0e-7 * (CoinMax(fabs(dj_[iSeq]), fabs(djval)));
2784                         if (fabs(djval - dj_[iSeq]) > comp)
2785                              printf("Bad dj %g for %d - true is %g\n",
2786                                     dj_[iSeq], iSeq, djval);
2787                         assert (fabs(djval) < 1.0e-3 || djval * dj_[iSeq] > 0.0);
2788                         rowArray_[1]->clear();
2789                    }
2790               }
2791          }
2792#endif
2793
2794          // we found a pivot column
2795          // update the incoming column
2796          unpackPacked(rowArray_[1]);
2797          // save reduced cost
2798          double saveDj = dualIn_;
2799          factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
2800          // Get extra rows
2801          matrix_->extendUpdated(this, rowArray_[1], 0);
2802          // do ratio test and re-compute dj
2803#ifdef CLP_USER_DRIVEN
2804          if (solveType_ != 2 || (moreSpecialOptions_ & 512) == 0) {
2805#endif
2806               primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
2807                         ifValuesPass);
2808#ifdef CLP_USER_DRIVEN
2809          } else {
2810               int status = eventHandler_->event(ClpEventHandler::pivotRow);
2811               if (status >= 0) {
2812                    problemStatus_ = 5;
2813                    secondaryStatus_ = ClpEventHandler::pivotRow;
2814                    break;
2815               }
2816          }
2817#endif
2818          if (ifValuesPass) {
2819               saveDj = dualIn_;
2820               //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
2821               if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
2822                    if(fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
2823                         // try other way
2824                         directionIn_ = -directionIn_;
2825                         primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
2826                                   0);
2827                    }
2828                    if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
2829                         if (solveType_ == 1) {
2830                              // reject it
2831                              char x = isColumn(sequenceIn_) ? 'C' : 'R';
2832                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
2833                                        << x << sequenceWithin(sequenceIn_)
2834                                        << CoinMessageEol;
2835                              setFlagged(sequenceIn_);
2836                              progress_.clearBadTimes();
2837                              lastBadIteration_ = numberIterations_; // say be more cautious
2838                              clearAll();
2839                              pivotRow_ = -1;
2840                         }
2841                         returnCode = -5;
2842                         break;
2843                    }
2844               }
2845          }
2846          // need to clear toIndex_ in gub
2847          // ? when can I clear stuff
2848          // Clean up any gub stuff
2849          matrix_->extendUpdated(this, rowArray_[1], 1);
2850          double checkValue = 1.0e-2;
2851          if (largestDualError_ > 1.0e-5)
2852               checkValue = 1.0e-1;
2853          double test2 = dualTolerance_;
2854          double test1 = 1.0e-20;
2855#if 0 //def FEB_TRY
2856          if (factorization_->pivots() < 1) {
2857               test1 = -1.0e-4;
2858               if ((saveDj < 0.0 && dualIn_ < -1.0e-5 * dualTolerance_) ||
2859                         (saveDj > 0.0 && dualIn_ > 1.0e-5 * dualTolerance_))
2860                    test2 = 0.0; // allow through
2861          }
2862#endif
2863          if (!ifValuesPass && solveType_ == 1 && (saveDj * dualIn_ < test1 ||
2864                    fabs(saveDj - dualIn_) > checkValue*(1.0 + fabs(saveDj)) ||
2865                    fabs(dualIn_) < test2)) {
2866               if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) >
2867                         1.0e5)) {
2868                    char x = isColumn(sequenceIn_) ? 'C' : 'R';
2869                    handler_->message(CLP_PRIMAL_DJ, messages_)
2870                              << x << sequenceWithin(sequenceIn_)
2871                              << saveDj << dualIn_
2872                              << CoinMessageEol;
2873                    if(lastGoodIteration_ != numberIterations_) {
2874                         clearAll();
2875                         pivotRow_ = -1; // say no weights update
2876                         returnCode = -4;
2877                         if(lastGoodIteration_ + 1 == numberIterations_) {
2878                              // not looking wonderful - try cleaning bounds
2879                              // put non-basics to bounds in case tolerance moved
2880                              nonLinearCost_->checkInfeasibilities(0.0);
2881                         }
2882                         sequenceOut_ = -1;
2883                         break;
2884                    } else {
2885                         // take on more relaxed criterion
2886                         if (saveDj * dualIn_ < test1 ||
2887                                   fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) ||
2888                                   fabs(dualIn_) < test2) {
2889                              // need to reject something
2890                              char x = isColumn(sequenceIn_) ? 'C' : 'R';
2891                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
2892                                        << x << sequenceWithin(sequenceIn_)
2893                                        << CoinMessageEol;
2894                              setFlagged(sequenceIn_);
2895#ifdef FEB_TRY
2896                              // Make safer?
2897                              factorization_->saferTolerances (1.0e-15, -1.03);
2898#endif
2899                              progress_.clearBadTimes();
2900                              lastBadIteration_ = numberIterations_; // say be more cautious
2901                              clearAll();
2902                              pivotRow_ = -1;
2903                              returnCode = -5;
2904                              sequenceOut_ = -1;
2905                              break;
2906                         }
2907                    }
2908               } else {
2909                    //printf("%d %g %g\n",numberIterations_,saveDj,dualIn_);
2910               }
2911          }
2912          if (pivotRow_ >= 0) {
2913 #ifdef CLP_USER_DRIVEN1
2914               // Got good pivot - may need to unflag stuff
2915               userChoiceWasGood(this);
2916 #endif
2917               if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
2918                    // **** Coding for user interface
2919                    // do ray
2920                    primalRay(rowArray_[1]);
2921                    // update duals
2922                    // as packed need to find pivot row
2923                    //assert (rowArray_[1]->packedMode());
2924                    //int i;
2925
2926                    //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2927                    CoinAssert (fabs(alpha_) > 1.0e-8);
2928                    double multiplier = dualIn_ / alpha_;
2929                    rowArray_[0]->insert(pivotRow_, multiplier);
2930                    factorization_->updateColumnTranspose(rowArray_[2], rowArray_[0]);
2931                    // put row of tableau in rowArray[0] and columnArray[0]
2932                    matrix_->transposeTimes(this, -1.0,
2933                                            rowArray_[0], columnArray_[1], columnArray_[0]);
2934                    // update column djs
2935                    int i;
2936                    int * index = columnArray_[0]->getIndices();
2937                    int number = columnArray_[0]->getNumElements();
2938                    double * element = columnArray_[0]->denseVector();
2939                    for (i = 0; i < number; i++) {
2940                         int ii = index[i];
2941                         dj_[ii] += element[ii];
2942                         reducedCost_[ii] = dj_[ii];
2943                         element[ii] = 0.0;
2944                    }
2945                    columnArray_[0]->setNumElements(0);
2946                    // and row djs
2947                    index = rowArray_[0]->getIndices();
2948                    number = rowArray_[0]->getNumElements();
2949                    element = rowArray_[0]->denseVector();
2950                    for (i = 0; i < number; i++) {
2951                         int ii = index[i];
2952                         dj_[ii+numberColumns_] += element[ii];
2953                         dual_[ii] = dj_[ii+numberColumns_];
2954                         element[ii] = 0.0;
2955                    }
2956                    rowArray_[0]->setNumElements(0);
2957                    // check incoming
2958                    CoinAssert (fabs(dj_[sequenceIn_]) < 1.0e-1);
2959               }
2960               // if stable replace in basis
2961               // If gub or odd then alpha and pivotRow may change
2962               int updateType = 0;
2963               int updateStatus = matrix_->generalExpanded(this, 3, updateType);
2964               if (updateType >= 0)
2965                    updateStatus = factorization_->replaceColumn(this,
2966                                   rowArray_[2],
2967                                   rowArray_[1],
2968                                   pivotRow_,
2969                                   alpha_,
2970                                   (moreSpecialOptions_ & 16) != 0);
2971
2972               // if no pivots, bad update but reasonable alpha - take and invert
2973               if (updateStatus == 2 &&
2974                         lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
2975                    updateStatus = 4;
2976               if (updateStatus == 1 || updateStatus == 4) {
2977                    // slight error
2978                    if (factorization_->pivots() > 5 || updateStatus == 4) {
2979                         returnCode = -3;
2980                    }
2981               } else if (updateStatus == 2) {
2982                    // major error
2983                    // better to have small tolerance even if slower
2984                    factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
2985                    int maxFactor = factorization_->maximumPivots();
2986                    if (maxFactor > 10) {
2987                         if (forceFactorization_ < 0)
2988                              forceFactorization_ = maxFactor;
2989                         forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
2990                    }
2991                    // later we may need to unwind more e.g. fake bounds
2992                    if(lastGoodIteration_ != numberIterations_) {
2993                         clearAll();
2994                         pivotRow_ = -1;
2995                         if (solveType_ == 1 || (moreSpecialOptions_ & 512) != 0) {
2996                              returnCode = -4;
2997                              break;
2998                         } else {
2999                              // user in charge - re-factorize
3000                              int lastCleaned = 0;
3001                              ClpSimplexProgress dummyProgress;
3002                              if (saveStatus_)
3003                                   statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3004                              else
3005                                   statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3006                              roundAgain = true;
3007                              continue;
3008                         }
3009                    } else {
3010                         // need to reject something
3011                         if (solveType_ == 1) {
3012                              char x = isColumn(sequenceIn_) ? 'C' : 'R';
3013                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3014                                        << x << sequenceWithin(sequenceIn_)
3015                                        << CoinMessageEol;
3016                              setFlagged(sequenceIn_);
3017                              progress_.clearBadTimes();
3018                         }
3019                         lastBadIteration_ = numberIterations_; // say be more cautious
3020                         clearAll();
3021                         pivotRow_ = -1;
3022                         sequenceOut_ = -1;
3023                         returnCode = -5;
3024                         break;
3025
3026                    }
3027               } else if (updateStatus == 3) {
3028                    // out of memory
3029                    // increase space if not many iterations
3030                    if (factorization_->pivots() <
3031                              0.5 * factorization_->maximumPivots() &&
3032                              factorization_->pivots() < 200)
3033                         factorization_->areaFactor(
3034                              factorization_->areaFactor() * 1.1);
3035                    returnCode = -2; // factorize now
3036               } else if (updateStatus == 5) {
3037                    problemStatus_ = -2; // factorize now
3038               }
3039               // here do part of steepest - ready for next iteration
3040               if (!ifValuesPass)
3041                    primalColumnPivot_->updateWeights(rowArray_[1]);
3042          } else {
3043               if (pivotRow_ == -1) {
3044                    // no outgoing row is valid
3045                    if (valueOut_ != COIN_DBL_MAX) {
3046                         double objectiveChange = 0.0;
3047                         theta_ = valueOut_ - valueIn_;
3048                         updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3049                         solution_[sequenceIn_] += theta_;
3050                    }
3051                    rowArray_[0]->clear();
3052 #ifdef CLP_USER_DRIVEN1
3053                    /* Note if valueOut_ < COIN_DBL_MAX and
3054                       theta_ reasonable then this may be a valid sub flip */
3055                    if(!userChoiceValid2(this)) {
3056                      if (factorization_->pivots()<5) {
3057                        // flag variable
3058                        char x = isColumn(sequenceIn_) ? 'C' : 'R';
3059                        handler_->message(CLP_SIMPLEX_FLAG, messages_)
3060                          << x << sequenceWithin(sequenceIn_)
3061                          << CoinMessageEol;
3062                        setFlagged(sequenceIn_);
3063                        progress_.clearBadTimes();
3064                        roundAgain = true;
3065                        continue;
3066                      } else {
3067                        // try refactorizing first
3068                        returnCode = 4; //say looks odd but has iterated
3069                        break;
3070                      }
3071                    }
3072 #endif
3073                    if (!factorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3074                         returnCode = 2; //say looks unbounded
3075                         // do ray
3076                         primalRay(rowArray_[1]);
3077                    } else if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
3078                         // refactorize
3079                         int lastCleaned = 0;
3080                         ClpSimplexProgress dummyProgress;
3081                         if (saveStatus_)
3082                              statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3083                         else
3084                              statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3085                         roundAgain = true;
3086                         continue;
3087                    } else {
3088                         acceptablePivot_ = 1.0e-8;
3089                         returnCode = 4; //say looks unbounded but has iterated
3090                    }
3091                    break;
3092               } else {
3093                    // flipping from bound to bound
3094               }
3095          }
3096
3097          double oldCost = 0.0;
3098          if (sequenceOut_ >= 0)
3099               oldCost = cost_[sequenceOut_];
3100          // update primal solution
3101
3102          double objectiveChange = 0.0;
3103          // after this rowArray_[1] is not empty - used to update djs
3104          // If pivot row >= numberRows then may be gub
3105          int savePivot = pivotRow_;
3106          if (pivotRow_ >= numberRows_)
3107               pivotRow_ = -1;
3108          updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3109          pivotRow_ = savePivot;
3110
3111          double oldValue = valueIn_;
3112          if (directionIn_ == -1) {
3113               // as if from upper bound
3114               if (sequenceIn_ != sequenceOut_) {
3115                    // variable becoming basic
3116                    valueIn_ -= fabs(theta_);
3117               } else {
3118                    valueIn_ = lowerIn_;
3119               }
3120          } else {
3121               // as if from lower bound
3122               if (sequenceIn_ != sequenceOut_) {
3123                    // variable becoming basic
3124                    valueIn_ += fabs(theta_);
3125               } else {
3126                    valueIn_ = upperIn_;
3127               }
3128          }
3129          objectiveChange += dualIn_ * (valueIn_ - oldValue);
3130          // outgoing
3131          if (sequenceIn_ != sequenceOut_) {
3132               if (directionOut_ > 0) {
3133                    valueOut_ = lowerOut_;
3134               } else {
3135                    valueOut_ = upperOut_;
3136               }
3137               if(valueOut_ < lower_[sequenceOut_] - primalTolerance_)
3138                    valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
3139               else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
3140                    valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
3141               // may not be exactly at bound and bounds may have changed
3142               // Make sure outgoing looks feasible
3143               directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
3144               // May have got inaccurate
3145               //if (oldCost!=cost_[sequenceOut_])
3146               //printf("costchange on %d from %g to %g\n",sequenceOut_,
3147               //       oldCost,cost_[sequenceOut_]);
3148               if (solveType_ != 2)
3149                    dj_[sequenceOut_] = cost_[sequenceOut_] - oldCost; // normally updated next iteration
3150               solution_[sequenceOut_] = valueOut_;
3151          }
3152          // change cost and bounds on incoming if primal
3153          nonLinearCost_->setOne(sequenceIn_, valueIn_);
3154          int whatNext = housekeeping(objectiveChange);
3155          //nonLinearCost_->validate();
3156#if CLP_DEBUG >1
3157          {
3158               double sum;
3159               int ninf = matrix_->checkFeasible(this, sum);
3160               if (ninf)
3161                    printf("infeas %d\n", ninf);
3162          }
3163#endif
3164          if (whatNext == 1) {
3165               returnCode = -2; // refactorize
3166          } else if (whatNext == 2) {
3167               // maximum iterations or equivalent
3168               returnCode = 3;
3169          } else if(numberIterations_ == lastGoodIteration_
3170                    + 2 * factorization_->maximumPivots()) {
3171               // done a lot of flips - be safe
3172               returnCode = -2; // refactorize
3173          }
3174          // Check event
3175          {
3176               int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3177               if (status >= 0) {
3178                    problemStatus_ = 5;
3179                    secondaryStatus_ = ClpEventHandler::endOfIteration;
3180                    returnCode = 3;
3181               }
3182          }
3183     }
3184     if ((solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) &&
3185               (returnCode == -2 || returnCode == -3)) {
3186          // refactorize here
3187          int lastCleaned = 0;
3188          ClpSimplexProgress dummyProgress;
3189          if (saveStatus_)
3190               statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3191          else
3192               statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3193          if (problemStatus_ == 5) {
3194               printf("Singular basis\n");
3195               problemStatus_ = -1;
3196               returnCode = 5;
3197          }
3198     }
3199#ifdef CLP_DEBUG
3200     {
3201          int i;
3202          // not [1] as may have information
3203          for (i = 0; i < 4; i++) {
3204               if (i != 1)
3205                    rowArray_[i]->checkClear();
3206          }
3207          for (i = 0; i < 2; i++) {
3208               columnArray_[i]->checkClear();
3209          }
3210     }
3211#endif
3212     return returnCode;
3213}
3214// Create primal ray
3215void
3216ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
3217{
3218     delete [] ray_;
3219     ray_ = new double [numberColumns_];
3220     CoinZeroN(ray_, numberColumns_);
3221     int number = rowArray->getNumElements();
3222     int * index = rowArray->getIndices();
3223     double * array = rowArray->denseVector();
3224     double way = -directionIn_;
3225     int i;
3226     double zeroTolerance = 1.0e-12;
3227     if (sequenceIn_ < numberColumns_)
3228          ray_[sequenceIn_] = directionIn_;
3229     if (!rowArray->packedMode()) {
3230          for (i = 0; i < number; i++) {
3231               int iRow = index[i];
3232               int iPivot = pivotVariable_[iRow];
3233               double arrayValue = array[iRow];
3234               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3235                    ray_[iPivot] = way * arrayValue;
3236          }
3237     } else {
3238          for (i = 0; i < number; i++) {
3239               int iRow = index[i];
3240               int iPivot = pivotVariable_[iRow];
3241               double arrayValue = array[i];
3242               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3243                    ray_[iPivot] = way * arrayValue;
3244          }
3245     }
3246}
3247/* Get next superbasic -1 if none,
3248   Normal type is 1
3249   If type is 3 then initializes sorted list
3250   if 2 uses list.
3251*/
3252int
3253ClpSimplexPrimal::nextSuperBasic(int superBasicType,
3254                                 CoinIndexedVector * columnArray)
3255{
3256     int returnValue = -1;
3257     bool finished = false;
3258     while (!finished) {
3259          returnValue = firstFree_;
3260          int iColumn = firstFree_ + 1;
3261          if (superBasicType > 1) {
3262               if (superBasicType > 2) {
3263                    // Initialize list
3264                    // Wild guess that lower bound more natural than upper
3265                    int number = 0;
3266                    double * work = columnArray->denseVector();
3267                    int * which = columnArray->getIndices();
3268                    for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
3269                         if (!flagged(iColumn)) {
3270                              if (getStatus(iColumn) == superBasic) {
3271                                   if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3272                                        solution_[iColumn] = lower_[iColumn];
3273                                        setStatus(iColumn, atLowerBound);
3274                                   } else if (fabs(solution_[iColumn] - upper_[iColumn])
3275                                              <= primalTolerance_) {
3276                                        solution_[iColumn] = upper_[iColumn];
3277                                        setStatus(iColumn, atUpperBound);
3278                                   } else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3279                                        setStatus(iColumn, isFree);
3280                                        break;
3281                                   } else if (!flagged(iColumn)) {
3282                                        // put ones near bounds at end after sorting
3283                                        work[number] = - CoinMin(0.1 * (solution_[iColumn] - lower_[iColumn]),
3284                                                                 upper_[iColumn] - solution_[iColumn]);
3285                                        which[number++] = iColumn;
3286                                   }
3287                              }
3288                         }
3289                    }
3290                    CoinSort_2(work, work + number, which);
3291                    columnArray->setNumElements(number);
3292                    CoinZeroN(work, number);
3293               }
3294               int * which = columnArray->getIndices();
3295               int number = columnArray->getNumElements();
3296               if (!number) {
3297                    // finished
3298                    iColumn = numberRows_ + numberColumns_;
3299                    returnValue = -1;
3300               } else {
3301                    number--;
3302                    returnValue = which[number];
3303                    iColumn = returnValue;
3304                    columnArray->setNumElements(number);
3305               }
3306          } else {
3307               for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
3308                    if (!flagged(iColumn)) {
3309                         if (getStatus(iColumn) == superBasic) {
3310                              if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3311                                   solution_[iColumn] = lower_[iColumn];
3312                                   setStatus(iColumn, atLowerBound);
3313                              } else if (fabs(solution_[iColumn] - upper_[iColumn])
3314                                         <= primalTolerance_) {
3315                                   solution_[iColumn] = upper_[iColumn];
3316                                   setStatus(iColumn, atUpperBound);
3317                              } else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3318                                   setStatus(iColumn, isFree);
3319                                   break;
3320                              } else {
3321                                   break;
3322                              }
3323                         }
3324                    }
3325               }
3326          }
3327          firstFree_ = iColumn;
3328          finished = true;
3329          if (firstFree_ == numberRows_ + numberColumns_)
3330               firstFree_ = -1;
3331          if (returnValue >= 0 && getStatus(returnValue) != superBasic && getStatus(returnValue) != isFree)
3332               finished = false; // somehow picked up odd one
3333     }
3334     return returnValue;
3335}
3336void
3337ClpSimplexPrimal::clearAll()
3338{
3339     // Clean up any gub stuff
3340     matrix_->extendUpdated(this, rowArray_[1], 1);
3341     int number = rowArray_[1]->getNumElements();
3342     int * which = rowArray_[1]->getIndices();
3343
3344     int iIndex;
3345     for (iIndex = 0; iIndex < number; iIndex++) {
3346
3347          int iRow = which[iIndex];
3348          clearActive(iRow);
3349     }
3350     rowArray_[1]->clear();
3351     // make sure any gub sets are clean
3352     matrix_->generalExpanded(this, 11, sequenceIn_);
3353
3354}
3355// Sort of lexicographic resolve
3356int
3357ClpSimplexPrimal::lexSolve()
3358{
3359     algorithm_ = +1;
3360     //specialOptions_ |= 4;
3361
3362     // save data
3363     ClpDataSave data = saveData();
3364     matrix_->refresh(this); // make sure matrix okay
3365
3366     // Save so can see if doing after dual
3367     int initialStatus = problemStatus_;
3368     int initialIterations = numberIterations_;
3369     int initialNegDjs = -1;
3370     // initialize - maybe values pass and algorithm_ is +1
3371     int ifValuesPass = 0;
3372#if 0
3373     // if so - put in any superbasic costed slacks
3374     // Start can skip some things in transposeTimes
3375     specialOptions_ |= 131072;
3376     if (ifValuesPass && specialOptions_ < 0x01000000) {
3377          // Get column copy
3378          const CoinPackedMatrix * columnCopy = matrix();
3379          const int * row = columnCopy->getIndices();
3380          const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3381          const int * columnLength = columnCopy->getVectorLengths();
3382          //const double * element = columnCopy->getElements();
3383          int n = 0;
3384          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3385               if (columnLength[iColumn] == 1) {
3386                    Status status = getColumnStatus(iColumn);
3387                    if (status != basic && status != isFree) {
3388                         double value = columnActivity_[iColumn];
3389                         if (fabs(value - columnLower_[iColumn]) > primalTolerance_ &&
3390                                   fabs(value - columnUpper_[iColumn]) > primalTolerance_) {
3391                              int iRow = row[columnStart[iColumn]];
3392                              if (getRowStatus(iRow) == basic) {
3393                                   setRowStatus(iRow, superBasic);
3394                                   setColumnStatus(iColumn, basic);
3395                                   n++;
3396                              }
3397                         }
3398                    }
3399               }
3400          }
3401          printf("%d costed slacks put in basis\n", n);
3402     }
3403#endif
3404     double * originalCost = NULL;
3405     double * originalLower = NULL;
3406     double * originalUpper = NULL;
3407     if (!startup(0, 0)) {
3408
3409          // Set average theta
3410          nonLinearCost_->setAverageTheta(1.0e3);
3411          int lastCleaned = 0; // last time objective or bounds cleaned up
3412
3413          // Say no pivot has occurred (for steepest edge and updates)
3414          pivotRow_ = -2;
3415
3416          // This says whether to restore things etc
3417          int factorType = 0;
3418          if (problemStatus_ < 0 && perturbation_ < 100) {
3419               perturb(0);
3420               // Can't get here if values pass
3421               assert (!ifValuesPass);
3422               gutsOfSolution(NULL, NULL);
3423               if (handler_->logLevel() > 2) {
3424                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
3425                              << numberIterations_ << objectiveValue();
3426                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3427                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3428                    handler_->printing(sumDualInfeasibilities_ > 0.0)
3429                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3430                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
3431                                       < numberDualInfeasibilities_)
3432                              << numberDualInfeasibilitiesWithoutFree_;
3433                    handler_->message() << CoinMessageEol;
3434               }
3435          }
3436          ClpSimplex * saveModel = NULL;
3437          int stopSprint = -1;
3438          int sprintPass = 0;
3439          int reasonableSprintIteration = 0;
3440          int lastSprintIteration = 0;
3441          double lastObjectiveValue = COIN_DBL_MAX;
3442          // Start check for cycles
3443          progress_.fillFromModel(this);
3444          progress_.startCheck();
3445          /*
3446            Status of problem:
3447            0 - optimal
3448            1 - infeasible
3449            2 - unbounded
3450            -1 - iterating
3451            -2 - factorization wanted
3452            -3 - redo checking without factorization
3453            -4 - looks infeasible
3454            -5 - looks unbounded
3455          */
3456          originalCost = CoinCopyOfArray(cost_, numberColumns_ + numberRows_);
3457          originalLower = CoinCopyOfArray(lower_, numberColumns_ + numberRows_);
3458          originalUpper = CoinCopyOfArray(upper_, numberColumns_ + numberRows_);
3459          while (problemStatus_ < 0) {
3460               int iRow, iColumn;
3461               // clear
3462               for (iRow = 0; iRow < 4; iRow++) {
3463                    rowArray_[iRow]->clear();
3464               }
3465
3466               for (iColumn = 0; iColumn < 2; iColumn++) {
3467                    columnArray_[iColumn]->clear();
3468               }
3469
3470               // give matrix (and model costs and bounds a chance to be
3471               // refreshed (normally null)
3472               matrix_->refresh(this);
3473               // If getting nowhere - why not give it a kick
3474#if 1
3475               if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
3476                         && initialStatus != 10) {
3477                    perturb(1);
3478                    matrix_->rhsOffset(this, true, false);
3479               }
3480#endif
3481               // If we have done no iterations - special
3482               if (lastGoodIteration_ == numberIterations_ && factorType)
3483                    factorType = 3;
3484               if (saveModel) {
3485                    // Doing sprint
3486                    if (sequenceIn_ < 0 || numberIterations_ >= stopSprint) {
3487                         problemStatus_ = -1;
3488                         originalModel(saveModel);
3489                         saveModel = NULL;
3490                         if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration &&
3491                                   sprintPass > 100)
3492                              primalColumnPivot_->switchOffSprint();
3493                         //lastSprintIteration=numberIterations_;
3494                         printf("End small model\n");
3495                    }
3496               }
3497
3498               // may factorize, checks if problem finished
3499               statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3500               if (initialStatus == 10) {
3501                    // cleanup phase
3502                    if(initialIterations != numberIterations_) {
3503                         if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
3504                              // getting worse - try perturbing
3505                              if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
3506                                   perturb(1);
3507                                   matrix_->rhsOffset(this, true, false);
3508                                   statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3509                              }
3510                         }
3511                    } else {
3512                         // save number of negative djs
3513                         if (!numberPrimalInfeasibilities_)
3514                              initialNegDjs = numberDualInfeasibilities_;
3515                         // make sure weight won't be changed
3516                         if (infeasibilityCost_ == 1.0e10)
3517                              infeasibilityCost_ = 1.000001e10;
3518                    }
3519               }
3520               // See if sprint says redo because of problems
3521               if (numberDualInfeasibilities_ == -776) {
3522                    // Need new set of variables
3523                    problemStatus_ = -1;
3524                    originalModel(saveModel);
3525                    saveModel = NULL;
3526                    //lastSprintIteration=numberIterations_;
3527                    printf("End small model after\n");
3528                    statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3529               }
3530               int numberSprintIterations = 0;
3531               int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
3532               if (problemStatus_ == 777) {
3533                    // problems so do one pass with normal
3534                    problemStatus_ = -1;
3535                    originalModel(saveModel);
3536                    saveModel = NULL;
3537                    // Skip factorization
3538                    //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3539                    statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3540               } else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
3541                    int numberSort = 0;
3542                    int numberFixed = 0;
3543                    int numberBasic = 0;
3544                    reasonableSprintIteration = numberIterations_ + 100;
3545                    int * whichColumns = new int[numberColumns_];
3546                    double * weight = new double[numberColumns_];
3547                    int numberNegative = 0;
3548                    double sumNegative = 0.0;
3549                    // now massage weight so all basic in plus good djs
3550                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3551                         double dj = dj_[iColumn];
3552                         switch(getColumnStatus(iColumn)) {
3553
3554                         case basic:
3555                              dj = -1.0e50;
3556                              numberBasic++;
3557                              break;
3558                         case atUpperBound:
3559                              dj = -dj;
3560                              break;
3561                         case isFixed:
3562                              dj = 1.0e50;
3563                              numberFixed++;
3564                              break;
3565                         case atLowerBound:
3566                              dj = dj;
3567                              break;
3568                         case isFree:
3569                              dj = -100.0 * fabs(dj);
3570                              break;
3571                         case superBasic:
3572                              dj = -100.0 * fabs(dj);
3573                              break;
3574                         }
3575                         if (dj < -dualTolerance_ && dj > -1.0e50) {
3576                              numberNegative++;
3577                              sumNegative -= dj;
3578                         }
3579                         weight[iColumn] = dj;
3580                         whichColumns[iColumn] = iColumn;
3581                    }
3582                    handler_->message(CLP_SPRINT, messages_)
3583                              << sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
3584                              << numberNegative
3585                              << CoinMessageEol;
3586                    sprintPass++;
3587                    lastSprintIteration = numberIterations_;
3588                    if (objectiveValue()*optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
3589                         // switch off
3590                         printf("Switching off sprint\n");
3591                         primalColumnPivot_->switchOffSprint();
3592                    } else {
3593                         lastObjectiveValue = objectiveValue() * optimizationDirection_;
3594                         // sort
3595                         CoinSort_2(weight, weight + numberColumns_, whichColumns);
3596                         numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
3597                         // Sort to make consistent ?
3598                         std::sort(whichColumns, whichColumns + numberSort);
3599                         saveModel = new ClpSimplex(this, numberSort, whichColumns);
3600                         delete [] whichColumns;
3601                         delete [] weight;
3602                         // Skip factorization
3603                         //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3604                         //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
3605                         stopSprint = numberIterations_ + numberSprintIterations;
3606                         printf("Sprint with %d columns for %d iterations\n",
3607                                numberSprintColumns, numberSprintIterations);
3608                    }
3609               }
3610
3611               // Say good factorization
3612               factorType = 1;
3613
3614               // Say no pivot has occurred (for steepest edge and updates)
3615               pivotRow_ = -2;
3616
3617               // exit if victory declared
3618               if (problemStatus_ >= 0) {
3619                    if (originalCost) {
3620                         // find number nonbasic with zero reduced costs
3621                         int numberDegen = 0;
3622                         int numberTotal = numberColumns_; //+numberRows_;
3623                         for (int i = 0; i < numberTotal; i++) {
3624                              cost_[i] = 0.0;
3625                              if (getStatus(i) == atLowerBound) {
3626                                   if (dj_[i] <= dualTolerance_) {
3627                                        cost_[i] = numberTotal - i + randomNumberGenerator_.randomDouble() * 0.5;
3628                                        numberDegen++;
3629                                   } else {
3630                                        // fix
3631                                        cost_[i] = 1.0e10; //upper_[i]=lower_[i];
3632                                   }
3633                              } else if (getStatus(i) == atUpperBound) {
3634                                   if (dj_[i] >= -dualTolerance_) {
3635                                        cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3636                                        numberDegen++;
3637                                   } else {
3638                                        // fix
3639                                        cost_[i] = -1.0e10; //lower_[i]=upper_[i];
3640                                   }
3641                              } else if (getStatus(i) == basic) {
3642                                   cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3643                              }
3644                         }
3645                         problemStatus_ = -1;
3646                         lastObjectiveValue = COIN_DBL_MAX;
3647                         // Start check for cycles
3648                         progress_.fillFromModel(this);
3649                         progress_.startCheck();
3650                         printf("%d degenerate after %d iterations\n", numberDegen,
3651                                numberIterations_);
3652                         if (!numberDegen) {
3653                              CoinMemcpyN(originalCost, numberTotal, cost_);
3654                              delete [] originalCost;
3655                              originalCost = NULL;
3656                              CoinMemcpyN(originalLower, numberTotal, lower_);
3657                              delete [] originalLower;
3658                              CoinMemcpyN(originalUpper, numberTotal, upper_);
3659                              delete [] originalUpper;
3660                         }
3661                         delete nonLinearCost_;
3662                         nonLinearCost_ = new ClpNonLinearCost(this);
3663                         progress_.endOddState();
3664                         continue;
3665                    } else {
3666                         printf("exiting after %d iterations\n", numberIterations_);
3667                         break;
3668                    }
3669               }
3670
3671               // test for maximum iterations
3672               if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
3673                    problemStatus_ = 3;
3674                    break;
3675               }
3676
3677               if (firstFree_ < 0) {
3678                    if (ifValuesPass) {
3679                         // end of values pass
3680                         ifValuesPass = 0;
3681                         int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
3682                         if (status >= 0) {
3683                              problemStatus_ = 5;
3684                              secondaryStatus_ = ClpEventHandler::endOfValuesPass;
3685                              break;
3686                         }
3687                    }
3688               }
3689               // Check event
3690               {
3691                    int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3692                    if (status >= 0) {
3693                         problemStatus_ = 5;
3694                         secondaryStatus_ = ClpEventHandler::endOfFactorization;
3695                         break;
3696                    }
3697               }
3698               // Iterate
3699               whileIterating(ifValuesPass ? 1 : 0);
3700          }
3701     }
3702     assert (!originalCost);
3703     // if infeasible get real values
3704     //printf("XXXXY final cost %g\n",infeasibilityCost_);
3705     progress_.initialWeight_ = 0.0;
3706     if (problemStatus_ == 1 && secondaryStatus_ != 6) {
3707          infeasibilityCost_ = 0.0;
3708          createRim(1 + 4);
3709          nonLinearCost_->checkInfeasibilities(0.0);
3710          sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
3711          numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
3712          // and get good feasible duals
3713          computeDuals(NULL);
3714     }
3715     // Stop can skip some things in transposeTimes
3716     specialOptions_ &= ~131072;
3717     // clean up
3718     unflag();
3719     finish(0);
3720     restoreData(data);
3721     return problemStatus_;
3722}
Note: See TracBrowser for help on using the repository browser.