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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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