source: stable/1.15/Clp/src/ClpSimplexPrimal.cpp @ 2006

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

changes for parallel and idiot

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