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

Last change on this file since 2327 was 2327, checked in by forrest, 12 months ago

stupid mistake

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