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

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

try and trap bad primal bases

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