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

Last change on this file since 2254 was 2254, checked in by forrest, 2 years ago

segfault in primal (and out some messages)

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