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

Last change on this file since 2080 was 2080, checked in by forrest, 5 years ago

try and take out null message and extra solve

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