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

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

more flexibility in pivoting

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