source: stable/1.15/Clp/src/ClpSimplexPrimal.cpp @ 1931

Last change on this file since 1931 was 1931, checked in by stefan, 7 years ago

sync with trunk rev 1930

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