source: stable/1.16/Clp/src/ClpSimplexPrimal.cpp @ 2163

Last change on this file since 2163 was 2163, checked in by forrest, 4 years ago

fix over-zealous saying infeasible

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