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

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

allow primal a bit more slack to reduce objective

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