# source:trunk/Clp/src/ClpSimplexDual.cpp@1710

Last change on this file since 1710 was 1710, checked in by forrest, 9 years ago

fix for initialSolve being called from odd place and fixes for difficult problems

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Id`
File size: 323.0 KB
Line
1/* \$Id: ClpSimplexDual.cpp 1710 2011-03-24 14:05:58Z forrest \$ */
5
6
7/* Notes on implementation of dual simplex algorithm.
8
9   When dual feasible:
10
11   If primal feasible, we are optimal.  Otherwise choose an infeasible
12   basic variable to leave basis (normally going to nearest bound) (B).  We
13   now need to find an incoming variable which will leave problem
14   dual feasible so we get the row of the tableau corresponding to
15   the basic variable (with the correct sign depending if basic variable
16   above or below feasibility region - as that affects whether reduced
17   cost on outgoing variable has to be positive or negative).
18
19   We now perform a ratio test to determine which incoming variable will
20   preserve dual feasibility (C).  If no variable found then problem
21   is infeasible (in primal sense).  If there is a variable, we then
22   perform pivot and repeat.  Trivial?
23
24   -------------------------------------------
25
26   A) How do we get dual feasible?  If all variables have bounds then
27   it is trivial to get feasible by putting non-basic variables to
28   correct bounds.  OSL did not have a phase 1/phase 2 approach but
29   instead effectively put fake bounds on variables and this is the
30   approach here, although I had hoped to make it cleaner.
31
32   If there is a weight of X on getting dual feasible:
33     Non-basic variables with negative reduced costs are put to
34     lesser of their upper bound and their lower bound + X.
35     Similarly, mutatis mutandis, for positive reduced costs.
36
37   Free variables should normally be in basis, otherwise I have
38   coding which may be able to come out (and may not be correct).
39
40   In OSL, this weight was changed heuristically, here at present
41   it is only increased if problem looks finished.  If problem is
42   feasible I check for unboundedness.  If not unbounded we
43   could play with going into primal.  As long as weights increase
44   any algorithm would be finite.
45
46   B) Which outgoing variable to choose is a virtual base class.
47   For difficult problems steepest edge is preferred while for
48   very easy (large) problems we will need partial scan.
49
50   C) Sounds easy, but this is hardest part of algorithm.
51      1) Instead of stopping at first choice, we may be able
52      to flip that variable to other bound and if objective
53      still improving choose again.  These mini iterations can
54      increase speed by orders of magnitude but we may need to
55      go to more of a bucket choice of variable rather than looking
56      at them one by one (for speed).
57      2) Accuracy.  Reduced costs may be of wrong sign but less than
58      tolerance.  Pivoting on these makes objective go backwards.
59      OSL modified cost so a zero move was made, Gill et al
60      (in primal analogue) modified so a strictly positive move was
61      made.  It is not quite as neat in dual but that is what we
62      try and do.  The two problems are that re-factorizations can
63      change reduced costs above and below tolerances and that when
64      finished we need to reset costs and try again.
65      3) Degeneracy.  Gill et al helps but may not be enough.  We
66      may need more.  Also it can improve speed a lot if we perturb
67      the costs significantly.
68
69  References:
70     Forrest and Goldfarb, Steepest-edge simplex algorithms for
71       linear programming - Mathematical Programming 1992
72     Forrest and Tomlin, Implementing the simplex method for
73       the Optimization Subroutine Library - IBM Systems Journal 1992
74     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
75       Procedure for Linear and Nonlinear Programming SOL report 1988
76
77
78  TODO:
79
80  a) Better recovery procedures.  At present I never check on forward
82     re-factorization frequency, but this is only on singular
83     factorizations.
84  b) Fast methods for large easy problems (and also the option for
85     the code to automatically choose which method).
86  c) We need to be able to stop in various ways for OSI - this
87     is fairly easy.
88
89 */
90#ifdef COIN_DEVELOP
91#undef COIN_DEVELOP
92#define COIN_DEVELOP 2
93#endif
94
95#include "CoinPragma.hpp"
96
97#include <math.h>
98
99#include "CoinHelperFunctions.hpp"
100#include "ClpHelperFunctions.hpp"
101#include "ClpSimplexDual.hpp"
102#include "ClpEventHandler.hpp"
103#include "ClpFactorization.hpp"
104#include "CoinPackedMatrix.hpp"
105#include "CoinIndexedVector.hpp"
106#include "CoinFloatEqual.hpp"
107#include "ClpDualRowDantzig.hpp"
108#include "ClpMessage.hpp"
109#include "ClpLinearObjective.hpp"
110#include <cfloat>
111#include <cassert>
112#include <string>
113#include <stdio.h>
114#include <iostream>
115//#define CLP_DEBUG 1
116// To force to follow another run put logfile name here and define
117//#define FORCE_FOLLOW
118#ifdef FORCE_FOLLOW
119static FILE * fpFollow = NULL;
120static char * forceFile = "old.log";
121static int force_in = -1;
122static int force_out = -1;
123static int force_iteration = 0;
124#endif
125//#define VUB
126#ifdef VUB
127extern int * vub;
128extern int * toVub;
129extern int * nextDescendent;
130#endif
131#ifdef NDEBUG
132#define NDEBUG_CLP
133#endif
134#ifndef CLP_INVESTIGATE
135#define NDEBUG_CLP
136#endif
137// dual
138
139/* *** Method
140   This is a vanilla version of dual simplex.
141
142   It tries to be a single phase approach with a weight of 1.0 being
143   given to getting optimal and a weight of dualBound_ being
144   given to getting dual feasible.  In this version I have used the
145   idea that this weight can be thought of as a fake bound.  If the
146   distance between the lower and upper bounds on a variable is less
147   than the feasibility weight then we are always better off flipping
148   to other bound to make dual feasible.  If the distance is greater
149   then we make up a fake bound dualBound_ away from one bound.
150   If we end up optimal or primal infeasible, we check to see if
151   bounds okay.  If so we have finished, if not we increase dualBound_
152   and continue (after checking if unbounded). I am undecided about
153   free variables - there is coding but I am not sure about it.  At
154   present I put them in basis anyway.
155
156   The code is designed to take advantage of sparsity so arrays are
157   seldom zeroed out from scratch or gone over in their entirety.
158   The only exception is a full scan to find outgoing variable.  This
159   will be changed to keep an updated list of infeasibilities (or squares
160   if steepest edge).  Also on easy problems we don't need full scan - just
161   pick first reasonable.
162
163   One problem is how to tackle degeneracy and accuracy.  At present
164   I am using the modification of costs which I put in OSL and which was
165   extended by Gill et al.  I am still not sure of the exact details.
166
167   The flow of dual is three while loops as follows:
168
169   while (not finished) {
170
171     while (not clean solution) {
172
173        Factorize and/or clean up solution by flipping variables so
174  dual feasible.  If looks finished check fake dual bounds.
175  Repeat until status is iterating (-1) or finished (0,1,2)
176
177     }
178
179     while (status==-1) {
180
181       Iterate until no pivot in or out or time to re-factorize.
182
183       Flow is:
184
185       choose pivot row (outgoing variable).  if none then
186 we are primal feasible so looks as if done but we need to
187 break and check bounds etc.
188
189 Get pivot row in tableau
190
191       Choose incoming column.  If we don't find one then we look
192 primal infeasible so break and check bounds etc.  (Also the
193 pivot tolerance is larger after any iterations so that may be
194 reason)
195
196       If we do find incoming column, we may have to adjust costs to
197 keep going forwards (anti-degeneracy).  Check pivot will be stable
198 and if unstable throw away iteration (we will need to implement
199 flagging of basic variables sometime) and break to re-factorize.
200 If minor error re-factorize after iteration.
201
202 Update everything (this may involve flipping variables to stay
203 dual feasible.
204
205     }
206
207   }
208
209   At present we never check we are going forwards.  I overdid that in
210   OSL so will try and make a last resort.
211
212   Needs partial scan pivot out option.
213   Needs dantzig, uninitialized and full steepest edge options (can still
214   use partial scan)
215
216   May need other anti-degeneracy measures, especially if we try and use
217   loose tolerances as a way to solve in fewer iterations.
218
219   I like idea of dynamic scaling.  This gives opportunity to decouple
220   different implications of scaling for accuracy, iteration count and
221   feasibility tolerance.
222
223*/
224#define CLEAN_FIXED 0
225// Startup part of dual (may be extended to other algorithms)
226int
227ClpSimplexDual::startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions)
228{
229     // If values pass then save given duals round check solution
230     // sanity check
231     // initialize - no values pass and algorithm_ is -1
232     // put in standard form (and make row copy)
233     // create modifiable copies of model rim and do optional scaling
234     // If problem looks okay
235     // Do initial factorization
236     // If user asked for perturbation - do it
237     numberFake_ = 0; // Number of variables at fake bounds
238     numberChanged_ = 0; // Number of variables with changed costs
239     if (!startup(0, startFinishOptions)) {
240          int usePrimal = 0;
241          // looks okay
242          // Superbasic variables not allowed
243          // If values pass then scale pi
244          if (ifValuesPass) {
245               if (problemStatus_ && perturbation_ < 100)
246                    usePrimal = perturb();
247               int i;
248               if (scalingFlag_ > 0) {
249                    for (i = 0; i < numberRows_; i++) {
250                         dual_[i] = saveDuals[i] * inverseRowScale_[i];
251                    }
252               } else {
253                    CoinMemcpyN(saveDuals, numberRows_, dual_);
254               }
255               // now create my duals
256               for (i = 0; i < numberRows_; i++) {
257                    // slack
258                    double value = dual_[i];
259                    value += rowObjectiveWork_[i];
260                    saveDuals[i+numberColumns_] = value;
261               }
262               CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
263               transposeTimes(-1.0, dual_, saveDuals);
264               // make reduced costs okay
265               for (i = 0; i < numberColumns_; i++) {
266                    if (getStatus(i) == atLowerBound) {
267                         if (saveDuals[i] < 0.0) {
268                              //if (saveDuals[i]<-1.0e-3)
269                              //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
270                              saveDuals[i] = 0.0;
271                         }
272                    } else if (getStatus(i) == atUpperBound) {
273                         if (saveDuals[i] > 0.0) {
274                              //if (saveDuals[i]>1.0e-3)
275                              //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
276                              saveDuals[i] = 0.0;
277                         }
278                    }
279               }
280               CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
281               // set up possible ones
282               for (i = 0; i < numberRows_ + numberColumns_; i++)
283                    clearPivoted(i);
284               int iRow;
285               for (iRow = 0; iRow < numberRows_; iRow++) {
286                    int iPivot = pivotVariable_[iRow];
287                    if (fabs(saveDuals[iPivot]) > dualTolerance_) {
288                         if (getStatus(iPivot) != isFree)
289                              setPivoted(iPivot);
290                    }
291               }
292          } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
293               // set up possible ones
294               for (int i = 0; i < numberRows_ + numberColumns_; i++)
295                    clearPivoted(i);
296               int iRow;
297               for (iRow = 0; iRow < numberRows_; iRow++) {
298                    int iPivot = pivotVariable_[iRow];
299                    if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
300                         setPivoted(iPivot);
301                    }
302               }
303          }
304
305          double objectiveChange;
306          assert (!numberFake_);
307          assert (numberChanged_ == 0);
308          if (!numberFake_) // if nonzero then adjust
309               changeBounds(1, NULL, objectiveChange);
310
311          if (!ifValuesPass) {
312               // Check optimal
313               if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
314                    problemStatus_ = 0;
315          }
316          if (problemStatus_ < 0 && perturbation_ < 100) {
317               bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
318               if (!inCbcOrOther)
319                    usePrimal = perturb();
320               // Can't get here if values pass
321               gutsOfSolution(NULL, NULL);
322#ifdef CLP_INVESTIGATE
323               if (numberDualInfeasibilities_)
324                    printf("ZZZ %d primal %d dual - sumdinf %g\n",
325                           numberPrimalInfeasibilities_,
326                           numberDualInfeasibilities_, sumDualInfeasibilities_);
327#endif
328               if (handler_->logLevel() > 2) {
329                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
330                              << numberIterations_ << objectiveValue();
331                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
332                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
333                    handler_->printing(sumDualInfeasibilities_ > 0.0)
334                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
335                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
336                                       < numberDualInfeasibilities_)
337                              << numberDualInfeasibilitiesWithoutFree_;
338                    handler_->message() << CoinMessageEol;
339               }
340               if (inCbcOrOther) {
341                    if (numberPrimalInfeasibilities_) {
342                         usePrimal = perturb();
343                         if (perturbation_ >= 101) {
344                              computeDuals(NULL);
345                              //gutsOfSolution(NULL,NULL);
346                              checkDualSolution(); // recompute objective
347                         }
348                    } else if (numberDualInfeasibilities_) {
349                         problemStatus_ = 10;
350                         if ((moreSpecialOptions_ & 32) != 0 && false)
351                              problemStatus_ = 0; // say optimal!!
352#if COIN_DEVELOP>2
353
354                         printf("returning at %d\n", __LINE__);
355#endif
356                         return 1; // to primal
357                    }
358               }
359          } else if (!ifValuesPass) {
360               gutsOfSolution(NULL, NULL);
361               // double check
362               if (numberDualInfeasibilities_ || numberPrimalInfeasibilities_)
363                    problemStatus_ = -1;
364          }
365          if (usePrimal) {
366               problemStatus_ = 10;
367#if COIN_DEVELOP>2
368               printf("returning to use primal (no obj) at %d\n", __LINE__);
369#endif
370          }
371          return usePrimal;
372     } else {
373          return 1;
374     }
375}
376void
377ClpSimplexDual::finishSolve(int startFinishOptions)
378{
379     assert (problemStatus_ || !sumPrimalInfeasibilities_);
380
381     // clean up
382     finish(startFinishOptions);
383}
384//#define CLP_REPORT_PROGRESS
385#ifdef CLP_REPORT_PROGRESS
386static int ixxxxxx = 0;
387static int ixxyyyy = 90;
388#endif
389#ifdef CLP_INVESTIGATE_SERIAL
390static int z_reason[7] = {0, 0, 0, 0, 0, 0, 0};
391static int z_thinks = -1;
392#endif
393void
394ClpSimplexDual::gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
395                           ClpDataSave & data)
396{
397#ifdef CLP_INVESTIGATE_SERIAL
398     z_reason[0]++;
399     z_thinks = -1;
400     int nPivots = 9999;
401#endif
402     double largestPrimalError = 0.0;
403     double largestDualError = 0.0;
404     // Start can skip some things in transposeTimes
405     specialOptions_ |= 131072;
406     int lastCleaned = 0; // last time objective or bounds cleaned up
407
408     // This says whether to restore things etc
409     // startup will have factorized so can skip
410     int factorType = 0;
411     // Start check for cycles
412     progress_.startCheck();
413     // Say change made on first iteration
415     // Say last objective infinite
416     //lastObjectiveValue_=-COIN_DBL_MAX;
417     progressFlag_ = 0;
418     /*
419       Status of problem:
420       0 - optimal
421       1 - infeasible
422       2 - unbounded
423       -1 - iterating
424       -2 - factorization wanted
425       -3 - redo checking without factorization
426       -4 - looks infeasible
427     */
428     while (problemStatus_ < 0) {
429          int iRow, iColumn;
430          // clear
431          for (iRow = 0; iRow < 4; iRow++) {
432               rowArray_[iRow]->clear();
433          }
434
435          for (iColumn = 0; iColumn < 2; iColumn++) {
436               columnArray_[iColumn]->clear();
437          }
438
439          // give matrix (and model costs and bounds a chance to be
440          // refreshed (normally null)
441          matrix_->refresh(this);
442          // If getting nowhere - why not give it a kick
443          // does not seem to work too well - do some more work
444          if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
445                    && initialStatus != 10) {
446               perturb();
447               // Can't get here if values pass
448               gutsOfSolution(NULL, NULL);
449               if (handler_->logLevel() > 2) {
450                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
451                              << numberIterations_ << objectiveValue();
452                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
453                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
454                    handler_->printing(sumDualInfeasibilities_ > 0.0)
455                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
456                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
457                                       < numberDualInfeasibilities_)
458                              << numberDualInfeasibilitiesWithoutFree_;
459                    handler_->message() << CoinMessageEol;
460               }
461          }
462          // see if in Cbc etc
463          bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
464#if 0
465          bool gotoPrimal = false;
466          if (inCbcOrOther && numberIterations_ > disasterArea_ + numberRows_ &&
467                    numberDualInfeasibilitiesWithoutFree_ && largestDualError_ > 1.0e-1) {
468               if (!disasterArea_) {
469                    printf("trying all slack\n");
470                    // try all slack basis
471                    allSlackBasis(true);
472                    disasterArea_ = 2 * numberRows_;
473               } else {
474                    printf("going to primal\n");
475                    // go to primal
476                    gotoPrimal = true;
477                    allSlackBasis(true);
478               }
479          }
480#endif
481          bool disaster = false;
482          if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
483               disasterArea_->saveInfo();
484               disaster = true;
485          }
486          // may factorize, checks if problem finished
487          statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
488                                ifValuesPass);
489          largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
490          largestDualError = CoinMax(largestDualError, largestDualError_);
491          if (disaster)
492               problemStatus_ = 3;
493          // If values pass then do easy ones on first time
494          if (ifValuesPass &&
495                    progress_.lastIterationNumber(0) < 0 && saveDuals) {
496               doEasyOnesInValuesPass(saveDuals);
497          }
498
499          // Say good factorization
500          factorType = 1;
501          if (data.sparseThreshold_) {
502               // use default at present
503               factorization_->sparseThreshold(0);
504               factorization_->goSparse();
505          }
506
507          // exit if victory declared
508          if (problemStatus_ >= 0)
509               break;
510
511          // test for maximum iterations
512          if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
513               problemStatus_ = 3;
514               break;
515          }
516          if (ifValuesPass && !saveDuals) {
517               // end of values pass
518               ifValuesPass = 0;
519               int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
520               if (status >= 0) {
521                    problemStatus_ = 5;
522                    secondaryStatus_ = ClpEventHandler::endOfValuesPass;
523                    break;
524               }
525          }
526          // Check event
527          {
528               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
529               if (status >= 0) {
530                    problemStatus_ = 5;
531                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
532                    break;
533               }
534          }
535          // Do iterations
536          int returnCode = whileIterating(saveDuals, ifValuesPass);
537#ifdef CLP_INVESTIGATE_SERIAL
538          nPivots = factorization_->pivots();
539#endif
540          if (!problemStatus_ && factorization_->pivots())
541               computeDuals(NULL); // need to compute duals
542          if (returnCode == -2)
543               factorType = 3;
544     }
545#ifdef CLP_INVESTIGATE_SERIAL
546     // NOTE - can fail if parallel
547     if (z_thinks != -1) {
548          assert (z_thinks < 4);
549          if ((!factorization_->pivots() && nPivots < 20) && z_thinks >= 0 && z_thinks < 2)
550               z_thinks += 4;
551          z_reason[1+z_thinks]++;
552     }
553     if ((z_reason[0] % 1000) == 0) {
554          printf("Reason");
555          for (int i = 0; i < 7; i++)
556               printf(" %d", z_reason[i]);
557          printf("\n");
558     }
559#endif
560     // Stop can skip some things in transposeTimes
561     specialOptions_ &= ~131072;
562     largestPrimalError_ = largestPrimalError;
563     largestDualError_ = largestDualError;
564}
565int
566ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
567{
568     bestObjectiveValue_ = -COIN_DBL_MAX;
569     algorithm_ = -1;
570     moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
571     // save data
572     ClpDataSave data = saveData();
573     double * saveDuals = NULL;
574     int saveDont = dontFactorizePivots_;
575     if ((specialOptions_ & 2048) == 0)
576          dontFactorizePivots_ = 0;
577     else if(!dontFactorizePivots_)
578          dontFactorizePivots_ = 20;
579     if (ifValuesPass) {
580          saveDuals = new double [numberRows_+numberColumns_];
581          CoinMemcpyN(dual_, numberRows_, saveDuals);
582     }
583     if (alphaAccuracy_ != -1.0)
584          alphaAccuracy_ = 1.0;
585     int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
586     // Save so can see if doing after primal
587     int initialStatus = problemStatus_;
588     if (!returnCode && !numberDualInfeasibilities_ &&
589               !numberPrimalInfeasibilities_ && perturbation_ < 101) {
590          returnCode = 1; // to skip gutsOfDual
591          problemStatus_ = 0;
592     }
593
594     if (!returnCode)
595          gutsOfDual(ifValuesPass, saveDuals, initialStatus, data);
596     if (!problemStatus_) {
597          // see if cutoff reached
598          double limit = 0.0;
599          getDblParam(ClpDualObjectiveLimit, limit);
600          if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
601                    limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
602               // actually infeasible on objective
603               problemStatus_ = 1;
604               secondaryStatus_ = 1;
605          }
606     }
607     // If infeasible but primal errors - try dual
608     if (problemStatus_==1 && numberPrimalInfeasibilities_) {
609       bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
610       double factor = (!inCbcOrOther) ? 1.0 : 0.3;
611       double averageInfeasibility = sumPrimalInfeasibilities_/
612         static_cast<double>(numberPrimalInfeasibilities_);
613       if (averageInfeasibility<factor*largestPrimalError_)
614         problemStatus_= 10;
615     }
616
617     if (problemStatus_ == 10)
618          startFinishOptions |= 1;
619     finishSolve(startFinishOptions);
620     delete [] saveDuals;
621
622     // Restore any saved stuff
623     restoreData(data);
624     dontFactorizePivots_ = saveDont;
625     if (problemStatus_ == 3)
626          objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
627     return problemStatus_;
628}
629// old way
630#if 0
631int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
632{
633     algorithm_ = -1;
634
635     // save data
636     ClpDataSave data = saveData();
637     // Save so can see if doing after primal
638     int initialStatus = problemStatus_;
639
640     // If values pass then save given duals round check solution
641     double * saveDuals = NULL;
642     if (ifValuesPass) {
643          saveDuals = new double [numberRows_+numberColumns_];
644          CoinMemcpyN(dual_, numberRows_, saveDuals);
645     }
646     // sanity check
647     // initialize - no values pass and algorithm_ is -1
648     // put in standard form (and make row copy)
649     // create modifiable copies of model rim and do optional scaling
650     // If problem looks okay
651     // Do initial factorization
652     // If user asked for perturbation - do it
653     if (!startup(0, startFinishOptions)) {
654          // looks okay
655          // Superbasic variables not allowed
656          // If values pass then scale pi
657          if (ifValuesPass) {
658               if (problemStatus_ && perturbation_ < 100)
659                    perturb();
660               int i;
661               if (scalingFlag_ > 0) {
662                    for (i = 0; i < numberRows_; i++) {
663                         dual_[i] = saveDuals[i] * inverseRowScale_[i];
664                    }
665               } else {
666                    CoinMemcpyN(saveDuals, numberRows_, dual_);
667               }
668               // now create my duals
669               for (i = 0; i < numberRows_; i++) {
670                    // slack
671                    double value = dual_[i];
672                    value += rowObjectiveWork_[i];
673                    saveDuals[i+numberColumns_] = value;
674               }
675               CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
676               transposeTimes(-1.0, dual_, saveDuals);
677               // make reduced costs okay
678               for (i = 0; i < numberColumns_; i++) {
679                    if (getStatus(i) == atLowerBound) {
680                         if (saveDuals[i] < 0.0) {
681                              //if (saveDuals[i]<-1.0e-3)
682                              //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
683                              saveDuals[i] = 0.0;
684                         }
685                    } else if (getStatus(i) == atUpperBound) {
686                         if (saveDuals[i] > 0.0) {
687                              //if (saveDuals[i]>1.0e-3)
688                              //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
689                              saveDuals[i] = 0.0;
690                         }
691                    }
692               }
693               CoinMemcpyN(saveDuals, numberColumns_ + numberRows_, dj_);
694               // set up possible ones
695               for (i = 0; i < numberRows_ + numberColumns_; i++)
696                    clearPivoted(i);
697               int iRow;
698               for (iRow = 0; iRow < numberRows_; iRow++) {
699                    int iPivot = pivotVariable_[iRow];
700                    if (fabs(saveDuals[iPivot]) > dualTolerance_) {
701                         if (getStatus(iPivot) != isFree)
702                              setPivoted(iPivot);
703                    }
704               }
705          } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
706               // set up possible ones
707               for (int i = 0; i < numberRows_ + numberColumns_; i++)
708                    clearPivoted(i);
709               int iRow;
710               for (iRow = 0; iRow < numberRows_; iRow++) {
711                    int iPivot = pivotVariable_[iRow];
712                    if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
713                         setPivoted(iPivot);
714                    }
715               }
716          }
717
718          double objectiveChange;
719          numberFake_ = 0; // Number of variables at fake bounds
720          numberChanged_ = 0; // Number of variables with changed costs
721          changeBounds(1, NULL, objectiveChange);
722
723          int lastCleaned = 0; // last time objective or bounds cleaned up
724
725          if (!ifValuesPass) {
726               // Check optimal
727               if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
728                    problemStatus_ = 0;
729          }
730          if (problemStatus_ < 0 && perturbation_ < 100) {
731               perturb();
732               // Can't get here if values pass
733               gutsOfSolution(NULL, NULL);
734               if (handler_->logLevel() > 2) {
735                    handler_->message(CLP_SIMPLEX_STATUS, messages_)
736                              << numberIterations_ << objectiveValue();
737                    handler_->printing(sumPrimalInfeasibilities_ > 0.0)
738                              << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
739                    handler_->printing(sumDualInfeasibilities_ > 0.0)
740                              << sumDualInfeasibilities_ << numberDualInfeasibilities_;
741                    handler_->printing(numberDualInfeasibilitiesWithoutFree_
742                                       < numberDualInfeasibilities_)
743                              << numberDualInfeasibilitiesWithoutFree_;
744                    handler_->message() << CoinMessageEol;
745               }
746          }
747
748          // This says whether to restore things etc
749          // startup will have factorized so can skip
750          int factorType = 0;
751          // Start check for cycles
752          progress_.startCheck();
753          // Say change made on first iteration
755          /*
756            Status of problem:
757            0 - optimal
758            1 - infeasible
759            2 - unbounded
760            -1 - iterating
761            -2 - factorization wanted
762            -3 - redo checking without factorization
763            -4 - looks infeasible
764          */
765          while (problemStatus_ < 0) {
766               int iRow, iColumn;
767               // clear
768               for (iRow = 0; iRow < 4; iRow++) {
769                    rowArray_[iRow]->clear();
770               }
771
772               for (iColumn = 0; iColumn < 2; iColumn++) {
773                    columnArray_[iColumn]->clear();
774               }
775
776               // give matrix (and model costs and bounds a chance to be
777               // refreshed (normally null)
778               matrix_->refresh(this);
779               // If getting nowhere - why not give it a kick
780               // does not seem to work too well - do some more work
781               if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
782                         && initialStatus != 10) {
783                    perturb();
784                    // Can't get here if values pass
785                    gutsOfSolution(NULL, NULL);
786                    if (handler_->logLevel() > 2) {
787                         handler_->message(CLP_SIMPLEX_STATUS, messages_)
788                                   << numberIterations_ << objectiveValue();
789                         handler_->printing(sumPrimalInfeasibilities_ > 0.0)
790                                   << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
791                         handler_->printing(sumDualInfeasibilities_ > 0.0)
792                                   << sumDualInfeasibilities_ << numberDualInfeasibilities_;
793                         handler_->printing(numberDualInfeasibilitiesWithoutFree_
794                                            < numberDualInfeasibilities_)
795                                   << numberDualInfeasibilitiesWithoutFree_;
796                         handler_->message() << CoinMessageEol;
797                    }
798               }
799               // may factorize, checks if problem finished
800               statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
801                                     ifValuesPass);
802               // If values pass then do easy ones on first time
803               if (ifValuesPass &&
804                         progress_.lastIterationNumber(0) < 0 && saveDuals) {
805                    doEasyOnesInValuesPass(saveDuals);
806               }
807
808               // Say good factorization
809               factorType = 1;
810               if (data.sparseThreshold_) {
811                    // use default at present
812                    factorization_->sparseThreshold(0);
813                    factorization_->goSparse();
814               }
815
816               // exit if victory declared
817               if (problemStatus_ >= 0)
818                    break;
819
820               // test for maximum iterations
821               if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
822                    problemStatus_ = 3;
823                    break;
824               }
825               if (ifValuesPass && !saveDuals) {
826                    // end of values pass
827                    ifValuesPass = 0;
828                    int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
829                    if (status >= 0) {
830                         problemStatus_ = 5;
831                         secondaryStatus_ = ClpEventHandler::endOfValuesPass;
832                         break;
833                    }
834               }
835               // Check event
836               {
837                    int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
838                    if (status >= 0) {
839                         problemStatus_ = 5;
840                         secondaryStatus_ = ClpEventHandler::endOfFactorization;
841                         break;
842                    }
843               }
844               // Do iterations
845               whileIterating(saveDuals, ifValuesPass);
846          }
847     }
848
849     assert (problemStatus_ || !sumPrimalInfeasibilities_);
850
851     // clean up
852     finish(startFinishOptions);
853     delete [] saveDuals;
854
855     // Restore any saved stuff
856     restoreData(data);
857     return problemStatus_;
858}
859#endif
860//#define CHECK_ACCURACY
861#ifdef CHECK_ACCURACY
862static double zzzzzz[100000];
863#endif
864/* Reasons to come out:
865   -1 iterations etc
866   -2 inaccuracy
867   -3 slight inaccuracy (and done iterations)
868   +0 looks optimal (might be unbounded - but we will investigate)
869   +1 looks infeasible
870   +3 max iterations
871 */
872int
873ClpSimplexDual::whileIterating(double * & givenDuals, int ifValuesPass)
874{
875#ifdef CLP_INVESTIGATE_SERIAL
876     z_thinks = -1;
877#endif
878#ifdef CLP_DEBUG
879     int debugIteration = -1;
880#endif
881     {
882          int i;
883          for (i = 0; i < 4; i++) {
884               rowArray_[i]->clear();
885          }
886          for (i = 0; i < 2; i++) {
887               columnArray_[i]->clear();
888          }
889     }
890#ifdef CLP_REPORT_PROGRESS
891     double * savePSol = new double [numberRows_+numberColumns_];
892     double * saveDj = new double [numberRows_+numberColumns_];
893     double * saveCost = new double [numberRows_+numberColumns_];
894     unsigned char * saveStat = new unsigned char [numberRows_+numberColumns_];
895#endif
896     // if can't trust much and long way from optimal then relax
897     if (largestPrimalError_ > 10.0)
898          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
899     else
900          factorization_->relaxAccuracyCheck(1.0);
901     // status stays at -1 while iterating, >=0 finished, -2 to invert
902     // status -3 to go to top without an invert
903     int returnCode = -1;
904     double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
905
906#if 0
907     // compute average infeasibility for backward test
908     double averagePrimalInfeasibility = sumPrimalInfeasibilities_ /
909                                         ((double ) (numberPrimalInfeasibilities_ + 1));
910#endif
911
912     // Get dubious weights
913     CoinBigIndex * dubiousWeights = NULL;
914#ifdef DUBIOUS_WEIGHTS
915     factorization_->getWeights(rowArray_[0]->getIndices());
916     dubiousWeights = matrix_->dubiousWeights(this, rowArray_[0]->getIndices());
917#endif
918     // If values pass then get list of candidates
919     int * candidateList = NULL;
920     int numberCandidates = 0;
921#ifdef CLP_DEBUG
922     bool wasInValuesPass = (givenDuals != NULL);
923#endif
924     int candidate = -1;
925     if (givenDuals) {
926          assert (ifValuesPass);
927          ifValuesPass = 1;
928          candidateList = new int[numberRows_];
929          // move reduced costs across
930          CoinMemcpyN(givenDuals, numberRows_ + numberColumns_, dj_);
931          int iRow;
932          for (iRow = 0; iRow < numberRows_; iRow++) {
933               int iPivot = pivotVariable_[iRow];
934               if (flagged(iPivot))
935                    continue;
936               if (fabs(dj_[iPivot]) > dualTolerance_) {
937                    // for now safer to ignore free ones
938                    if (lower_[iPivot] > -1.0e50 || upper_[iPivot] < 1.0e50)
939                         if (pivoted(iPivot))
940                              candidateList[numberCandidates++] = iRow;
941               } else {
942                    clearPivoted(iPivot);
943               }
944          }
945          // and set first candidate
946          if (!numberCandidates) {
947               delete [] candidateList;
948               delete [] givenDuals;
949               givenDuals = NULL;
950               candidateList = NULL;
951               int iRow;
952               for (iRow = 0; iRow < numberRows_; iRow++) {
953                    int iPivot = pivotVariable_[iRow];
954                    clearPivoted(iPivot);
955               }
956          }
957     } else {
958          assert (!ifValuesPass);
959     }
960#ifdef CHECK_ACCURACY
961     {
962          if (numberIterations_) {
963               int il = -1;
964               double largest = 1.0e-1;
965               int ilnb = -1;
966               double largestnb = 1.0e-8;
967               for (int i = 0; i < numberRows_ + numberColumns_; i++) {
968                    double diff = fabs(solution_[i] - zzzzzz[i]);
969                    if (diff > largest) {
970                         largest = diff;
971                         il = i;
972                    }
973                    if (getColumnStatus(i) != basic) {
974                         if (diff > largestnb) {
975                              largestnb = diff;
976                              ilnb = i;
977                         }
978                    }
979               }
980               if (il >= 0 && ilnb < 0)
981                    printf("largest diff of %g at %d, nonbasic %g at %d\n",
982                           largest, il, largestnb, ilnb);
983          }
984     }
985#endif
986     while (problemStatus_ == -1) {
987          //if (numberIterations_>=101624)
988          //resetFakeBounds(-1);
989#ifdef CLP_DEBUG
990          if (givenDuals) {
991               double value5 = 0.0;
992               int i;
993               for (i = 0; i < numberRows_ + numberColumns_; i++) {
994                    if (dj_[i] < -1.0e-6)
995                         if (upper_[i] < 1.0e20)
996                              value5 += dj_[i] * upper_[i];
997                         else
998                              printf("bad dj %g on %d with large upper status %d\n",
999                                     dj_[i], i, status_[i] & 7);
1000                    else if (dj_[i] > 1.0e-6)
1001                         if (lower_[i] > -1.0e20)
1002                              value5 += dj_[i] * lower_[i];
1003                         else
1004                              printf("bad dj %g on %d with large lower status %d\n",
1005                                     dj_[i], i, status_[i] & 7);
1006               }
1007               printf("Values objective Value %g\n", value5);
1008          }
1009          if ((handler_->logLevel() & 32) && wasInValuesPass) {
1010               double value5 = 0.0;
1011               int i;
1012               for (i = 0; i < numberRows_ + numberColumns_; i++) {
1013                    if (dj_[i] < -1.0e-6)
1014                         if (upper_[i] < 1.0e20)
1015                              value5 += dj_[i] * upper_[i];
1016                         else if (dj_[i] > 1.0e-6)
1017                              if (lower_[i] > -1.0e20)
1018                                   value5 += dj_[i] * lower_[i];
1019               }
1020               printf("Values objective Value %g\n", value5);
1021               {
1022                    int i;
1023                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
1024                         int iSequence = i;
1025                         double oldValue;
1026
1027                         switch(getStatus(iSequence)) {
1028
1029                         case basic:
1030                         case ClpSimplex::isFixed:
1031                              break;
1032                         case isFree:
1033                         case superBasic:
1034                              abort();
1035                              break;
1036                         case atUpperBound:
1037                              oldValue = dj_[iSequence];
1038                              //assert (oldValue<=tolerance);
1039                              assert (fabs(solution_[iSequence] - upper_[iSequence]) < 1.0e-7);
1040                              break;
1041                         case atLowerBound:
1042                              oldValue = dj_[iSequence];
1043                              //assert (oldValue>=-tolerance);
1044                              assert (fabs(solution_[iSequence] - lower_[iSequence]) < 1.0e-7);
1045                              break;
1046                         }
1047                    }
1048               }
1049          }
1050#endif
1051#ifdef CLP_DEBUG
1052          {
1053               int i;
1054               for (i = 0; i < 4; i++) {
1055                    rowArray_[i]->checkClear();
1056               }
1057               for (i = 0; i < 2; i++) {
1058                    columnArray_[i]->checkClear();
1059               }
1060          }
1061#endif
1062#if CLP_DEBUG>2
1063          // very expensive
1064          if (numberIterations_ > 3063 && numberIterations_ < 30700) {
1065               //handler_->setLogLevel(63);
1066               double saveValue = objectiveValue_;
1067               double * saveRow1 = new double[numberRows_];
1068               double * saveRow2 = new double[numberRows_];
1069               CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
1070               CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
1071               double * saveColumn1 = new double[numberColumns_];
1072               double * saveColumn2 = new double[numberColumns_];
1073               CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
1074               CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
1075               gutsOfSolution(NULL, NULL);
1076               printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
1077                      numberIterations_,
1078                      saveValue, objectiveValue_, sumDualInfeasibilities_);
1079               if (saveValue > objectiveValue_ + 1.0e-2)
1081               CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
1082               CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
1083               CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
1084               CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
1085               delete [] saveRow1;
1086               delete [] saveRow2;
1087               delete [] saveColumn1;
1088               delete [] saveColumn2;
1089               objectiveValue_ = saveValue;
1090          }
1091#endif
1092#if 0
1093          //    if (factorization_->pivots()){
1094          {
1095               int iPivot;
1096               double * array = rowArray_[3]->denseVector();
1097               int i;
1098               for (iPivot = 0; iPivot < numberRows_; iPivot++) {
1099                    int iSequence = pivotVariable_[iPivot];
1100                    unpack(rowArray_[3], iSequence);
1101                    factorization_->updateColumn(rowArray_[2], rowArray_[3]);
1102                    assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
1103                    array[iPivot] = 0.0;
1104                    for (i = 0; i < numberRows_; i++)
1105                         assert (fabs(array[i]) < 1.0e-4);
1106                    rowArray_[3]->clear();
1107               }
1108          }
1109#endif
1110#ifdef CLP_DEBUG
1111          {
1112               int iSequence, number = numberRows_ + numberColumns_;
1113               for (iSequence = 0; iSequence < number; iSequence++) {
1114                    double lowerValue = lower_[iSequence];
1115                    double upperValue = upper_[iSequence];
1116                    double value = solution_[iSequence];
1117                    if(getStatus(iSequence) != basic && getStatus(iSequence) != isFree) {
1118                         assert(lowerValue > -1.0e20);
1119                         assert(upperValue < 1.0e20);
1120                    }
1121                    switch(getStatus(iSequence)) {
1122
1123                    case basic:
1124                         break;
1125                    case isFree:
1126                    case superBasic:
1127                         break;
1128                    case atUpperBound:
1129                         assert (fabs(value - upperValue) <= primalTolerance_) ;
1130                         break;
1131                    case atLowerBound:
1132                    case ClpSimplex::isFixed:
1133                         assert (fabs(value - lowerValue) <= primalTolerance_) ;
1134                         break;
1135                    }
1136               }
1137          }
1138          if(numberIterations_ == debugIteration) {
1139               printf("dodgy iteration coming up\n");
1140          }
1141#endif
1142#if 0
1143          printf("checking nz\n");
1144          for (int i = 0; i < 3; i++) {
1145               if (!rowArray_[i]->getNumElements())
1146                    rowArray_[i]->checkClear();
1147          }
1148#endif
1149          // choose row to go out
1150          // dualRow will go to virtual row pivot choice algorithm
1151          // make sure values pass off if it should be
1152          if (numberCandidates)
1153               candidate = candidateList[--numberCandidates];
1154          else
1155               candidate = -1;
1156          dualRow(candidate);
1157          if (pivotRow_ >= 0) {
1158               // we found a pivot row
1159               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
1160                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
1161                              << pivotRow_
1162                              << CoinMessageEol;
1163               }
1164               // check accuracy of weights
1165               dualRowPivot_->checkAccuracy();
1166               // Get good size for pivot
1167               // Allow first few iterations to take tiny
1168               double acceptablePivot = 1.0e-1 * acceptablePivot_;
1169               if (numberIterations_ > 100)
1170                    acceptablePivot = acceptablePivot_;
1171               if (factorization_->pivots() > 10 ||
1172                         (factorization_->pivots() && saveSumDual))
1173                    acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1174               else if (factorization_->pivots() > 5)
1175                    acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1176               else if (factorization_->pivots())
1177                    acceptablePivot = acceptablePivot_; // relax
1178               // But factorizations complain if <1.0e-8
1179               //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
1180               double bestPossiblePivot = 1.0;
1181               // get sign for finding row of tableau
1182               if (candidate < 0) {
1183                    // normal iteration
1184                    // create as packed
1185                    double direction = directionOut_;
1186                    rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1187                    factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1188                    // Allow to do dualColumn0
1190                         spareIntArray_[0] = 1;
1191                    spareDoubleArray_[0] = acceptablePivot;
1192                    rowArray_[3]->clear();
1193                    sequenceIn_ = -1;
1194                    // put row of tableau in rowArray[0] and columnArray[0]
1195                    assert (!rowArray_[1]->getNumElements());
1196                    if (!scaledMatrix_) {
1197                         if ((moreSpecialOptions_ & 8) != 0 && !rowScale_)
1198                              spareIntArray_[0] = 1;
1199                         matrix_->transposeTimes(this, -1.0,
1200                                                 rowArray_[0], rowArray_[1], columnArray_[0]);
1201                    } else {
1202                         double * saveR = rowScale_;
1203                         double * saveC = columnScale_;
1204                         rowScale_ = NULL;
1205                         columnScale_ = NULL;
1206                         if ((moreSpecialOptions_ & 8) != 0)
1207                              spareIntArray_[0] = 1;
1208                         scaledMatrix_->transposeTimes(this, -1.0,
1209                                                       rowArray_[0], rowArray_[1], columnArray_[0]);
1210                         rowScale_ = saveR;
1211                         columnScale_ = saveC;
1212                    }
1213#ifdef CLP_REPORT_PROGRESS
1214                    memcpy(savePSol, solution_, (numberColumns_ + numberRows_)*sizeof(double));
1215                    memcpy(saveDj, dj_, (numberColumns_ + numberRows_)*sizeof(double));
1216                    memcpy(saveCost, cost_, (numberColumns_ + numberRows_)*sizeof(double));
1217                    memcpy(saveStat, status_, (numberColumns_ + numberRows_)*sizeof(char));
1218#endif
1219                    // do ratio test for normal iteration
1220                    bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
1221                                                   columnArray_[1], acceptablePivot, dubiousWeights);
1222               } else {
1223                    // Make sure direction plausible
1224                    CoinAssert (upperOut_ < 1.0e50 || lowerOut_ > -1.0e50);
1225                    // If in integer cleanup do direction using duals
1226                    // may be wrong way round
1227                    if(ifValuesPass == 2) {
1228                         if (dual_[pivotRow_] > 0.0) {
1229                              // this will give a -1 in pivot row (as slacks are -1.0)
1230                              directionOut_ = 1;
1231                         } else {
1232                              directionOut_ = -1;
1233                         }
1234                    }
1235                    if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
1236                         if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
1237                              directionOut_ = 1;
1238                    } else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
1239                         if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
1240                              directionOut_ = -1;
1241                    }
1242                    double direction = directionOut_;
1243                    rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1244                    factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1245                    // put row of tableau in rowArray[0] and columnArray[0]
1246                    if (!scaledMatrix_) {
1247                         matrix_->transposeTimes(this, -1.0,
1248                                                 rowArray_[0], rowArray_[3], columnArray_[0]);
1249                    } else {
1250                         double * saveR = rowScale_;
1251                         double * saveC = columnScale_;
1252                         rowScale_ = NULL;
1253                         columnScale_ = NULL;
1254                         scaledMatrix_->transposeTimes(this, -1.0,
1255                                                       rowArray_[0], rowArray_[3], columnArray_[0]);
1256                         rowScale_ = saveR;
1257                         columnScale_ = saveC;
1258                    }
1259                    acceptablePivot *= 10.0;
1260                    // do ratio test
1261                    if (ifValuesPass == 1) {
1262                         checkPossibleValuesMove(rowArray_[0], columnArray_[0],
1263                                                 acceptablePivot);
1264                    } else {
1265                         checkPossibleCleanup(rowArray_[0], columnArray_[0],
1266                                              acceptablePivot);
1267                         if (sequenceIn_ < 0) {
1268                              rowArray_[0]->clear();
1269                              columnArray_[0]->clear();
1270                              continue; // can't do anything
1271                         }
1272                    }
1273
1274                    // recompute true dualOut_
1275                    if (directionOut_ < 0) {
1276                         dualOut_ = valueOut_ - upperOut_;
1277                    } else {
1278                         dualOut_ = lowerOut_ - valueOut_;
1279                    }
1280                    // check what happened if was values pass
1281                    // may want to move part way i.e. movement
1282                    bool normalIteration = (sequenceIn_ != sequenceOut_);
1283
1284                    clearPivoted(sequenceOut_);  // make sure won't be done again
1285                    // see if end of values pass
1286                    if (!numberCandidates) {
1287                         int iRow;
1288                         delete [] candidateList;
1289                         delete [] givenDuals;
1290                         candidate = -2; // -2 signals end
1291                         givenDuals = NULL;
1292                         candidateList = NULL;
1293                         ifValuesPass = 1;
1294                         for (iRow = 0; iRow < numberRows_; iRow++) {
1295                              int iPivot = pivotVariable_[iRow];
1296                              //assert (fabs(dj_[iPivot]),1.0e-5);
1297                              clearPivoted(iPivot);
1298                         }
1299                    }
1300                    if (!normalIteration) {
1301                         //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1302                         //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1303                         updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1304                         if (candidate == -2)
1305                              problemStatus_ = -2;
1306                         continue; // skip rest of iteration
1307                    } else {
1308                         // recompute dualOut_
1309                         if (directionOut_ < 0) {
1310                              dualOut_ = valueOut_ - upperOut_;
1311                         } else {
1312                              dualOut_ = lowerOut_ - valueOut_;
1313                         }
1314                    }
1315               }
1316               if (sequenceIn_ >= 0) {
1317                    // normal iteration
1318                    // update the incoming column
1319                    double btranAlpha = -alpha_ * directionOut_; // for check
1320                    unpackPacked(rowArray_[1]);
1321                    // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1322                    // and update dual weights (can do in parallel - with extra array)
1323                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
1324                                                          rowArray_[2],
1325                                                          rowArray_[3],
1326                                                          rowArray_[1]);
1327                    // see if update stable
1328#ifdef CLP_DEBUG
1329                    if ((handler_->logLevel() & 32))
1330                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
1331#endif
1332                    double checkValue = 1.0e-7;
1333                    // if can't trust much and long way from optimal then relax
1334                    if (largestPrimalError_ > 10.0)
1335                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1336                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1337                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1338                         handler_->message(CLP_DUAL_CHECK, messages_)
1339                                   << btranAlpha
1340                                   << alpha_
1341                                   << CoinMessageEol;
1342                         if (factorization_->pivots()) {
1343                              dualRowPivot_->unrollWeights();
1344                              problemStatus_ = -2; // factorize now
1345                              rowArray_[0]->clear();
1346                              rowArray_[1]->clear();
1347                              columnArray_[0]->clear();
1348                              returnCode = -2;
1349                              break;
1350                         } else {
1351                              // take on more relaxed criterion
1352                              double test;
1353                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
1354                                   test = 1.0e-1 * fabs(alpha_);
1355                              else
1356                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
1357                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1358                                        fabs(btranAlpha - alpha_) > test) {
1359                                   dualRowPivot_->unrollWeights();
1360                                   // need to reject something
1361                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
1362                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
1363                                             << x << sequenceWithin(sequenceOut_)
1364                                             << CoinMessageEol;
1365#ifdef COIN_DEVELOP
1366                                   printf("flag a %g %g\n", btranAlpha, alpha_);
1367#endif
1368                                   //#define FEB_TRY
1369#ifdef FEB_TRY
1370                                   // Make safer?
1371                                   factorization_->saferTolerances (1.0e-15, -1.03);
1372#endif
1373                                   setFlagged(sequenceOut_);
1375                                   lastBadIteration_ = numberIterations_; // say be more cautious
1376                                   rowArray_[0]->clear();
1377                                   rowArray_[1]->clear();
1378                                   columnArray_[0]->clear();
1379                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
1380                                        //printf("I think should declare infeasible\n");
1381                                        problemStatus_ = 1;
1382                                        returnCode = 1;
1383                                        break;
1384                                   }
1385                                   continue;
1386                              }
1387                         }
1388                    }
1389                    // update duals BEFORE replaceColumn so can do updateColumn
1390                    double objectiveChange = 0.0;
1391                    // do duals first as variables may flip bounds
1392                    // rowArray_[0] and columnArray_[0] may have flips
1393                    // so use rowArray_[3] for work array from here on
1394                    int nswapped = 0;
1395                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1396                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1397                    if (candidate == -1) {
1398                         // make sure incoming doesn't count
1399                         Status saveStatus = getStatus(sequenceIn_);
1400                         setStatus(sequenceIn_, basic);
1401                         nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
1402                                                      rowArray_[2], theta_,
1403                                                      objectiveChange, false);
1404                         setStatus(sequenceIn_, saveStatus);
1405                    } else {
1406                         updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1407                    }
1408                    double oldDualOut = dualOut_;
1409                    // which will change basic solution
1410                    if (nswapped) {
1411                         if (rowArray_[2]->getNumElements()) {
1412                              factorization_->updateColumn(rowArray_[3], rowArray_[2]);
1413                              dualRowPivot_->updatePrimalSolution(rowArray_[2],
1414                                                                  1.0, objectiveChange);
1415                         }
1416                         // recompute dualOut_
1417                         valueOut_ = solution_[sequenceOut_];
1418                         if (directionOut_ < 0) {
1419                              dualOut_ = valueOut_ - upperOut_;
1420                         } else {
1421                              dualOut_ = lowerOut_ - valueOut_;
1422                         }
1423#if 0
1424                         if (dualOut_ < 0.0) {
1425#ifdef CLP_DEBUG
1426                              if (handler_->logLevel() & 32) {
1427                                   printf(" dualOut_ %g %g save %g\n", dualOut_, averagePrimalInfeasibility, saveDualOut);
1428                                   printf("values %g %g %g %g %g %g %g\n", lowerOut_, valueOut_, upperOut_,
1429                                          objectiveChange,);
1430                              }
1431#endif
1432                              if (upperOut_ == lowerOut_)
1433                                   dualOut_ = 0.0;
1434                         }
1435                         if(dualOut_ < -CoinMax(1.0e-12 * averagePrimalInfeasibility, 1.0e-8)
1436                                   && factorization_->pivots() > 100 &&
1437                                   getStatus(sequenceIn_) != isFree) {
1438                              // going backwards - factorize
1439                              dualRowPivot_->unrollWeights();
1440                              problemStatus_ = -2; // factorize now
1441                              returnCode = -2;
1442                              break;
1443                         }
1444#endif
1445                    }
1446                    // amount primal will move
1447                    double movement = -dualOut_ * directionOut_ / alpha_;
1448                    double movementOld = oldDualOut * directionOut_ / alpha_;
1449                    // so objective should increase by fabs(dj)*movement
1450                    // but we already have objective change - so check will be good
1451                    if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
1452#ifdef CLP_DEBUG
1453                         if (handler_->logLevel() & 32)
1454                              printf("movement %g, swap change %g, rest %g  * %g\n",
1455                                     objectiveChange + fabs(movement * dualIn_),
1456                                     objectiveChange, movement, dualIn_);
1457#endif
1458                         if(factorization_->pivots()) {
1459                              // going backwards - factorize
1460                              dualRowPivot_->unrollWeights();
1461                              problemStatus_ = -2; // factorize now
1462                              returnCode = -2;
1463                              break;
1464                         }
1465                    }
1466                    // if stable replace in basis
1468                                       rowArray_[2],
1469                                       rowArray_[1],
1470                                       pivotRow_,
1471                                       alpha_,
1472                                       (moreSpecialOptions_ & 16) != 0,
1473                                       acceptablePivot);
1474                    // If looks like bad pivot - refactorize
1475                    if (fabs(dualOut_) > 1.0e50)
1477                    // if no pivots, bad update but reasonable alpha - take and invert
1478                    if (updateStatus == 2 &&
1479                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
1482                         // slight error
1483                         if (factorization_->pivots() > 5 || updateStatus == 4) {
1484                              problemStatus_ = -2; // factorize now
1485                              returnCode = -3;
1486                         }
1487                    } else if (updateStatus == 2) {
1488                         // major error
1489                         dualRowPivot_->unrollWeights();
1490                         // later we may need to unwind more e.g. fake bounds
1491                         if (factorization_->pivots() &&
1492                                   ((moreSpecialOptions_ & 16) == 0 || factorization_->pivots() > 4)) {
1493                              problemStatus_ = -2; // factorize now
1494                              returnCode = -2;
1495                              moreSpecialOptions_ |= 16;
1496                              break;
1497                         } else {
1498                              // need to reject something
1499                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
1500                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
1501                                        << x << sequenceWithin(sequenceOut_)
1502                                        << CoinMessageEol;
1503#ifdef COIN_DEVELOP
1504                              printf("flag b %g\n", alpha_);
1505#endif
1506                              setFlagged(sequenceOut_);
1508                              lastBadIteration_ = numberIterations_; // say be more cautious
1509                              rowArray_[0]->clear();
1510                              rowArray_[1]->clear();
1511                              columnArray_[0]->clear();
1512                              // make sure dual feasible
1513                              // look at all rows and columns
1514                              double objectiveChange = 0.0;
1515                              updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
1516                                                0.0, objectiveChange, true);
1517                              rowArray_[1]->clear();
1518                              columnArray_[0]->clear();
1519                              continue;
1520                         }
1521                    } else if (updateStatus == 3) {
1522                         // out of memory
1523                         // increase space if not many iterations
1524                         if (factorization_->pivots() <
1525                                   0.5 * factorization_->maximumPivots() &&
1526                                   factorization_->pivots() < 200)
1527                              factorization_->areaFactor(
1528                                   factorization_->areaFactor() * 1.1);
1529                         problemStatus_ = -2; // factorize now
1530                    } else if (updateStatus == 5) {
1531                         problemStatus_ = -2; // factorize now
1532                    }
1533                    // update primal solution
1534                    if (theta_ < 0.0 && candidate == -1) {
1535#ifdef CLP_DEBUG
1536                         if (handler_->logLevel() & 32)
1537                              printf("negative theta %g\n", theta_);
1538#endif
1539                         theta_ = 0.0;
1540                    }
1541                    // do actual flips
1542                    flipBounds(rowArray_[0], columnArray_[0]);
1543                    //rowArray_[1]->expand();
1544                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
1545                                                        movement,
1546                                                        objectiveChange);
1547#ifdef CLP_DEBUG
1548                    double oldobj = objectiveValue_;
1549#endif
1550                    // modify dualout
1551                    dualOut_ /= alpha_;
1552                    dualOut_ *= -directionOut_;
1553                    //setStatus(sequenceIn_,basic);
1554                    dj_[sequenceIn_] = 0.0;
1555                    double oldValue = valueIn_;
1556                    if (directionIn_ == -1) {
1557                         // as if from upper bound
1558                         valueIn_ = upperIn_ + dualOut_;
1559                    } else {
1560                         // as if from lower bound
1561                         valueIn_ = lowerIn_ + dualOut_;
1562                    }
1563                    objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
1564                    // outgoing
1565                    // set dj to zero unless values pass
1566                    if (directionOut_ > 0) {
1567                         valueOut_ = lowerOut_;
1568                         if (candidate == -1)
1569                              dj_[sequenceOut_] = theta_;
1570                    } else {
1571                         valueOut_ = upperOut_;
1572                         if (candidate == -1)
1573                              dj_[sequenceOut_] = -theta_;
1574                    }
1575                    solution_[sequenceOut_] = valueOut_;
1576                    int whatNext = housekeeping(objectiveChange);
1577#ifdef CLP_REPORT_PROGRESS
1578                    if (ixxxxxx > ixxyyyy - 5) {
1579                         handler_->setLogLevel(63);
1580                         int nTotal = numberColumns_ + numberRows_;
1581                         double oldObj = 0.0;
1582                         double newObj = 0.0;
1583                         for (int i = 0; i < nTotal; i++) {
1584                              if (savePSol[i])
1585                                   oldObj += savePSol[i] * saveCost[i];
1586                              if (solution_[i])
1587                                   newObj += solution_[i] * cost_[i];
1588                              bool printIt = false;
1589                              if (cost_[i] != saveCost[i])
1590                                   printIt = true;
1591                              if (status_[i] != saveStat[i])
1592                                   printIt = true;
1593                              if (printIt)
1594                                   printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
1595                                          i, saveStat[i], saveCost[i], savePSol[i],
1596                                          status_[i], cost_[i], solution_[i]);
1597                              // difference
1598                              savePSol[i] = solution_[i] - savePSol[i];
1599                         }
1600                         printf("pivots %d, old obj %g new %g\n",
1601                                factorization_->pivots(),
1602                                oldObj, newObj);
1603                         memset(saveDj, 0, numberRows_ * sizeof(double));
1604                         times(1.0, savePSol, saveDj);
1605                         double largest = 1.0e-6;
1606                         int k = -1;
1607                         for (int i = 0; i < numberRows_; i++) {
1608                              saveDj[i] -= savePSol[i+numberColumns_];
1609                              if (fabs(saveDj[i]) > largest) {
1610                                   largest = fabs(saveDj[i]);
1611                                   k = i;
1612                              }
1613                         }
1614                         if (k >= 0)
1615                              printf("Not null %d %g\n", k, largest);
1616                    }
1617#endif
1618#ifdef VUB
1619                    {
1620                         if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) || toVub[sequenceIn_] >= 0 ||
1621                                   (sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) || toVub[sequenceOut_] >= 0) {
1622                              int inSequence = sequenceIn_;
1623                              int inVub = -1;
1624                              if (sequenceIn_ < numberColumns_)
1625                                   inVub = vub[sequenceIn_];
1626                              int inBack = toVub[inSequence];
1627                              int inSlack = -1;
1628                              if (inSequence >= numberColumns_ && inBack >= 0) {
1629                                   inSlack = inSequence - numberColumns_;
1630                                   inSequence = inBack;
1631                                   inBack = toVub[inSequence];
1632                              }
1633                              if (inVub >= 0)
1634                                   printf("Vub %d in ", inSequence);
1635                              if (inBack >= 0 && inSlack < 0)
1636                                   printf("%d (descendent of %d) in ", inSequence, inBack);
1637                              if (inSlack >= 0)
1638                                   printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
1639                              int outSequence = sequenceOut_;
1640                              int outVub = -1;
1641                              if (sequenceOut_ < numberColumns_)
1642                                   outVub = vub[sequenceOut_];
1643                              int outBack = toVub[outSequence];
1644                              int outSlack = -1;
1645                              if (outSequence >= numberColumns_ && outBack >= 0) {
1646                                   outSlack = outSequence - numberColumns_;
1647                                   outSequence = outBack;
1648                                   outBack = toVub[outSequence];
1649                              }
1650                              if (outVub >= 0)
1651                                   printf("Vub %d out ", outSequence);
1652                              if (outBack >= 0 && outSlack < 0)
1653                                   printf("%d (descendent of %d) out ", outSequence, outBack);
1654                              if (outSlack >= 0)
1655                                   printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
1656                              printf("\n");
1657                         }
1658                    }
1659#endif
1660#if 0
1661                    if (numberIterations_ > 206033)
1662                         handler_->setLogLevel(63);
1663                    if (numberIterations_ > 210567)
1664                         exit(77);
1665#endif
1666                    if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
1667                         handler_->message(CLP_END_VALUES_PASS, messages_)
1668                                   << numberIterations_;
1669                         whatNext = 1;
1670                    }
1671#ifdef CHECK_ACCURACY
1672                    if (whatNext) {
1673                         CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
1674                    }
1675#endif
1676                    //if (numberIterations_==1890)
1677                    //whatNext=1;
1678                    //if (numberIterations_>2000)
1679                    //exit(77);
1680                    // and set bounds correctly
1681                    originalBound(sequenceIn_);
1682                    changeBound(sequenceOut_);
1683#ifdef CLP_DEBUG
1684                    if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1685                         printf("obj backwards %g %g\n", objectiveValue_, oldobj);
1686#endif
1687#if 0
1688                    {
1689                         for (int i = 0; i < numberRows_ + numberColumns_; i++) {
1690                              FakeBound bound = getFakeBound(i);
1691                              if (bound == ClpSimplexDual::upperFake) {
1692                                   assert (upper_[i] < 1.0e20);
1693                              } else if (bound == ClpSimplexDual::lowerFake) {
1694                                   assert (lower_[i] > -1.0e20);
1695                              } else if (bound == ClpSimplexDual::bothFake) {
1696                                   assert (upper_[i] < 1.0e20);
1697                                   assert (lower_[i] > -1.0e20);
1698                              }
1699                         }
1700                    }
1701#endif
1702                    if (whatNext == 1 || candidate == -2) {
1703                         problemStatus_ = -2; // refactorize
1704                    } else if (whatNext == 2) {
1705                         // maximum iterations or equivalent
1706                         problemStatus_ = 3;
1707                         returnCode = 3;
1708                         break;
1709                    }
1710                    // Check event
1711                    {
1712                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1713                         if (status >= 0) {
1714                              problemStatus_ = 5;
1715                              secondaryStatus_ = ClpEventHandler::endOfIteration;
1716                              returnCode = 4;
1717                              break;
1718                         }
1719                    }
1720               } else {
1721#ifdef CLP_INVESTIGATE_SERIAL
1722                    z_thinks = 1;
1723#endif
1724                    // no incoming column is valid
1725                    pivotRow_ = -1;
1726#ifdef CLP_DEBUG
1727                    if (handler_->logLevel() & 32)
1728                         printf("** no column pivot\n");
1729#endif
1730                    if (factorization_->pivots() < 2 && acceptablePivot_ <= 1.0e-8) {
1731                         //&&goodAccuracy()) {
1732                         // If not in branch and bound etc save ray
1733                         delete [] ray_;
1734                         if ((specialOptions_&(1024 | 4096)) == 0 || (specialOptions_ & 32) != 0) {
1735                              // create ray anyway
1736                              ray_ = new double [ numberRows_];
1737                              rowArray_[0]->expand(); // in case packed
1738                              CoinMemcpyN(rowArray_[0]->denseVector(), numberRows_, ray_);
1739                         } else {
1740                              ray_ = NULL;
1741                         }
1742                         // If we have just factorized and infeasibility reasonable say infeas
1743                         double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
1744                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
1745                              double testValue = 1.0e-4;
1746                              if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
1747                                   testValue = 1.0e-6;
1748                              if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue
1749                                        || (specialOptions_ & 64) == 0) {
1750                                   // say infeasible
1751                                   problemStatus_ = 1;
1752                                   // unless primal feasible!!!!
1753                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
1754                                   //   numberDualInfeasibilities_,sumDualInfeasibilities_);
1755                                   //#define TEST_CLP_NODE
1756#ifndef TEST_CLP_NODE
1757                                   // Should be correct - but ...
1758                                   int numberFake = numberAtFakeBound();
1759                                   double sumPrimal =  (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
1760                                   if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 ||
1761                                             (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
1762                                        if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
1763                                             problemStatus_ = -4;
1764#ifdef COIN_DEVELOP
1765                                             printf("status to -4 at %d - primalinf %g pivots %d\n",
1766                                                    __LINE__, sumPrimalInfeasibilities_,
1767                                                    factorization_->pivots());
1768#endif
1769                                        } else {
1770                                             problemStatus_ = 10;
1771#if COIN_DEVELOP>1
1772                                             printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
1773                                                    __LINE__, numberPrimalInfeasibilities_,
1774                                                    sumPrimalInfeasibilities_,
1775                                                    numberDualInfeasibilities_, sumDualInfeasibilities_,
1776                                                    numberFake_, dualBound_, factorization_->pivots(),
1777                                                    (specialOptions_ & 1024) != 0 ? 1 : 0,
1778                                                    (specialOptions_ & 2048) != 0 ? 1 : 0,
1779                                                    (specialOptions_ & 4096) != 0 ? 1 : 0,
1780                                                    (specialOptions_ & 8192) != 0 ? 1 : 0,
1781                                                    (specialOptions_ & 16384) != 0 ? 1 : 0
1782                                                   );
1783#endif
1784                                             // Get rid of objective
1785                                             if ((specialOptions_ & 16384) == 0)
1786                                                  objective_ = new ClpLinearObjective(NULL, numberColumns_);
1787                                        }
1788                                   }
1789#else
1790                                   if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) {
1791#ifdef COIN_DEVELOP
1792                                        printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1793                                               __LINE__, numberPrimalInfeasibilities_,
1794                                               sumPrimalInfeasibilities_,
1795                                               numberDualInfeasibilities_, sumDualInfeasibilities_,
1796                                               numberFake_, dualBound_, factorization_->pivots());
1797#endif
1798                                        if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
1799                                             problemStatus_ = 10;
1800#if COIN_DEVELOP>1
1801                                             printf("returning at %d\n", __LINE__);
1802#endif
1803                                             // Get rid of objective
1804                                             if ((specialOptions_ & 16384) == 0)
1805                                                  objective_ = new ClpLinearObjective(NULL, numberColumns_);
1806                                        }
1807                                   }
1808#endif
1809                                   rowArray_[0]->clear();
1810                                   columnArray_[0]->clear();
1811                                   returnCode = 1;
1812                                   break;
1813                              }
1814                         }
1815                         // If special option set - put off as long as possible
1816                         if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) {
1817                              if (factorization_->pivots() == 0)
1818                                   problemStatus_ = -4; //say looks infeasible
1819                         } else {
1820                              // flag
1821                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
1822                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
1823                                        << x << sequenceWithin(sequenceOut_)
1824                                        << CoinMessageEol;
1825#ifdef COIN_DEVELOP
1826                              printf("flag c\n");
1827#endif
1828                              setFlagged(sequenceOut_);
1829                              if (!factorization_->pivots()) {
1830                                   rowArray_[0]->clear();
1831                                   columnArray_[0]->clear();
1832                                   continue;
1833                              }
1834                         }
1835                    }
1836                    if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
1837                         acceptablePivot_ = 1.0e-8;
1838                    rowArray_[0]->clear();
1839                    columnArray_[0]->clear();
1840                    returnCode = 1;
1841                    break;
1842               }
1843          } else {
1844#ifdef CLP_INVESTIGATE_SERIAL
1845               z_thinks = 0;
1846#endif
1847               // no pivot row
1848#ifdef CLP_DEBUG
1849               if (handler_->logLevel() & 32)
1850                    printf("** no row pivot\n");
1851#endif
1852               // If in branch and bound try and get rid of fixed variables
1853               if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
1854                    assert (!candidateList);
1855                    candidateList = new int[numberRows_];
1856                    int iRow;
1857                    for (iRow = 0; iRow < numberRows_; iRow++) {
1858                         int iPivot = pivotVariable_[iRow];
1859                         if (flagged(iPivot) || !pivoted(iPivot))
1860                              continue;
1861                         assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
1862                         candidateList[numberCandidates++] = iRow;
1863                    }
1864                    // and set first candidate
1865                    if (!numberCandidates) {
1866                         delete [] candidateList;
1867                         candidateList = NULL;
1868                         int iRow;
1869                         for (iRow = 0; iRow < numberRows_; iRow++) {
1870                              int iPivot = pivotVariable_[iRow];
1871                              clearPivoted(iPivot);
1872                         }
1873                    } else {
1874                         ifValuesPass = 2;
1875                         continue;
1876                    }
1877               }
1878               int numberPivots = factorization_->pivots();
1879               bool specialCase;
1880               int useNumberFake;
1881               returnCode = 0;
1882               if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
1883                         (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101)
1884                         && dualBound_ >= 1.0e8) {
1885                    specialCase = true;
1886                    // as dual bound high - should be okay
1887                    useNumberFake = 0;
1888               } else {
1889                    specialCase = false;
1890                    useNumberFake = numberFake_;
1891               }
1892               if (!numberPivots || specialCase) {
1893                    // may have crept through - so may be optimal
1894                    // check any flagged variables
1895                    int iRow;
1896                    for (iRow = 0; iRow < numberRows_; iRow++) {
1897                         int iPivot = pivotVariable_[iRow];
1898                         if (flagged(iPivot))
1899                              break;
1900                    }
1901                    if (iRow < numberRows_ && numberPivots) {
1902                         // try factorization
1903                         returnCode = -2;
1904                    }
1905
1906                    if (useNumberFake || numberDualInfeasibilities_) {
1907                         // may be dual infeasible
1908                         if ((specialOptions_ & 1024) == 0)
1909                              problemStatus_ = -5;
1910                         else if (!useNumberFake && numberPrimalInfeasibilities_
1911                                   && !numberPivots)
1912                              problemStatus_ = 1;
1913                    } else {
1914                         if (iRow < numberRows_) {
1915#ifdef COIN_DEVELOP
1916                              std::cout << "Flagged variables at end - infeasible?" << std::endl;
1917                              printf("Probably infeasible - pivot was %g\n", alpha_);
1918#endif
1919                              //if (fabs(alpha_)<1.0e-4) {
1920                              //problemStatus_=1;
1921                              //} else {
1922#ifdef CLP_DEBUG
1923                              abort();
1924#endif
1925                              //}
1926                              problemStatus_ = -5;
1927                         } else {
1928                              problemStatus_ = 0;
1929#ifndef CLP_CHECK_NUMBER_PIVOTS
1930#define CLP_CHECK_NUMBER_PIVOTS 10
1931#endif
1932#if CLP_CHECK_NUMBER_PIVOTS < 20
1933                              if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
1934#ifndef NDEBUG_CLP
1935                                   int nTotal = numberRows_ + numberColumns_;
1936                                   double * comp = CoinCopyOfArray(solution_, nTotal);
1937#endif
1938                                   computePrimals(rowActivityWork_, columnActivityWork_);
1939#ifndef NDEBUG_CLP
1940                                   double largest = 1.0e-5;
1942                                   for (int i = 0; i < nTotal; i++) {
1943                                        double value = solution_[i];
1944                                        double larger = CoinMax(fabs(value), fabs(comp[i]));
1945                                        double tol = 1.0e-5 + 1.0e-5 * larger;
1946                                        double diff = fabs(value - comp[i]);
1947                                        if (diff - tol > largest) {
1949                                             largest = diff - tol;
1950                                        }
1951                                   }
1954#endif
1955                                   checkPrimalSolution(rowActivityWork_, columnActivityWork_);
1956                                   if (numberPrimalInfeasibilities_) {
1957#ifdef CLP_INVESTIGATE
1958                                        printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
1959                                               sumPrimalInfeasibilities_);
1960#endif
1961                                        problemStatus_ = -1;
1962                                        returnCode = -2;
1963                                   }
1964#ifndef NDEBUG_CLP
1965                                   memcpy(solution_, comp, nTotal * sizeof(double));
1966                                   delete [] comp;
1967#endif
1968                              }
1969#endif
1970                              if (!problemStatus_) {
1971                                   // make it look OK
1972                                   numberPrimalInfeasibilities_ = 0;
1973                                   sumPrimalInfeasibilities_ = 0.0;
1974                                   numberDualInfeasibilities_ = 0;
1975                                   sumDualInfeasibilities_ = 0.0;
1976                                   // May be perturbed
1977                                   if (perturbation_ == 101 || numberChanged_) {
1978                                        numberChanged_ = 0; // Number of variables with changed costs
1979                                        perturbation_ = 102; // stop any perturbations
1980                                        //double changeCost;
1981                                        //changeBounds(1,NULL,changeCost);
1982                                        createRim4(false);
1983                                        // make sure duals are current
1984                                        computeDuals(givenDuals);
1985                                        checkDualSolution();
1986                                        progress_.modifyObjective(-COIN_DBL_MAX);
1987                                        if (numberDualInfeasibilities_) {
1988                                             problemStatus_ = 10; // was -3;
1989                                        } else {
1990                                             computeObjectiveValue(true);
1991                                        }
1992                                   } else if (numberPivots) {
1993                                        computeObjectiveValue(true);
1994                                   }
1995                                   if (numberPivots < -1000) {
1996                                        // objective may be wrong
1997                                        objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
1998                                        objectiveValue_ += objective_->nonlinearOffset();
1999                                        objectiveValue_ /= (objectiveScale_ * rhsScale_);
2000                                        if ((specialOptions_ & 16384) == 0) {
2001                                             // and dual_ may be wrong (i.e. for fixed or basic)
2002                                             CoinIndexedVector * arrayVector = rowArray_[1];
2003                                             arrayVector->clear();
2004                                             int iRow;
2005                                             double * array = arrayVector->denseVector();
2006                                             /* Use dual_ instead of array
2007                                                Even though dual_ is only numberRows_ long this is
2008                                                okay as gets permuted to longer rowArray_[2]
2009                                             */
2010                                             arrayVector->setDenseVector(dual_);
2011                                             int * index = arrayVector->getIndices();
2012                                             int number = 0;
2013                                             for (iRow = 0; iRow < numberRows_; iRow++) {
2014                                                  int iPivot = pivotVariable_[iRow];
2015                                                  double value = cost_[iPivot];
2016                                                  dual_[iRow] = value;
2017                                                  if (value) {
2018                                                       index[number++] = iRow;
2019                                                  }
2020                                             }
2021                                             arrayVector->setNumElements(number);
2022                                             // Extended duals before "updateTranspose"
2023                                             matrix_->dualExpanded(this, arrayVector, NULL, 0);
2024                                             // Btran basic costs
2025                                             rowArray_[2]->clear();
2026                                             factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2027                                             // and return vector
2028                                             arrayVector->setDenseVector(array);
2029                                        }
2030                                   }
2031                                   sumPrimalInfeasibilities_ = 0.0;
2032                              }
2033                              if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
2034                                   CoinIndexedVector * arrayVector = rowArray_[1];
2035                                   arrayVector->clear();
2036                                   double * rhs = arrayVector->denseVector();
2037                                   times(1.0, solution_, rhs);
2038#ifdef CHECK_ACCURACY
2040#endif
2042                                   int i;
2043                                   for ( i = 0; i < numberRows_; i++) {
2044                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
2045                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2047#ifdef CHECK_ACCURACY
2048                                             printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2049                                                    rowLowerWork_[i], rowUpperWork_[i],
2050                                                    rhs[i], rowActivityWork_[i]);
2051#endif
2052                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2053#ifdef CHECK_ACCURACY
2055                                             printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2056#endif
2057                                        }
2058                                        rhs[i] = 0.0;
2059                                   }
2060                                   for ( i = 0; i < numberColumns_; i++) {
2061                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
2062                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2064#ifdef CHECK_ACCURACY
2065                                             printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2066                                                    columnLowerWork_[i], columnUpperWork_[i],
2067                                                    solution_[i], columnActivityWork_[i]);
2068#endif
2069                                        }
2070                                   }
2072                                        problemStatus_ = -3;
2073                                        returnCode = -2;
2074                                        // Force to re-factorize early next time
2075                                        int numberPivots = factorization_->pivots();
2076                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2077                                   }
2078                              }
2079                         }
2080                    }
2081               } else {
2082                    problemStatus_ = -3;
2083                    returnCode = -2;
2084                    // Force to re-factorize early next time
2085                    int numberPivots = factorization_->pivots();
2086                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2087               }
2088               break;
2089          }
2090     }
2091     if (givenDuals) {
2092          CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2093          // get rid of any values pass array
2094          delete [] candidateList;
2095     }
2096     delete [] dubiousWeights;
2097#ifdef CLP_REPORT_PROGRESS
2098     if (ixxxxxx > ixxyyyy - 5) {
2099          int nTotal = numberColumns_ + numberRows_;
2100          double oldObj = 0.0;
2101          double newObj = 0.0;
2102          for (int i = 0; i < nTotal; i++) {
2103               if (savePSol[i])
2104                    oldObj += savePSol[i] * saveCost[i];
2105               if (solution_[i])
2106                    newObj += solution_[i] * cost_[i];
2107               bool printIt = false;
2108               if (cost_[i] != saveCost[i])
2109                    printIt = true;
2110               if (status_[i] != saveStat[i])
2111                    printIt = true;
2112               if (printIt)
2113                    printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2114                           i, saveStat[i], saveCost[i], savePSol[i],
2115                           status_[i], cost_[i], solution_[i]);
2116               // difference
2117               savePSol[i] = solution_[i] - savePSol[i];
2118          }
2119          printf("exit pivots %d, old obj %g new %g\n",
2120                 factorization_->pivots(),
2121                 oldObj, newObj);
2122          memset(saveDj, 0, numberRows_ * sizeof(double));
2123          times(1.0, savePSol, saveDj);
2124          double largest = 1.0e-6;
2125          int k = -1;
2126          for (int i = 0; i < numberRows_; i++) {
2127               saveDj[i] -= savePSol[i+numberColumns_];
2128               if (fabs(saveDj[i]) > largest) {
2129                    largest = fabs(saveDj[i]);
2130                    k = i;
2131               }
2132          }
2133          if (k >= 0)
2134               printf("Not null %d %g\n", k, largest);
2135     }
2136     delete [] savePSol ;
2137     delete [] saveDj ;
2138     delete [] saveCost ;
2139     delete [] saveStat ;
2140#endif
2141     return returnCode;
2142}
2143/* The duals are updated by the given arrays.
2144   Returns number of infeasibilities.
2145   rowArray and columnarray will have flipped
2146   The output vector has movement (row length array) */
2147int
2148ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
2149                                  CoinIndexedVector * columnArray,
2150                                  CoinIndexedVector * outputArray,
2151                                  double theta,
2152                                  double & objectiveChange,
2153                                  bool fullRecompute)
2154{
2155
2156     outputArray->clear();
2157
2158
2159     int numberInfeasibilities = 0;
2160     int numberRowInfeasibilities = 0;
2161
2162     // get a tolerance
2163     double tolerance = dualTolerance_;
2164     // we can't really trust infeasibilities if there is dual error
2165     double error = CoinMin(1.0e-2, largestDualError_);
2166     // allow tolerance at least slightly bigger than standard
2167     tolerance = tolerance +  error;
2168
2169     double changeObj = 0.0;
2170
2171     // Coding is very similar but we can save a bit by splitting
2172     // Do rows
2173     if (!fullRecompute) {
2174          int i;
2175          double * COIN_RESTRICT reducedCost = djRegion(0);
2176          const double * COIN_RESTRICT lower = lowerRegion(0);
2177          const double * COIN_RESTRICT upper = upperRegion(0);
2178          const double * COIN_RESTRICT cost = costRegion(0);
2179          double * COIN_RESTRICT work;
2180          int number;
2181          int * COIN_RESTRICT which;
2182          const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
2183          assert(rowArray->packedMode());
2184          work = rowArray->denseVector();
2185          number = rowArray->getNumElements();
2186          which = rowArray->getIndices();
2187          double multiplier[] = { -1.0, 1.0};
2188          for (i = 0; i < number; i++) {
2189               int iSequence = which[i];
2190               double alphaI = work[i];
2191               work[i] = 0.0;
2192               int iStatus = (statusArray[iSequence] & 3) - 1;
2193               if (iStatus) {
2194                    double value = reducedCost[iSequence] - theta * alphaI;
2195                    reducedCost[iSequence] = value;
2196                    double mult = multiplier[iStatus-1];
2197                    value *= mult;
2198                    if (value < -tolerance) {
2199                         // flipping bounds
2200                         double movement = mult * (lower[iSequence] - upper[iSequence]);
2201                         which[numberInfeasibilities++] = iSequence;
2202#ifndef NDEBUG
2203                         if (fabs(movement) >= 1.0e30)
2204                              resetFakeBounds(-1000 - iSequence);
2205#endif
2206#ifdef CLP_DEBUG
2207                         if ((handler_->logLevel() & 32))
2208                              printf("%d %d, new dj %g, alpha %g, movement %g\n",
2209                                     0, iSequence, value, alphaI, movement);
2210#endif
2211                         changeObj -= movement * cost[iSequence];
2213                    }
2214               }
2215          }
2216          // Do columns
2217          reducedCost = djRegion(1);
2218          lower = lowerRegion(1);
2219          upper = upperRegion(1);
2220          cost = costRegion(1);
2221          // set number of infeasibilities in row array
2222          numberRowInfeasibilities = numberInfeasibilities;
2223          rowArray->setNumElements(numberInfeasibilities);
2224          numberInfeasibilities = 0;
2225          work = columnArray->denseVector();
2226          number = columnArray->getNumElements();
2227          which = columnArray->getIndices();
2228          if ((moreSpecialOptions_ & 8) != 0) {
2229               const unsigned char * COIN_RESTRICT statusArray = status_;
2230               for (i = 0; i < number; i++) {
2231                    int iSequence = which[i];
2232                    double alphaI = work[i];
2233                    work[i] = 0.0;
2234
2235                    int iStatus = (statusArray[iSequence] & 3) - 1;
2236                    if (iStatus) {
2237                         double value = reducedCost[iSequence] - theta * alphaI;
2238                         reducedCost[iSequence] = value;
2239                         double mult = multiplier[iStatus-1];
2240                         value *= mult;
2241                         if (value < -tolerance) {
2242                              // flipping bounds
2243                              double movement = mult * (upper[iSequence] - lower[iSequence]);
2244                              which[numberInfeasibilities++] = iSequence;
2245#ifndef NDEBUG
2246                              if (fabs(movement) >= 1.0e30)
2247                                   resetFakeBounds(-1000 - iSequence);
2248#endif
2249#ifdef CLP_DEBUG
2250                              if ((handler_->logLevel() & 32))
2251                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2252                                          1, iSequence, value, alphaI, movement);
2253#endif
2254                              changeObj += movement * cost[iSequence];
2256                         }
2257                    }
2258               }
2259          } else {
2260               for (i = 0; i < number; i++) {
2261                    int iSequence = which[i];
2262                    double alphaI = work[i];
2263                    work[i] = 0.0;
2264
2265                    Status status = getStatus(iSequence);
2266                    if (status == atLowerBound) {
2267                         double value = reducedCost[iSequence] - theta * alphaI;
2268                         reducedCost[iSequence] = value;
2269                         double movement = 0.0;
2270
2271                         if (value < -tolerance) {
2272                              // to upper bound
2273                              which[numberInfeasibilities++] = iSequence;
2274                              movement = upper[iSequence] - lower[iSequence];
2275#ifndef NDEBUG
2276                              if (fabs(movement) >= 1.0e30)
2277                                   resetFakeBounds(-1000 - iSequence);
2278#endif
2279#ifdef CLP_DEBUG
2280                              if ((handler_->logLevel() & 32))
2281                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2282                                          1, iSequence, value, alphaI, movement);
2283#endif
2284                              changeObj += movement * cost[iSequence];
2286                         }
2287                    } else if (status == atUpperBound) {
2288                         double value = reducedCost[iSequence] - theta * alphaI;
2289                         reducedCost[iSequence] = value;
2290                         double movement = 0.0;
2291
2292                         if (value > tolerance) {
2293                              // to lower bound (if swap)
2294                              which[numberInfeasibilities++] = iSequence;
2295                              movement = lower[iSequence] - upper[iSequence];
2296#ifndef NDEBUG
2297                              if (fabs(movement) >= 1.0e30)
2298                                   resetFakeBounds(-1000 - iSequence);
2299#endif
2300#ifdef CLP_DEBUG
2301                              if ((handler_->logLevel() & 32))
2302                                   printf("%d %d, new dj %g, alpha %g, movement %g\n",
2303                                          1, iSequence, value, alphaI, movement);
2304#endif
2305                              changeObj += movement * cost[iSequence];
2307                         }
2308                    } else if (status == isFree) {
2309                         double value = reducedCost[iSequence] - theta * alphaI;
2310                         reducedCost[iSequence] = value;
2311                    }
2312               }
2313          }
2314     } else  {
2315          double * COIN_RESTRICT solution = solutionRegion(0);
2316          double * COIN_RESTRICT reducedCost = djRegion(0);
2317          const double * COIN_RESTRICT lower = lowerRegion(0);
2318          const double * COIN_RESTRICT upper = upperRegion(0);
2319          const double * COIN_RESTRICT cost = costRegion(0);
2320          int * COIN_RESTRICT which;
2321          which = rowArray->getIndices();
2322          int iSequence;
2323          for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2324               double value = reducedCost[iSequence];
2325
2326               Status status = getStatus(iSequence + numberColumns_);
2327               // more likely to be at upper bound ?
2328               if (status == atUpperBound) {
2329                    double movement = 0.0;
2330                    //#define NO_SWAP7
2331                    if (value > tolerance) {
2332                         // to lower bound (if swap)
2333                         // put back alpha
2334                         which[numberInfeasibilities++] = iSequence;
2335                         movement = lower[iSequence] - upper[iSequence];
2336                         changeObj += movement * cost[iSequence];
2338#ifndef NDEBUG
2339                         if (fabs(movement) >= 1.0e30)
2340                              resetFakeBounds(-1000 - iSequence);
2341#endif
2342#ifndef NO_SWAP7
2343                    } else if (value > -tolerance) {
2344                         // at correct bound but may swap
2345                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2346                         if (bound == ClpSimplexDual::upperFake) {
2347                              movement = lower[iSequence] - upper[iSequence];
2348#ifndef NDEBUG
2349                              if (fabs(movement) >= 1.0e30)
2350                                   resetFakeBounds(-1000 - iSequence);
2351#endif
2352                              setStatus(iSequence + numberColumns_, atLowerBound);
2353                              solution[iSequence] = lower[iSequence];
2354                              changeObj += movement * cost[iSequence];
2355                              //numberFake_--;
2356                              //setFakeBound(iSequence+numberColumns_,noFake);
2357                         }
2358#endif
2359                    }
2360               } else if (status == atLowerBound) {
2361                    double movement = 0.0;
2362
2363                    if (value < -tolerance) {
2364                         // to upper bound
2365                         // put back alpha
2366                         which[numberInfeasibilities++] = iSequence;
2367                         movement = upper[iSequence] - lower[iSequence];
2368#ifndef NDEBUG
2369                         if (fabs(movement) >= 1.0e30)
2370                              resetFakeBounds(-1000 - iSequence);
2371#endif
2372                         changeObj += movement * cost[iSequence];
2374#ifndef NO_SWAP7
2375                    } else if (value < tolerance) {
2376                         // at correct bound but may swap
2377                         FakeBound bound = getFakeBound(iSequence + numberColumns_);
2378                         if (bound == ClpSimplexDual::lowerFake) {
2379                              movement = upper[iSequence] - lower[iSequence];
2380#ifndef NDEBUG
2381                              if (fabs(movement) >= 1.0e30)
2382                                   resetFakeBounds(-1000 - iSequence);
2383#endif
2384                              setStatus(iSequence + numberColumns_, atUpperBound);
2385                              solution[iSequence] = upper[iSequence];
2386                              changeObj += movement * cost[iSequence];
2387                              //numberFake_--;
2388                              //setFakeBound(iSequence+numberColumns_,noFake);
2389                         }
2390#endif
2391                    }
2392               }
2393          }
2394          // Do columns
2395          solution = solutionRegion(1);
2396          reducedCost = djRegion(1);
2397          lower = lowerRegion(1);
2398          upper = upperRegion(1);
2399          cost = costRegion(1);
2400          // set number of infeasibilities in row array
2401          numberRowInfeasibilities = numberInfeasibilities;
2402          rowArray->setNumElements(numberInfeasibilities);
2403          numberInfeasibilities = 0;
2404          which = columnArray->getIndices();
2405          for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2406               double value = reducedCost[iSequence];
2407
2408               Status status = getStatus(iSequence);
2409               if (status == atLowerBound) {
2410                    double movement = 0.0;
2411
2412                    if (value < -tolerance) {
2413                         // to upper bound
2414                         // put back alpha
2415                         which[numberInfeasibilities++] = iSequence;
2416                         movement = upper[iSequence] - lower[iSequence];
2417#ifndef NDEBUG
2418                         if (fabs(movement) >= 1.0e30)
2419                              resetFakeBounds(-1000 - iSequence);
2420#endif
2421                         changeObj += movement * cost[iSequence];
2423#ifndef NO_SWAP7
2424                    } else if (value < tolerance) {
2425                         // at correct bound but may swap
2426                         FakeBound bound = getFakeBound(iSequence);
2427                         if (bound == ClpSimplexDual::lowerFake) {
2428                              movement = upper[iSequence] - lower[iSequence];
2429#ifndef NDEBUG
2430                              if (fabs(movement) >= 1.0e30)
2431                                   resetFakeBounds(-1000 - iSequence);
2432#endif
2433                              setStatus(iSequence, atUpperBound);
2434                              solution[iSequence] = upper[iSequence];
2435                              changeObj += movement * cost[iSequence];
2436                              //numberFake_--;
2437                              //setFakeBound(iSequence,noFake);
2438                         }
2439#endif
2440                    }
2441               } else if (status == atUpperBound) {
2442                    double movement = 0.0;
2443
2444                    if (value > tolerance) {
2445                         // to lower bound (if swap)
2446                         // put back alpha
2447                         which[numberInfeasibilities++] = iSequence;
2448                         movement = lower[iSequence] - upper[iSequence];
2449#ifndef NDEBUG
2450                         if (fabs(movement) >= 1.0e30)
2451                              resetFakeBounds(-1000 - iSequence);
2452#endif
2453                         changeObj += movement * cost[iSequence];
2455#ifndef NO_SWAP7
2456                    } else if (value > -tolerance) {
2457                         // at correct bound but may swap
2458                         FakeBound bound = getFakeBound(iSequence);
2459                         if (bound == ClpSimplexDual::upperFake) {
2460                              movement = lower[iSequence] - upper[iSequence];
2461#ifndef NDEBUG
2462                              if (fabs(movement) >= 1.0e30)
2463                                   resetFakeBounds(-1000 - iSequence);
2464#endif
2465                              setStatus(iSequence, atLowerBound);
2466                              solution[iSequence] = lower[iSequence];
2467                              changeObj += movement * cost[iSequence];
2468                              //numberFake_--;
2469                              //setFakeBound(iSequence,noFake);
2470                         }
2471#endif
2472                    }
2473               }
2474          }
2475     }
2476
2477#ifdef CLP_DEBUG
2478     if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2479          printf("%d fake after full update\n", numberFake_);
2480#endif
2481     // set number of infeasibilities
2482     columnArray->setNumElements(numberInfeasibilities);
2483     numberInfeasibilities += numberRowInfeasibilities;
2484     if (fullRecompute) {
2485          // do actual flips
2486          flipBounds(rowArray, columnArray);
2487     }
2488     objectiveChange += changeObj;
2489     return numberInfeasibilities;
2490}
2491void
2492ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2493                                        CoinIndexedVector * columnArray,
2494                                        double theta)
2495{
2496
2497     // use a tighter tolerance except for all being okay
2498     double tolerance = dualTolerance_;
2499
2500     // Coding is very similar but we can save a bit by splitting
2501     // Do rows
2502     {
2503          int i;
2504          double * reducedCost = djRegion(0);
2505          double * work;
2506          int number;
2507          int * which;
2508          work = rowArray->denseVector();
2509          number = rowArray->getNumElements();
2510          which = rowArray->getIndices();
2511          for (i = 0; i < number; i++) {
2512               int iSequence = which[i];
2513               double alphaI = work[i];
2514               double value = reducedCost[iSequence] - theta * alphaI;
2515               work[i] = 0.0;
2516               reducedCost[iSequence] = value;
2517
2518               Status status = getStatus(iSequence + numberColumns_);
2519               // more likely to be at upper bound ?
2520               if (status == atUpperBound) {
2521
2522                    if (value > tolerance)
2523                         reducedCost[iSequence] = 0.0;
2524               } else if (status == atLowerBound) {
2525
2526                    if (value < -tolerance) {
2527                         reducedCost[iSequence] = 0.0;
2528                    }
2529               }
2530          }
2531     }
2532     rowArray->setNumElements(0);
2533
2534     // Do columns
2535     {
2536          int i;
2537          double * reducedCost = djRegion(1);
2538          double * work;
2539          int number;
2540          int * which;
2541          work = columnArray->denseVector();
2542          number = columnArray->getNumElements();
2543          which = columnArray->getIndices();
2544
2545          for (i = 0; i < number; i++) {
2546               int iSequence = which[i];
2547               double alphaI = work[i];
2548               double value = reducedCost[iSequence] - theta * alphaI;
2549               work[i] = 0.0;
2550               reducedCost[iSequence] = value;
2551
2552               Status status = getStatus(iSequence);
2553               if (status == atLowerBound) {
2554                    if (value < -tolerance)
2555                         reducedCost[iSequence] = 0.0;
2556               } else if (status == atUpperBound) {
2557                    if (value > tolerance)
2558                         reducedCost[iSequence] = 0.0;
2559               }
2560          }
2561     }
2562     columnArray->setNumElements(0);
2563}
2564/*
2565   Chooses dual pivot row
2566   Would be faster with separate region to scan
2567   and will have this (with square of infeasibility) when steepest
2568   For easy problems we can just choose one of the first rows we look at
2569*/
2570void
2572{
2573     // get pivot row using whichever method it is
2574     int chosenRow = -1;
2575#ifdef FORCE_FOLLOW
2576     bool forceThis = false;
2577     if (!fpFollow && strlen(forceFile)) {
2579          assert (fpFollow);
2580     }
2581     if (fpFollow) {
2582          if (numberIterations_ <= force_iteration) {
2583               // read to next Clp0102
2584               char temp[300];
2585               while (fgets(temp, 250, fpFollow)) {
2586                    if (strncmp(temp, "Clp0102", 7))
2587                         continue;
2588                    char cin, cout;
2589                    sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
2590                           &force_iteration, &cin, &force_in, &cout, &force_out);
2591                    if (cin == 'R')
2592                         force_in += numberColumns_;
2593                    if (cout == 'R')
2594                         force_out += numberColumns_;
2595                    forceThis = true;
2596                    assert (numberIterations_ == force_iteration - 1);
2597                    printf("Iteration %d will force %d out and %d in\n",
2598                           force_iteration, force_out, force_in);
2600                    break;
2601               }
2602          } else {
2603               // use old
2604               forceThis = true;
2605          }
2606          if (!forceThis) {
2607               fclose(fpFollow);
2609               forceFile = "";
2610          }
2611     }
2612#endif
2613     double freeAlpha = 0.0;
2614     if (alreadyChosen < 0) {
2615          // first see if any free variables and put them in basis
2616          int nextFree = nextSuperBasic();
2617          //nextFree=-1; //off
2618          if (nextFree >= 0) {
2619               // unpack vector and find a good pivot
2620               unpack(rowArray_[1], nextFree);
2621               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2622
2623               double * work = rowArray_[1]->denseVector();
2624               int number = rowArray_[1]->getNumElements();
2625               int * which = rowArray_[1]->getIndices();
2626               double bestFeasibleAlpha = 0.0;
2627               int bestFeasibleRow = -1;
2628               double bestInfeasibleAlpha = 0.0;
2629               int bestInfeasibleRow = -1;
2630               int i;
2631
2632               for (i = 0; i < number; i++) {
2633                    int iRow = which[i];
2634                    double alpha = fabs(work[iRow]);
2635                    if (alpha > 1.0e-3) {
2636                         int iSequence = pivotVariable_[iRow];
2637                         double value = solution_[iSequence];
2638                         double lower = lower_[iSequence];
2639                         double upper = upper_[iSequence];
2640                         double infeasibility = 0.0;
2641                         if (value > upper)
2642                              infeasibility = value - upper;
2643                         else if (value < lower)
2644                              infeasibility = lower - value;
2645                         if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2646                              if (!flagged(iSequence)) {
2647                                   bestInfeasibleAlpha = infeasibility * alpha;
2648                                   bestInfeasibleRow = iRow;
2649                              }
2650                         }
2651                         if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
2652                              bestFeasibleAlpha = alpha;
2653                              bestFeasibleRow = iRow;
2654                         }
2655                    }
2656               }
2657               if (bestInfeasibleRow >= 0)
2658                    chosenRow = bestInfeasibleRow;
2659               else if (bestFeasibleAlpha > 1.0e-2)
2660                    chosenRow = bestFeasibleRow;
2661               if (chosenRow >= 0) {
2662                    pivotRow_ = chosenRow;
2663                    freeAlpha = work[chosenRow];
2664               }
2665               rowArray_[1]->clear();
2666          }
2667     } else {
2668          // in values pass
2670#ifdef FORCE_FOLLOW
2671          if(forceThis) {
2673               chosenRow = -1;
2674               for (int i = 0; i < numberRows_; i++) {
2675                    if (pivotVariable_[i] == force_out) {
2676                         chosenRow = i;
2677                         break;
2678                    }
2679               }
2680               assert (chosenRow >= 0);
2681          }
2682#endif
2683          pivotRow_ = chosenRow;
2684     }
2685     if (chosenRow < 0)
2686          pivotRow_ = dualRowPivot_->pivotRow();
2687
2688     if (pivotRow_ >= 0) {
2689          sequenceOut_ = pivotVariable_[pivotRow_];
2690          valueOut_ = solution_[sequenceOut_];
2691          lowerOut_ = lower_[sequenceOut_];
2692          upperOut_ = upper_[sequenceOut_];
2693          if (alreadyChosen < 0) {
2694               // if we have problems we could try other way and hope we get a
2695               // zero pivot?
2696               if (valueOut_ > upperOut_) {
2697                    directionOut_ = -1;
2698                    dualOut_ = valueOut_ - upperOut_;
2699               } else if (valueOut_ < lowerOut_) {
2700                    directionOut_ = 1;
2701                    dualOut_ = lowerOut_ - valueOut_;
2702               } else {
2703#if 1
2704                    // odd (could be free) - it's feasible - go to nearest
2705                    if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2706                         directionOut_ = 1;
2707                         dualOut_ = lowerOut_ - valueOut_;
2708                    } else {
2709                         directionOut_ = -1;
2710                         dualOut_ = valueOut_ - upperOut_;
2711                    }
2712#else
2713                    // odd (could be free) - it's feasible - improve obj
2714                    printf("direction from alpha of %g is %d\n",
2715                           freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2716                    if (valueOut_ - lowerOut_ > 1.0e20)
2717                         freeAlpha = 1.0;
2718                    else if(upperOut_ - valueOut_ > 1.0e20)
2719                         freeAlpha = -1.0;
2720                    //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2721                    if (freeAlpha < 0.0) {
2722                         directionOut_ = 1;
2723                         dualOut_ = lowerOut_ - valueOut_;
2724                    } else {
2725                         directionOut_ = -1;
2726                         dualOut_ = valueOut_ - upperOut_;
2727                    }
2728                    printf("direction taken %d - bounds %g %g %g\n",
2729                           directionOut_, lowerOut_, valueOut_, upperOut_);
2730#endif
2731               }
2732#ifdef CLP_DEBUG
2733               assert(dualOut_ >= 0.0);
2734#endif
2735          } else {
2736               // in values pass so just use sign of dj
2737               // We don't want to go through any barriers so set dualOut low
2738               // free variables will never be here
2739               dualOut_ = 1.0e-6;
2740               if (dj_[sequenceOut_] > 0.0) {
2741                    // this will give a -1 in pivot row (as slacks are -1.0)
2742                    directionOut_ = 1;
2743               } else {
2744                    directionOut_ = -1;
2745               }
2746          }
2747     }
2748     return ;
2749}
2750// Checks if any fake bounds active - if so returns number and modifies
2751// dualBound_ and everything.
2752// Free variables will be left as free
2753// Returns number of bounds changed if >=0
2754// Returns -1 if not initialize and no effect
2755// Fills in changeVector which can be used to see if unbounded
2756// and cost of change vector
2757int
2758ClpSimplexDual::changeBounds(int initialize,
2759                             CoinIndexedVector * outputArray,
2760                             double & changeCost)
2761{
2762     numberFake_ = 0;
2763     if (!initialize) {
2764          int numberInfeasibilities;
2765          double newBound;
2766          newBound = 5.0 * dualBound_;
2767          numberInfeasibilities = 0;
2768          changeCost = 0.0;
2769          // put back original bounds and then check
2770          createRim1(false);
2771          int iSequence;
2772          // bounds will get bigger - just look at ones at bounds
2773          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2774               double lowerValue = lower_[iSequence];
2775               double upperValue = upper_[iSequence];
2776               double value = solution_[iSequence];
2777               setFakeBound(iSequence, ClpSimplexDual::noFake);
2778               switch(getStatus(iSequence)) {
2779
2780               case basic:
2781               case ClpSimplex::isFixed:
2782                    break;
2783               case isFree:
2784               case superBasic:
2785                    break;
2786               case atUpperBound:
2787                    if (fabs(value - upperValue) > primalTolerance_)
2788                         numberInfeasibilities++;
2789                    break;
2790               case atLowerBound:
2791                    if (fabs(value - lowerValue) > primalTolerance_)
2792                         numberInfeasibilities++;
2793                    break;
2794               }
2795          }
2796          // If dual infeasible then carry on
2797          if (numberInfeasibilities) {
2798               handler_->message(CLP_DUAL_CHECKB, messages_)
2799                         << newBound
2800                         << CoinMessageEol;
2801               int iSequence;
2802               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2803                    double lowerValue = lower_[iSequence];
2804                    double upperValue = upper_[iSequence];
2805                    double newLowerValue;
2806                    double newUpperValue;
2807                    Status status = getStatus(iSequence);
2808                    if (status == atUpperBound ||
2809                              status == atLowerBound) {
2810                         double value = solution_[iSequence];
2811                         if (value - lowerValue <= upperValue - value) {
2812                              newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
2813                              newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
2814                         } else {
2815                              newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
2816                              newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
2817                         }
2818                         lower_[iSequence] = newLowerValue;
2819                         upper_[iSequence] = newUpperValue;
2820                         if (newLowerValue > lowerValue) {
2821                              if (newUpperValue < upperValue) {
2822                                   setFakeBound(iSequence, ClpSimplexDual::bothFake);
2823#ifdef CLP_INVESTIGATE
2824                                   abort(); // No idea what should happen here - I have never got here
2825#endif
2826                                   numberFake_++;
2827                              } else {
2828                                   setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2829                                   numberFake_++;
2830                              }
2831                         } else {
2832                              if (newUpperValue < upperValue) {
2833                                   setFakeBound(iSequence, ClpSimplexDual::upperFake);
2834                                   numberFake_++;
2835                              }
2836                         }
2837                         if (status == atUpperBound)
2838                              solution_[iSequence] = newUpperValue;
2839                         else
2840                              solution_[iSequence] = newLowerValue;
2841                         double movement = solution_[iSequence] - value;
2842                         if (movement && outputArray) {
2843                              if (iSequence >= numberColumns_) {
2845                                   changeCost += movement * cost_[iSequence];
2846                              } else {
2848                                   changeCost += movement * cost_[iSequence];
2849                              }
2850                         }
2851                    }
2852               }
2853               dualBound_ = newBound;
2854          } else {
2855               numberInfeasibilities = -1;
2856          }
2857          return numberInfeasibilities;
2858     } else if (initialize == 1 || initialize == 3) {
2859          int iSequence;
2860          if (initialize == 3) {
2861               for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2862                    setFakeBound(iSequence, ClpSimplexDual::noFake);
2863               }
2864          }
2865          double testBound = 0.999999 * dualBound_;
2866          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2867               Status status = getStatus(iSequence);
2868               if (status == atUpperBound ||
2869                         status == atLowerBound) {
2870                    double lowerValue = lower_[iSequence];
2871                    double upperValue = upper_[iSequence];
2872                    double value = solution_[iSequence];
2873                    if (lowerValue > -largeValue_ || upperValue < largeValue_) {
2874                         if (true || lowerValue - value > -0.5 * dualBound_ ||
2875                                   upperValue - value < 0.5 * dualBound_) {
2876                              if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
2877                                   if (upperValue > lowerValue + testBound) {
2878                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2879                                             numberFake_++;
2880                                        upper_[iSequence] = lowerValue + dualBound_;
2881                                        setFakeBound(iSequence, ClpSimplexDual::upperFake);
2882                                   }
2883                              } else {
2884                                   if (lowerValue < upperValue - testBound) {
2885                                        if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2886                                             numberFake_++;
2887                                        lower_[iSequence] = upperValue - dualBound_;
2888                                        setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2889                                   }
2890                              }
2891                         } else {
2892                              if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2893                                   numberFake_++;
2894                              lower_[iSequence] = -0.5 * dualBound_;
2895                              upper_[iSequence] = 0.5 * dualBound_;
2896                              setFakeBound(iSequence, ClpSimplexDual::bothFake);
2897                              abort();
2898                         }
2899                         if (status == atUpperBound)
2900                              solution_[iSequence] = upper_[iSequence];
2901                         else
2902                              solution_[iSequence] = lower_[iSequence];
2903                    } else {
2904                         // set non basic free variables to fake bounds
2905                         // I don't think we should ever get here
2906                         CoinAssert(!("should not be here"));
2907                         lower_[iSequence] = -0.5 * dualBound_;
2908                         upper_[iSequence] = 0.5 * dualBound_;
2909                         setFakeBound(iSequence, ClpSimplexDual::bothFake);
2910                         numberFake_++;
2911                         setStatus(iSequence, atUpperBound);
2912                         solution_[iSequence] = 0.5 * dualBound_;
2913                    }
2914               } else if (status == basic) {
2915                    // make sure not at fake bound and bounds correct
2916                    setFakeBound(iSequence, ClpSimplexDual::noFake);
2917                    double gap = upper_[iSequence] - lower_[iSequence];
2918                    if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
2919                         if (iSequence < numberColumns_) {
2920                              if (columnScale_) {
2921                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2922                                   // lower
2923                                   double value = columnLower_[iSequence];
2924                                   if (value > -1.0e30) {
2925                                        value *= multiplier;
2926                                   }
2927                                   lower_[iSequence] = value;
2928                                   // upper
2929                                   value = columnUpper_[iSequence];
2930                                   if (value < 1.0e30) {
2931                                        value *= multiplier;
2932                                   }
2933                                   upper_[iSequence] = value;
2934                              } else {
2935                                   lower_[iSequence] = columnLower_[iSequence];;
2936                                   upper_[iSequence] = columnUpper_[iSequence];;
2937                              }
2938                         } else {
2939                              int iRow = iSequence - numberColumns_;
2940                              if (rowScale_) {
2941                                   // lower
2942                                   double multiplier = rhsScale_ * rowScale_[iRow];
2943                                   double value = rowLower_[iRow];
2944                                   if (value > -1.0e30) {
2945                                        value *= multiplier;
2946                                   }
2947                                   lower_[iSequence] = value;
2948                                   // upper
2949                                   value = rowUpper_[iRow];
2950                                   if (value < 1.0e30) {
2951                                        value *= multiplier;
2952                                   }
2953                                   upper_[iSequence] = value;
2954                              } else {
2955                                   lower_[iSequence] = rowLower_[iRow];;
2956                                   upper_[iSequence] = rowUpper_[iRow];;
2957                              }
2958                         }
2959                    }
2960               }
2961          }
2962
2963          return 1;
2964     } else {
2965          // just reset changed ones
2966          if (columnScale_) {
2967               int iSequence;
2968               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2969                    FakeBound fakeStatus = getFakeBound(iSequence);
2970                    if (fakeStatus != noFake) {
2971                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
2972                              // lower
2973                              double value = columnLower_[iSequence];
2974                              if (value > -1.0e30) {
2975                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2976                                   value *= multiplier;
2977                              }
2978                              columnLowerWork_[iSequence] = value;
2979                         }
2980                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
2981                              // upper
2982                              double value = columnUpper_[iSequence];
2983                              if (value < 1.0e30) {
2984                                   double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2985                                   value *= multiplier;
2986                              }
2987                              columnUpperWork_[iSequence] = value;
2988                         }
2989                    }
2990               }
2991               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2992                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
2993                    if (fakeStatus != noFake) {
2994                         if ((static_cast<int> (fakeStatus) & 1) != 0) {
2995                              // lower
2996                              double value = rowLower_[iSequence];
2997                              if (value > -1.0e30) {
2998                                   double multiplier = rhsScale_ * rowScale_[iSequence];
2999                                   value *= multiplier;
3000                              }
3001                              rowLowerWork_[iSequence] = value;
3002                         }
3003                         if ((static_cast<int> (fakeStatus) & 2) != 0) {
3004                              // upper
3005                              double value = rowUpper_[iSequence];
3006                              if (value < 1.0e30) {
3007                                   double multiplier = rhsScale_ * rowScale_[iSequence];
3008                                   value *= multiplier;
3009                              }
3010                              rowUpperWork_[iSequence] = value;
3011                         }
3012                    }
3013               }
3014          } else {
3015               int iSequence;
3016               for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3017                    FakeBound fakeStatus = getFakeBound(iSequence);
3018                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3019                         // lower
3020                         columnLowerWork_[iSequence] = columnLower_[iSequence];
3021                    }
3022                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3023                         // upper
3024                         columnUpperWork_[iSequence] = columnUpper_[iSequence];
3025                    }
3026               }
3027               for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3028                    FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3029                    if ((static_cast<int> (fakeStatus) & 1) != 0) {
3030                         // lower
3031                         rowLowerWork_[iSequence] = rowLower_[iSequence];
3032                    }
3033                    if ((static_cast<int> (fakeStatus) & 2) != 0) {
3034                         // upper
3035                         rowUpperWork_[iSequence] = rowUpper_[iSequence];
3036                    }
3037               }
3038          }
3039          return 0;
3040     }
3041}
3042int
3043ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3044                            const CoinIndexedVector * columnArray,
3045                            CoinIndexedVector * spareArray,
3046                            double acceptablePivot,
3047                            double & upperReturn, double &bestReturn, double & badFree)
3048{
3049     // do first pass to get possibles
3050     double * spare = spareArray->denseVector();
3051     int * index = spareArray->getIndices();
3052     const double * work;
3053     int number;
3054     const int * which;
3055     const double * reducedCost;
3056     // We can also see if infeasible or pivoting on free
3057     double tentativeTheta = 1.0e15;
3058     double upperTheta = 1.0e31;
3059     double freePivot = acceptablePivot;
3060     double bestPossible = 0.0;
3061     int numberRemaining = 0;
3062     int i;
3064     if ((moreSpecialOptions_ & 8) != 0) {
3065          // No free or super basic
3066          double multiplier[] = { -1.0, 1.0};
3067          double dualT = - dualTolerance_;
3068          for (int iSection = 0; iSection < 2; iSection++) {
3069
3071               unsigned char * statusArray;
3072               if (!iSection) {
3073                    work = rowArray->denseVector();
3074                    number = rowArray->getNumElements();
3075                    which = rowArray->getIndices();
3076                    reducedCost = rowReducedCost_;
3078                    statusArray = status_ + numberColumns_;
3079               } else {
3080                    work = columnArray->denseVector();
3081                    number = columnArray->getNumElements();
3082                    which = columnArray->getIndices();
3083                    reducedCost = reducedCostWork_;
3085                    statusArray = status_;
3086               }
3087
3088               for (i = 0; i < number; i++) {
3089                    int iSequence = which[i];
3090                    double alpha;
3091                    double oldValue;
3092                    double value;
3093
3094                    assert (getStatus(iSequence + addSequence) != isFree
3095                            && getStatus(iSequence + addSequence) != superBasic);
3096                    int iStatus = (statusArray[iSequence] & 3) - 1;
3097                    if (iStatus) {
3098                         double mult = multiplier[iStatus-1];
3099                         alpha = work[i] * mult;
3100                         if (alpha > 0.0) {
3101                              oldValue = reducedCost[iSequence] * mult;
3102                              value = oldValue - tentativeTheta * alpha;
3103                              if (value < dualT) {
3104                                   bestPossible = CoinMax(bestPossible, alpha);
3105                                   value = oldValue - upperTheta * alpha;
3106                                   if (value < dualT && alpha >= acceptablePivot) {
3107                                        upperTheta = (oldValue - dualT) / alpha;
3108                                        //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3109                                   }
3111                                   spare[numberRemaining] = alpha * mult;
3112                                   index[numberRemaining++] = iSequence + addSequence;
3113                              }
3114                         }
3115                    }
3116               }
3117          }
3118     } else {
3119          // some free or super basic
3120          for (int iSection = 0; iSection < 2; iSection++) {
3121
3123
3124               if (!iSection) {
3125                    work = rowArray->denseVector();
3126                    number = rowArray->getNumElements();
3127                    which = rowArray->getIndices();
3128                    reducedCost = rowReducedCost_;
3130               } else {
3131                    work = columnArray->denseVector();
3132                    number = columnArray->getNumElements();
3133                    which = columnArray->getIndices();
3134                    reducedCost = reducedCostWork_;
3136               }
3137
3138               for (i = 0; i < number; i++) {
3139                    int iSequence = which[i];
3140                    double alpha;
3141                    double oldValue;
3142                    double value;
3143                    bool keep;
3144
3146
3147                    case basic:
3148                    case ClpSimplex::isFixed:
3149                         break;
3150                    case isFree:
3151                    case superBasic:
3152                         alpha = work[i];
3153                         bestPossible = CoinMax(bestPossible, fabs(alpha));
3154                         oldValue = reducedCost[iSequence];
3155                         // If free has to be very large - should come in via dualRow
3157                         //break;
3158                         if (oldValue > dualTolerance_) {
3159                              keep = true;
3160                         } else if (oldValue < -dualTolerance_) {
3161                              keep = true;
3162                         } else {
3163                              if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3164                                   keep = true;
3165                              } else {
3166                                   keep = false;
3168                              }
3169                         }
3170                         if (keep) {
3171                              // free - choose largest
3172                              if (fabs(alpha) > freePivot) {
3173                                   freePivot = fabs(alpha);
3174                                   sequenceIn_ = iSequence + addSequence;
3175                                   theta_ = oldValue / alpha;
3176                                   alpha_ = alpha;
3177                              }
3178                         }
3179                         break;
3180                    case atUpperBound:
3181                         alpha = work[i];
3182                         oldValue = reducedCost[iSequence];
3183                         value = oldValue - tentativeTheta * alpha;
3184                         //assert (oldValue<=dualTolerance_*1.0001);
3185                         if (value > dualTolerance_) {
3186                              bestPossible = CoinMax(bestPossible, -alpha);
3187                              value = oldValue - upperTheta * alpha;
3188                              if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3189                                   upperTheta = (oldValue - dualTolerance_) / alpha;
3190                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3191                              }
3193                              spare[numberRemaining] = alpha;
3194                              index[numberRemaining++] = iSequence + addSequence;
3195                         }
3196                         break;
3197                    case atLowerBound:
3198                         alpha = work[i];
3199                         oldValue = reducedCost[iSequence];
3200                         value = oldValue - tentativeTheta * alpha;
3201                         //assert (oldValue>=-dualTolerance_*1.0001);
3202                         if (value < -dualTolerance_) {
3203                              bestPossible = CoinMax(bestPossible, alpha);
3204                              value = oldValue - upperTheta * alpha;
3205                              if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3206                                   upperTheta = (oldValue + dualTolerance_) / alpha;
3207                                   //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3208                              }
3210                              spare[numberRemaining] = alpha;
3211                              index[numberRemaining++] = iSequence + addSequence;
3212                         }
3213                         break;
3214                    }
3215               }
3216          }
3217     }
3218     upperReturn = upperTheta;
3219     bestReturn = bestPossible;
3220     return numberRemaining;
3221}
3222/*
3223   Row array has row part of pivot row (as duals so sign may be switched)
3224   Column array has column part.
3225   This chooses pivot column.
3226   Spare array will be needed when we start getting clever.
3227   We will check for basic so spare array will never overflow.
3228   If necessary will modify costs
3229*/
3230double
3231ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3232                           CoinIndexedVector * columnArray,
3233                           CoinIndexedVector * spareArray,
3234                           CoinIndexedVector * spareArray2,
3235                           double acceptablePivot,
3236                           CoinBigIndex * /*dubiousWeights*/)
3237{
3238     int numberPossiblySwapped = 0;
3239     int numberRemaining = 0;
3240
3241     double totalThru = 0.0; // for when variables flip
3242     //double saveAcceptable=acceptablePivot;
3243     //acceptablePivot=1.0e-9;
3244
3245     double bestEverPivot = acceptablePivot;
3246     int lastSequence = -1;
3247     double lastPivot = 0.0;
3248     double upperTheta;
3249     double newTolerance = dualTolerance_;
3250     //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3251     // will we need to increase tolerance
3252     bool thisIncrease = false;
3253     // If we think we need to modify costs (not if something from broad sweep)
3254     bool modifyCosts = false;
3255     // Increase in objective due to swapping bounds (may be negative)
3256     double increaseInObjective = 0.0;
3257
3258     // use spareArrays to put ones looked at in
3259     // we are going to flip flop between
3260     int iFlip = 0;
3261     // Possible list of pivots
3262     int interesting[2];
3263     // where possible swapped ones are
3264     int swapped[2];
3265     // for zeroing out arrays after
3266     int marker[2][2];
3267     // pivot elements
3268     double * array[2], * spare, * spare2;
3269     // indices
3270     int * indices[2], * index, * index2;
3271     spareArray2->clear();
3272     array[0] = spareArray->denseVector();
3273     indices[0] = spareArray->getIndices();
3274     spare = array[0];
3275     index = indices[0];
3276     array[1] = spareArray2->denseVector();
3277     indices[1] = spareArray2->getIndices();
3278     int i;
3279
3280     // initialize lists
3281     for (i = 0; i < 2; i++) {
3282          interesting[i] = 0;
3283          swapped[i] = numberColumns_;
3284          marker[i][0] = 0;
3285          marker[i][1] = numberColumns_;
3286     }
3287     /*
3288       First we get a list of possible pivots.  We can also see if the
3289       problem looks infeasible or whether we want to pivot in free variable.
3290       This may make objective go backwards but can only happen a finite
3291       number of times and I do want free variables basic.
3292
3293       Then we flip back and forth.  At the start of each iteration
3294       interesting[iFlip] should have possible candidates and swapped[iFlip]
3295       will have pivots if we decide to take a previous pivot.
3296       At end of each iteration interesting[1-iFlip] should have
3297       candidates if we go through this theta and swapped[1-iFlip]
3298       pivots if we don't go through.
3299
3300       At first we increase theta and see what happens.  We start
3301       theta at a reasonable guess.  If in right area then we do bit by bit.
3302
3303      */
3304
3305     // do first pass to get possibles
3306     upperTheta = 1.0e31;
3307     double bestPossible = 0.0;
3309     alpha_ = 0.0;
3310     if (spareIntArray_[0] >= 0) {
3311          numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3313     } else {
3315          numberRemaining = spareArray->getNumElements();
3316          spareArray->setNumElements(0);
3317          upperTheta = spareDoubleArray_[0];
3318          bestPossible = spareDoubleArray_[1];
3319          if (spareIntArray_[0] == -1) {
3320               theta_ = spareDoubleArray_[2];
3321               alpha_ = spareDoubleArray_[3];
3322               sequenceIn_ = spareIntArray_[1];
3323          } else {
3324#if 0
3325               int n = numberRemaining;
3326               double u = upperTheta;
3327               double b = bestPossible;
3328               upperTheta = 1.0e31;
3329               bestPossible = 0.0;
3330               numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3332               assert (n == numberRemaining);
3333               assert (fabs(b - bestPossible) < 1.0e-7);
3334               assert (fabs(u - upperTheta) < 1.0e-7);
3335#endif
3336          }
3337     }
3338     // switch off
3339     spareIntArray_[0] = 0;
3340     // We can also see if infeasible or pivoting on free
3341     double tentativeTheta = 1.0e25;
3342     interesting[0] = numberRemaining;
3343     marker[0][0] = numberRemaining;
3344
3345     if (!numberRemaining && sequenceIn_ < 0)
3346          return 0.0; // Looks infeasible
3347
3348     // If sum of bad small pivots too much
3349#define MORE_CAREFUL
3350#ifdef MORE_CAREFUL
3352#endif
3353     if (sequenceIn_ >= 0) {
3354          // free variable - always choose
3355     } else {
3356
3357          theta_ = 1.0e50;
3358          // now flip flop between spare arrays until reasonable theta
3359          tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3360
3361          // loops increasing tentative theta until can't go through
3362
3363          while (tentativeTheta < 1.0e22) {
3364               double thruThis = 0.0;
3365
3366               double bestPivot = acceptablePivot;
3367               int bestSequence = -1;
3368
3369               numberPossiblySwapped = numberColumns_;
3370               numberRemaining = 0;
3371
3372               upperTheta = 1.0e50;
3373
3374               spare = array[iFlip];
3375               index = indices[iFlip];
3376               spare2 = array[1-iFlip];
3377               index2 = indices[1-iFlip];
3378
3379               // try 3 different ways
3380               // 1 bias increase by ones with slightly wrong djs
3381               // 2 bias by all
3382               // 3 bias by all - tolerance
3383#define TRYBIAS 3
3384
3385
3386               double increaseInThis = 0.0; //objective increase in this loop
3387
3388               for (i = 0; i < interesting[iFlip]; i++) {
3389                    int iSequence = index[i];
3390                    double alpha = spare[i];
3391                    double oldValue = dj_[iSequence];
3392                    double value = oldValue - tentativeTheta * alpha;
3393
3394                    if (alpha < 0.0) {
3395                         //at upper bound
3396                         if (value > newTolerance) {
3397                              double range = upper_[iSequence] - lower_[iSequence];
3398                              thruThis -= range * alpha;
3399#if TRYBIAS==1
3400                              if (oldValue > 0.0)
3401                                   increaseInThis -= oldValue * range;
3402#elif TRYBIAS==2
3403                              increaseInThis -= oldValue * range;
3404#else
3405                              increaseInThis -= (oldValue + dualTolerance_) * range;
3406#endif
3407                              // goes on swapped list (also means candidates if too many)
3408                              spare2[--numberPossiblySwapped] = alpha;
3409                              index2[numberPossiblySwapped] = iSequence;
3410                              if (fabs(alpha) > bestPivot) {
3411                                   bestPivot = fabs(alpha);
3412                                   bestSequence = numberPossiblySwapped;
3413                              }
3414                         } else {
3415                              value = oldValue - upperTheta * alpha;
3416                              if (value > newTolerance && -alpha >= acceptablePivot)
3417                                   upperTheta = (oldValue - newTolerance) / alpha;
3418                              spare2[numberRemaining] = alpha;
3419                              index2[numberRemaining++] = iSequence;
3420                         }
3421                    } else {
3422                         // at lower bound
3423                         if (value < -newTolerance) {
3424                              double range = upper_[iSequence] - lower_[iSequence];
3425                              thruThis += range * alpha;
3426                              //?? is this correct - and should we look at good ones
3427#if TRYBIAS==1
3428                              if (oldValue < 0.0)
3429                                   increaseInThis += oldValue * range;
3430#elif TRYBIAS==2
3431                              increaseInThis += oldValue * range;
3432#else
3433                              increaseInThis += (oldValue - dualTolerance_) * range;
3434#endif
3435                              // goes on swapped list (also means candidates if too many)
3436                              spare2[--numberPossiblySwapped] = alpha;
3437                              index2[numberPossiblySwapped] = iSequence;
3438                              if (fabs(alpha) > bestPivot) {
3439                                   bestPivot = fabs(alpha);
3440                                   bestSequence = numberPossiblySwapped;
3441                              }
3442                         } else {
3443                              value = oldValue - upperTheta * alpha;
3444                              if (value < -newTolerance && alpha >= acceptablePivot)
3445                                   upperTheta = (oldValue + newTolerance) / alpha;
3446                              spare2[numberRemaining] = alpha;
3447                              index2[numberRemaining++] = iSequence;
3448                         }
3449                    }
3450               }
3451               swapped[1-iFlip] = numberPossiblySwapped;
3452               interesting[1-iFlip] = numberRemaining;
3453               marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3454               marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3455
3456               if (totalThru + thruThis >= fabs(dualOut_) ||
3457                         increaseInObjective + increaseInThis < 0.0) {
3458                    // We should be pivoting in this batch
3459                    // so compress down to this lot
3460                    numberRemaining = 0;
3461                    for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3462                         spare[numberRemaining] = spare2[i];
3463                         index[numberRemaining++] = index2[i];
3464                    }
3465                    interesting[iFlip] = numberRemaining;
3466                    int iTry;
3467#define MAXTRY 100
3468                    // first get ratio with tolerance
3469                    for (iTry = 0; iTry < MAXTRY; iTry++) {
3470
3471                         upperTheta = 1.0e50;
3472                         numberPossiblySwapped = numberColumns_;
3473                         numberRemaining = 0;
3474
3475                         increaseInThis = 0.0; //objective increase in this loop
3476
3477                         thruThis = 0.0;
3478
3479                         spare = array[iFlip];
3480                         index = indices[iFlip];
3481                         spare2 = array[1-iFlip];
3482                         index2 = indices[1-iFlip];
3483                         for (i = 0; i < interesting[iFlip]; i++) {
3484                              int iSequence = index[i];
3485                              double alpha = spare[i];
3486                              double oldValue = dj_[iSequence];
3487                              double value = oldValue - upperTheta * alpha;
3488
3489                              if (alpha < 0.0) {
3490                                   //at upper bound
3491                                   if (value > newTolerance) {
3492                                        if (-alpha >= acceptablePivot) {
3493                                             upperTheta = (oldValue - newTolerance) / alpha;
3494                                        }
3495                                   }
3496                              } else {
3497                                   // at lower bound
3498                                   if (value < -newTolerance) {
3499                                        if (alpha >= acceptablePivot) {
3500                                             upperTheta = (oldValue + newTolerance) / alpha;
3501                                        }
3502                                   }
3503                              }
3504                         }
3505                         bestPivot = acceptablePivot;
3506                         sequenceIn_ = -1;
3507#ifdef DUBIOUS_WEIGHTS
3508                         double bestWeight = COIN_DBL_MAX;
3509#endif
3510                         double largestPivot = acceptablePivot;
3511                         // now choose largest and sum all ones which will go through
3512                         //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3513                         // Sum of bad small pivots
3514#ifdef MORE_CAREFUL
3517#endif
3518                         // Make sure upperTheta will work (-O2 and above gives problems)
3519                         upperTheta *= 1.0000000001;
3520                         for (i = 0; i < interesting[iFlip]; i++) {
3521                              int iSequence = index[i];
3522                              double alpha = spare[i];
3523                              double value = dj_[iSequence] - upperTheta * alpha;
3525
3527
3528                              if (alpha < 0.0) {
3529                                   //at upper bound
3530                                   if (value >= 0.0) {
3532#if TRYBIAS==1
3534#elif TRYBIAS==2
3536#else
3537                                        badDj = -dj_[iSequence] - dualTolerance_;
3538#endif
3539                                   }
3540                              } else {
3541                                   // at lower bound
3542                                   if (value <= 0.0) {
3544#if TRYBIAS==1
3546#elif TRYBIAS==2
3548#else
3549                                        badDj = dj_[iSequence] - dualTolerance_;
3550#endif
3551                                   }
3552                              }
3554                                   // add to list of remaining
3555                                   spare2[numberRemaining] = alpha;
3556                                   index2[numberRemaining++] = iSequence;
3557                              } else {
3558                                   // add to list of swapped
3559                                   spare2[--numberPossiblySwapped] = alpha;
3560                                   index2[numberPossiblySwapped] = iSequence;
3561                                   // select if largest pivot
3562                                   bool take = false;
3563                                   double absAlpha = fabs(alpha);
3564#ifdef DUBIOUS_WEIGHTS
3565                                   // User could do anything to break ties here
3566                                   double weight;
3567                                   if (dubiousWeights)
3568                                        weight = dubiousWeights[iSequence];
3569                                   else
3570                                        weight = 1.0;
3571                                   weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3572                                   if (absAlpha > 2.0 * bestPivot) {
3573                                        take = true;
3574                                   } else if (absAlpha > largestPivot) {
3575                                        // could multiply absAlpha and weight
3576                                        if (weight * bestPivot < bestWeight * absAlpha)
3577                                             take = true;
3578                                   }
3579#else
3580                                   if (absAlpha > bestPivot)
3581                                        take = true;
3582#endif
3583#ifdef MORE_CAREFUL
3584                                   if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3585                                        if (alpha < 0.0) {
3586                                             //at upper bound
3587                                             if (value > dualTolerance_) {
3588                                                  double gap = upper_[iSequence] - lower_[iSequence];
3589                                                  if (gap < 1.0e20)
3590                                                       sumBadPivots += value * gap;
3591                                                  else
3593                                                  //printf("bad %d alpha %g dj at upper %g\n",
3594                                                  //     iSequence,alpha,value);
3595                                             }
3596                                        } else {
3597                                             //at lower bound
3598                                             if (value < -dualTolerance_) {
3599                                                  double gap = upper_[iSequence] - lower_[iSequence];
3600                                                  if (gap < 1.0e20)
3601                                                       sumBadPivots -= value * gap;
3602                                                  else
3604                                                  //printf("bad %d alpha %g dj at lower %g\n",
3605                                                  //     iSequence,alpha,value);
3606                                             }
3607                                        }
3608                                   }
3609#endif
3610#ifdef FORCE_FOLLOW
3611                                   if (iSequence == force_in) {
3612                                        printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3613                                        take = true;
3614                                   }
3615#endif
3616                                   if (take) {
3617                                        sequenceIn_ = numberPossiblySwapped;
3618                                        bestPivot =  absAlpha;
3619                                        theta_ = dj_[iSequence] / alpha;
3620                                        largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3621#ifdef DUBIOUS_WEIGHTS
3622                                        bestWeight = weight;
3623#endif
3624                                        //printf(" taken seq %d alpha %g weight %d\n",
3625                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3626                                   } else {
3627                                        //printf(" not taken seq %d alpha %g weight %d\n",
3628                                        //   iSequence,absAlpha,dubiousWeights[iSequence]);
3629                                   }
3630                                   double range = upper_[iSequence] - lower_[iSequence];
3631                                   thruThis += range * fabs(alpha);
3632                                   increaseInThis += badDj * range;
3633                              }
3634                         }
3635                         marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3636                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3637#ifdef MORE_CAREFUL
3638                         // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3639                         if (sumBadPivots > 1.0e4) {
3640                              if (handler_->logLevel() > 1)
3641                                   printf("maybe forcing re-factorization - sum %g  %d pivots\n", sumBadPivots,
3642                                          factorization_->pivots());
3643                              if(factorization_->pivots() > 3) {
3645                                   break;
3646                              }
3647                         }
3648#endif
3649                         swapped[1-iFlip] = numberPossiblySwapped;
3650                         interesting[1-iFlip] = numberRemaining;
3651                         // If we stop now this will be increase in objective (I think)
3652                         double increase = (fabs(dualOut_) - totalThru) * theta_;
3653                         increase += increaseInObjective;
3654                         if (theta_ < 0.0)
3655                              thruThis += fabs(dualOut_); // force using this one
3656                         if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3657                              // back
3658                              // We may need to be more careful - we could do by
3659                              // switch so we always do fine grained?
3660                              bestPivot = 0.0;
3661                         } else {
3663                              totalThru += thruThis;
3664                              increaseInObjective += increaseInThis;
3665                         }
3666                         if (bestPivot < 0.1 * bestEverPivot &&
3667                                   bestEverPivot > 1.0e-6 &&
3668                                   (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
3669                              // back to previous one
3670                              sequenceIn_ = lastSequence;
3671                              // swap regions
3672                              iFlip = 1 - iFlip;
3673                              break;
3674                         } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3675                              if (lastPivot > acceptablePivot) {
3676                                   // back to previous one
3677                                   sequenceIn_ = lastSequence;
3678                                   // swap regions
3679                                   iFlip = 1 - iFlip;
3680                              } else {
3681                                   // can only get here if all pivots too small
3682                              }
3683                              break;
3684                         } else if (totalThru >= fabs(dualOut_)) {
3685                              modifyCosts = true; // fine grain - we can modify costs
3686                              break; // no point trying another loop
3687                         } else {
3688                              lastSequence = sequenceIn_;
3689                              if (bestPivot > bestEverPivot)
3690                                   bestEverPivot = bestPivot;
3691                              iFlip = 1 - iFlip;
3692                              modifyCosts = true; // fine grain - we can modify costs
3693                         }
3694                    }
3695                    if (iTry == MAXTRY)
3696                         iFlip = 1 - iFlip; // flip back
3697                    break;
3698               } else {
3699                    // skip this lot
3700                    if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
3701                         bestEverPivot = bestPivot;
3702                         lastSequence = bestSequence;
3703                    } else {
3704                         // keep old swapped
3705                         CoinMemcpyN(array[iFlip] + swapped[iFlip],
3706                                     numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3707                         CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3708                                     numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3709                         marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3710                         swapped[1-iFlip] = swapped[iFlip];
3711                    }
3712                    increaseInObjective += increaseInThis;
3713                    iFlip = 1 - iFlip; // swap regions
3714                    tentativeTheta = 2.0 * upperTheta;
3715                    totalThru += thruThis;
3716               }
3717          }
3718
3719          // can get here without sequenceIn_ set but with lastSequence
3720          if (sequenceIn_ < 0 && lastSequence >= 0) {
3721               // back to previous one
3722               sequenceIn_ = lastSequence;
3723               // swap regions
3724               iFlip = 1 - iFlip;
3725          }
3726
3727#define MINIMUMTHETA 1.0e-18
3728          // Movement should be minimum for anti-degeneracy - unless
3729          // fixed variable out
3730          double minimumTheta;
3731          if (upperOut_ > lowerOut_)
3732               minimumTheta = MINIMUMTHETA;
3733          else
3734               minimumTheta = 0.0;
3735          if (sequenceIn_ >= 0) {
3736               // at this stage sequenceIn_ is just pointer into index array
3737               // flip just so we can use iFlip
3738               iFlip = 1 - iFlip;
3739               spare = array[iFlip];
3740               index = indices[iFlip];
3741               double oldValue;
3742               alpha_ = spare[sequenceIn_];
3743               sequenceIn_ = indices[iFlip][sequenceIn_];
3744               oldValue = dj_[sequenceIn_];
3745               theta_ = CoinMax(oldValue / alpha_, 0.0);
3746               if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3747                    // can't pivot to zero
3748#if 0
3749                    if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3750                         theta_ = minimumTheta;
3751                    } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3752                         theta_ = minimumTheta;
3753                         thisIncrease = true;
3754                    } else {
3755                         theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3756                         thisIncrease = true;
3757                    }
3758#else
3759                    theta_ = minimumTheta;
3760#endif
3761               }
3762               // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3763               //int costOffset = numberRows_+numberColumns_;
3764               if (modifyCosts && !badSumPivots) {
3765                    int i;
3766                    for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3767                         int iSequence = index[i];
3768                         double alpha = spare[i];
3769                         double value = dj_[iSequence] - theta_ * alpha;
3770
3771                         // can't be free here
3772
3773                         if (alpha < 0.0) {
3774                              //at upper bound
3775                              if (value > dualTolerance_) {
3776                                   thisIncrease = true;
3777#define MODIFYCOST 2
3778#if MODIFYCOST
3779                                   // modify cost to hit new tolerance
3780                                   double modification = alpha * theta_ - dj_[iSequence]
3781                                                         + newTolerance;
3782                                   if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3783                                        if ((specialOptions_ & 16384) != 0) {
3784                                             if (fabs(modification) < 1.0e-8)
3785                                                  modification = 0.0;
3786                                        } else if ((specialOptions_ & 2048) != 0) {
3787                                             if (fabs(modification) < 1.0e-10)
3788                                                  modification = 0.0;
3789                                        } else {
3790                                             if (fabs(modification) < 1.0e-12)
3791                                                  modification = 0.0;
3792                                        }
3793                                   }
3794                                   dj_[iSequence] += modification;
3795                                   cost_[iSequence] +=  modification;
3796                                   if (modification)
3797                                        numberChanged_ ++; // Say changed costs
3798                                   //cost_[iSequence+costOffset] += modification; // save change
3799#endif
3800                              }
3801                         } else {
3802                              // at lower bound
3803                              if (-value > dualTolerance_) {
3804                                   thisIncrease = true;
3805#if MODIFYCOST
3806                                   // modify cost to hit new tolerance
3807                                   double modification = alpha * theta_ - dj_[iSequence]
3808                                                         - newTolerance;
3809                                   //modification = CoinMax(modification,-dualTolerance_);
3810                                   //assert (fabs(modification)<1.0e-7);
3811                                   if ((specialOptions_&(2048 + 4096)) != 0) {
3812                                        if ((specialOptions_ & 2048) != 0) {
3813                                             if (fabs(modification) < 1.0e-10)
3814                                                  modification = 0.0;
3815                                        } else {
3816                                             if (fabs(modification) < 1.0e-12)
3817                                                  modification = 0.0;
3818                                        }
3819                                   }
3820                                   dj_[iSequence] += modification;
3821                                   cost_[iSequence] +=  modification;
3822                                   if (modification)
3823                                        numberChanged_ ++; // Say changed costs
3824                                   //cost_[iSequence+costOffset] += modification; // save change
3825#endif
3826                              }
3827                         }
3828                    }
3829               }
3830          }
3831     }
3832
3833#ifdef MORE_CAREFUL
3834     // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3836               fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3837          if (handler_->logLevel() > 1)
3838               printf("forcing re-factorization\n");
3839          sequenceIn_ = -1;
3840     }
3841#endif
3842     if (sequenceIn_ >= 0) {
3843          lowerIn_ = lower_[sequenceIn_];
3844          upperIn_ = upper_[sequenceIn_];
3845          valueIn_ = solution_[sequenceIn_];
3846          dualIn_ = dj_[sequenceIn_];
3847
3848          if (numberTimesOptimal_) {
3849               // can we adjust cost back closer to original
3851          }
3852#if MODIFYCOST>1
3853          // modify cost to hit zero exactly
3854          // so (dualIn_+modification)==theta_*alpha_
3855          double modification = theta_ * alpha_ - dualIn_;
3856          // But should not move objective too much ??
3857#define DONT_MOVE_OBJECTIVE
3858#ifdef DONT_MOVE_OBJECTIVE
3859          double moveObjective = fabs(modification * solution_[sequenceIn_]);
3860          double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3861          if (moveObjective > smallMove) {
3862               if (handler_->logLevel() > 1)
3863                    printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3864                           modification, solution_[sequenceIn_]);
3865               modification *= smallMove / moveObjective;
3866          }
3867#endif
3869               modification = 0.0;
3870          if ((specialOptions_&(2048 + 4096)) != 0) {
3871               if ((specialOptions_ & 16384) != 0) {
3872                    // in fast dual
3873                    if (fabs(modification) < 1.0e-7)
3874                         modification = 0.0;
3875               } else if ((specialOptions_ & 2048) != 0) {
3876                    if (fabs(modification) < 1.0e-10)
3877                         modification = 0.0;
3878               } else {
3879                    if (fabs(modification) < 1.0e-12)
3880                         modification = 0.0;
3881               }
3882          }
3883          dualIn_ += modification;
3884          dj_[sequenceIn_] = dualIn_;
3885          cost_[sequenceIn_] += modification;
3886          if (modification)
3887               numberChanged_ ++; // Say changed costs
3888          //int costOffset = numberRows_+numberColumns_;
3889          //cost_[sequenceIn_+costOffset] += modification; // save change
3890          //assert (fabs(modification)<1.0e-6);
3891#ifdef CLP_DEBUG
3892          if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3893               printf("exact %d new cost %g, change %g\n", sequenceIn_,
3894                      cost_[sequenceIn_], modification);
3895#endif
3896#endif
3897
3898          if (alpha_ < 0.0) {
3899               // as if from upper bound
3900               directionIn_ = -1;
3901               upperIn_ = valueIn_;
3902          } else {
3903               // as if from lower bound
3904               directionIn_ = 1;
3905               lowerIn_ = valueIn_;
3906          }
3907     } else {
3908          // no pivot
3909          bestPossible = 0.0;
3910          alpha_ = 0.0;
3911     }
3912     //if (thisIncrease)
3913     //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3914
3915     // clear arrays
3916
3917     for (i = 0; i < 2; i++) {
3918          CoinZeroN(array[i], marker[i][0]);
3919          CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3920     }
3921     return bestPossible;
3922}
3923#ifdef CLP_ALL_ONE_FILE
3924#undef MAXTRY
3925#endif
3926/* Checks if tentative optimal actually means unbounded
3927   Returns -3 if not, 2 if is unbounded */
3928int
3929ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3930                               CoinIndexedVector * spare,
3931                               double changeCost)
3932{
3933     int status = 2; // say unbounded
3934     factorization_->updateColumn(spare, ray);
3935     // get reduced cost
3936     int i;
3937     int number = ray->getNumElements();
3938     int * index = ray->getIndices();
3939     double * array = ray->denseVector();
3940     for (i = 0; i < number; i++) {
3941          int iRow = index[i];
3942          int iPivot = pivotVariable_[iRow];
3943          changeCost -= cost(iPivot) * array[iRow];
3944     }
3945     double way;
3946     if (changeCost > 0.0) {
3947          //try going down
3948          way = 1.0;
3949     } else if (changeCost < 0.0) {
3950          //try going up
3951          way = -1.0;
3952     } else {
3953#ifdef CLP_DEBUG
3954          printf("can't decide on up or down\n");
3955#endif
3956          way = 0.0;
3957          status = -3;
3958     }
3959     double movement = 1.0e10 * way; // some largish number
3960     double zeroTolerance = 1.0e-14 * dualBound_;
3961     for (i = 0; i < number; i++) {
3962          int iRow = index[i];
3963          int iPivot = pivotVariable_[iRow];
3964          double arrayValue = array[iRow];
3965          if (fabs(arrayValue) < zeroTolerance)
3966               arrayValue = 0.0;
3967          double newValue = solution(iPivot) + movement * arrayValue;
3968          if (newValue > upper(iPivot) + primalTolerance_ ||
3969                    newValue < lower(iPivot) - primalTolerance_)
3970               status = -3; // not unbounded
3971     }
3972     if (status == 2) {
3973          // create ray
3974          delete [] ray_;
3975          ray_ = new double [numberColumns_];
3976          CoinZeroN(ray_, numberColumns_);
3977          for (i = 0; i < number; i++) {
3978               int iRow = index[i];
3979               int iPivot = pivotVariable_[iRow];
3980               double arrayValue = array[iRow];
3981               if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3982                    ray_[iPivot] = way * array[iRow];
3983          }
3984     }
3985     ray->clear();
3986     return status;
3987}
3988//static int count_alpha=0;
3989/* Checks if finished.  Updates status */
3990void
3991ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
3992                                      double * givenDuals, ClpDataSave & saveData,
3993                                      int ifValuesPass)
3994{
3995#ifdef CLP_INVESTIGATE_SERIAL
3996     if (z_thinks > 0 && z_thinks < 2)
3997          z_thinks += 2;
3998#endif
3999     bool arraysNotCreated = (type==0);
4000     // If lots of iterations then adjust costs if large ones
4001     if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4002          double largest = 0.0;
4003          for (int i = 0; i < numberRows_; i++) {
4004               int iColumn = pivotVariable_[i];
4005               largest = CoinMax(largest, fabs(cost_[iColumn]));
4006          }
4007          if (largest > 1.0e6) {
4008               objectiveScale_ = 1.0e6 / largest;
4009               for (int i = 0; i < numberRows_ + numberColumns_; i++)
4010                    cost_[i] *= objectiveScale_;
4011          }
4012     }
4013     int numberPivots = factorization_->pivots();
4014     double realDualInfeasibilities = 0.0;
4015     if (type == 2) {
4016          if (alphaAccuracy_ != -1.0)
4017               alphaAccuracy_ = -2.0;
4018          // trouble - restore solution
4019          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4020          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4021                      numberRows_, rowActivityWork_);
4022          CoinMemcpyN(savedSolution_ ,
4023                      numberColumns_, columnActivityWork_);
4024          // restore extra stuff
4025          int dummy;
4026          matrix_->generalExpanded(this, 6, dummy);
4027          forceFactorization_ = 1; // a bit drastic but ..
4028          changeMade_++; // say something changed
4029          // get correct bounds on all variables
4030          resetFakeBounds(0);
4031     }
4032     int tentativeStatus = problemStatus_;
4033     double changeCost;
4034     bool unflagVariables = true;
4035     bool weightsSaved = false;
4036     bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4037     int dontFactorizePivots = dontFactorizePivots_;
4038     if (type == 3) {
4039          type = 1;
4040          dontFactorizePivots = 1;
4041     }
4042     if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4043          if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4044               // factorize
4045               // later on we will need to recover from singularities
4046               // also we could skip if first time
4047               // save dual weights
4048               dualRowPivot_->saveWeights(this, 1);
4049               weightsSaved = true;
4050               if (type) {
4051                    // is factorization okay?
4052                    if (internalFactorize(1)) {
4053                         // no - restore previous basis
4054                         unflagVariables = false;
4055                         assert (type == 1);
4056                         changeMade_++; // say something changed
4057                         // Keep any flagged variables
4058                         int i;
4059                         for (i = 0; i < numberRows_ + numberColumns_; i++) {
4060                              if (flagged(i))
4061                                   saveStatus_[i] |= 64; //say flagged
4062                         }
4063                         CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4064                         CoinMemcpyN(savedSolution_ + numberColumns_ ,
4065                                     numberRows_, rowActivityWork_);
4066                         CoinMemcpyN(savedSolution_ ,
4067                                     numberColumns_, columnActivityWork_);
4068                         // restore extra stuff
4069                         int dummy;
4070                         matrix_->generalExpanded(this, 6, dummy);
4071                         // get correct bounds on all variables
4072                         resetFakeBounds(1);
4073                         // need to reject something
4074                         char x = isColumn(sequenceOut_) ? 'C' : 'R';
4075                         handler_->message(CLP_SIMPLEX_FLAG, messages_)
4076                                   << x << sequenceWithin(sequenceOut_)
4077                                   << CoinMessageEol;
4078#ifdef COIN_DEVELOP
4079                         printf("flag d\n");
4080#endif
4081                         setFlagged(sequenceOut_);
4083
4084                         // Go to safe
4085                         factorization_->pivotTolerance(0.99);
4086                         forceFactorization_ = 1; // a bit drastic but ..
4087                         type = 2;
4088                         //assert (internalFactorize(1)==0);
4089                         if (internalFactorize(1)) {
4090                              CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4091                              CoinMemcpyN(savedSolution_ + numberColumns_ ,
4092                                          numberRows_, rowActivityWork_);
4093                              CoinMemcpyN(savedSolution_ ,
4094                                          numberColumns_, columnActivityWork_);
4095                              // restore extra stuff
4096                              int dummy;
4097                              matrix_->generalExpanded(this, 6, dummy);
4098                              // debug
4099                              int returnCode = internalFactorize(1);
4100                              while (returnCode) {
4101                                   // ouch
4102                                   // switch off dense
4103                                   int saveDense = factorization_->denseThreshold();
4104                                   factorization_->setDenseThreshold(0);
4105                                   // Go to safe
4106                                   factorization_->pivotTolerance(0.99);
4107                                   // make sure will do safe factorization
4108                                   pivotVariable_[0] = -1;
4109                                   returnCode = internalFactorize(2);
4110                                   factorization_->setDenseThreshold(saveDense);
4111                              }
4112                              // get correct bounds on all variables
4113                              resetFakeBounds(1);
4114                         }
4115                    }
4116               }
4117               if (problemStatus_ != -4 || numberPivots > 10)
4118                    problemStatus_ = -3;
4119          }
4120     } else {
4121          //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4122          //count_alpha++;
4123          //if ((count_alpha%5000)==0)
4124          //printf("count alpha %d\n",count_alpha);
4125     }
4126     // at this stage status is -3 or -4 if looks infeasible
4127     // get primal and dual solutions
4128#if 0
4129     {
4130          int numberTotal = numberRows_ + numberColumns_;
4131          double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4132          double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4133          double tolerance = type ? 1.0e-4 : 1.0e-8;
4134          // always if values pass
4135          double saveObj = objectiveValue_;
4136          double sumPrimal = sumPrimalInfeasibilities_;
4137          int numberPrimal = numberPrimalInfeasibilities_;
4138          double sumDual = sumDualInfeasibilities_;
4139          int numberDual = numberDualInfeasibilities_;
4140          gutsOfSolution(givenDuals, NULL);
4141          int j;
4142          double largestPrimal = tolerance;
4143          int iPrimal = -1;
4144          for (j = 0; j < numberTotal; j++) {
4145               double difference = solution_[j] - saveSol[j];
4146               if (fabs(difference) > largestPrimal) {
4147                    iPrimal = j;
4148                    largestPrimal = fabs(difference);
4149               }
4150          }
4151          double largestDual = tolerance;
4152          int iDual = -1;
4153          for (j = 0; j < numberTotal; j++) {
4154               double difference = dj_[j] - saveDj[j];
4155               if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4156                    iDual = j;
4157                    largestDual = fabs(difference);
4158               }
4159          }
4160          if (!type) {
4161               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4162                         numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4163                         fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4164                         numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4165                         fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4166                    printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4167                           type, numberIterations_, numberPivots,
4168                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4169                           largestPrimal, iPrimal,
4170                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4171                           largestDual, iDual,
4172                           saveObj, objectiveValue_);
4173          } else {
4174               if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4175                         numberPrimalInfeasibilities_ || iPrimal >= 0 ||
4176                         numberDualInfeasibilities_ || iDual >= 0)
4177                    printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4178                           type, numberIterations_, numberPivots,
4179                           numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4180                           largestPrimal, iPrimal,
4181                           numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4182                           largestDual, iDual,
4183                           saveObj, objectiveValue_);
4184          }
4185          delete [] saveSol;
4186          delete [] saveDj;
4187     }
4188#else
4189     if (type || ifValuesPass)
4190          gutsOfSolution(givenDuals, NULL);
4191#endif
4192     // If bad accuracy treat as singular
4193     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
4194          // restore previous basis
4195          unflagVariables = false;
4196          changeMade_++; // say something changed
4197          // Keep any flagged variables
4198          int i;
4199          for (i = 0; i < numberRows_ + numberColumns_; i++) {
4200               if (flagged(i))
4201                    saveStatus_[i] |= 64; //say flagged
4202          }
4203          CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4204          CoinMemcpyN(savedSolution_ + numberColumns_ ,
4205                      numberRows_, rowActivityWork_);
4206          CoinMemcpyN(savedSolution_ ,
4207                      numberColumns_, columnActivityWork_);
4208          // restore extra stuff
4209          int dummy;
4210          matrix_->generalExpanded(this, 6, dummy);
4211          // get correct bounds on all variables
4212          resetFakeBounds(1);
4213          // need to reject something
4214          char x = isColumn(sequenceOut_) ? 'C' : 'R';
4215          handler_->message(CLP_SIMPLEX_FLAG, messages_)
4216                    << x << sequenceWithin(sequenceOut_)
4217                    << CoinMessageEol;
4218#ifdef COIN_DEVELOP
4219          printf("flag e\n");
4220#endif
4221          setFlagged(sequenceOut_);
4223
4224          // Go to safer
4225          double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4226          factorization_->pivotTolerance(newTolerance);
4227          forceFactorization_ = 1; // a bit drastic but ..
4228          if (alphaAccuracy_ != -1.0)
4229               alphaAccuracy_ = -2.0;
4230          type = 2;
4231          //assert (internalFactorize(1)==0);
4232          if (internalFactorize(1)) {
4233               CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4234               CoinMemcpyN(savedSolution_ + numberColumns_ ,
4235                           numberRows_, rowActivityWork_);
4236               CoinMemcpyN(savedSolution_ ,
4237                           numberColumns_, columnActivityWork_);
4238               // restore extra stuff
4239               int dummy;
4240               matrix_->generalExpanded(this, 6, dummy);
4241               // debug
4242               int returnCode = internalFactorize(1);
4243               while (returnCode) {
4244                    // ouch
4245                    // switch off dense
4246                    int saveDense = factorization_->denseThreshold();
4247                    factorization_->setDenseThreshold(0);
4248                    // Go to safe
4249                    factorization_->pivotTolerance(0.99);
4250                    // make sure will do safe factorization
4251                    pivotVariable_[0] = -1;
4252                    returnCode = internalFactorize(2);
4253                    factorization_->setDenseThreshold(saveDense);
4254               }
4255               // get correct bounds on all variables
4256               resetFakeBounds(1);
4257          }
4258          // get primal and dual solutions
4259          gutsOfSolution(givenDuals, NULL);
4260     } else if (goodAccuracy()) {
4261          // Can reduce tolerance
4262          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4263          factorization_->pivotTolerance(newTolerance);
4264     }
4265     bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4266                                   objectiveValue_ - bestPossibleImprovement_);
4268     // Double check infeasibility if no action
4269     if (progress_.lastIterationNumber(0) == numberIterations_) {
4270          if (dualRowPivot_->looksOptimal()) {
4271               numberPrimalInfeasibilities_ = 0;
4272               sumPrimalInfeasibilities_ = 0.0;
4273          }
4274#if 1
4275     } else {
4276          double thisObj = objectiveValue_ - bestPossibleImprovement_;
4277#ifdef CLP_INVESTIGATE
4278          assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4279          if (bestPossibleImprovement_)
4280               printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4281                      thisObj);
4282#endif
4283          double lastObj = progress_.lastObjective(0);
4284#ifndef NDEBUG
4285#ifdef COIN_DEVELOP
4286          resetFakeBounds(-1);
4287#endif
4288#endif
4289#ifdef CLP_REPORT_PROGRESS
4290          ixxxxxx++;
4291          if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4292               char temp[20];
4293               sprintf(temp, "sol%d.out", ixxxxxx);
4294               printf("sol%d.out\n", ixxxxxx);
4295               FILE * fp = fopen(temp, "w");
4296               int nTotal = numberRows_ + numberColumns_;
4297               for (int i = 0; i < nTotal; i++)
4298                    fprintf(fp, "%d %d %g %g %g %g %g\n",
4299                            i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4300               fclose(fp);
4301          }
4302#endif
4303          if(!ifValuesPass && firstFree_ < 0) {
4304               double testTol = 5.0e-3;
4305               if (progress_.timesFlagged() > 10) {
4306                    testTol *= pow(2.0, progress_.timesFlagged() - 8);
4307               } else if (progress_.timesFlagged() > 5) {
4308                    testTol *= 5.0;
4309               }
4310               if (lastObj > thisObj +
4311                         testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4312                    int maxFactor = factorization_->maximumPivots();
4313                    if ((specialOptions_ & 1048576) == 0) {
4314                         if (progress_.timesFlagged() > 10)
4316                         if (maxFactor > 10 - 9) {
4317#ifdef COIN_DEVELOP
4318                              printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4319#endif
4320                              //if (forceFactorization_<0)
4321                              //forceFactorization_= maxFactor;
4322                              //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4323                              if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4324                                        largestPrimalError_ < 1.0e2) {
4325                                   // Just save costs
4326                                   // save extra copy of cost_
4327                                   int nTotal = numberRows_ + numberColumns_;
4328                                   double * temp = new double [2*nTotal];
4329                                   memcpy(temp, cost_, nTotal * sizeof(double));
4330                                   memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4331                                   delete [] cost_;
4332                                   cost_ = temp;
4333                                   objectiveWork_ = cost_;
4334                                   rowObjectiveWork_ = cost_ + numberColumns_;
4335                                   progressFlag_ |= 4;
4336                              } else {
4337                                   forceFactorization_ = 1;
4338#ifdef COIN_DEVELOP
4339                                   printf("Reducing factorization frequency - bad backwards\n");
4340#endif
4341#if 1
4342                                   unflagVariables = false;
4343                                   changeMade_++; // say something changed
4344                                   int nTotal = numberRows_ + numberColumns_;
4345                                   CoinMemcpyN(saveStatus_, nTotal, status_);
4346                                   CoinMemcpyN(savedSolution_ + numberColumns_ ,
4347                                               numberRows_, rowActivityWork_);
4348                                   CoinMemcpyN(savedSolution_ ,
4349                                               numberColumns_, columnActivityWork_);
4350                                   if ((progressFlag_ & 4) == 0) {
4351                                        // save extra copy of cost_
4352                                        double * temp = new double [2*nTotal];
4353                                        memcpy(temp, cost_, nTotal * sizeof(double));
4354                                        memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4355                                        delete [] cost_;
4356                                        cost_ = temp;
4357                                        objectiveWork_ = cost_;
4358                                        rowObjectiveWork_ = cost_ + numberColumns_;
4359                                        progressFlag_ |= 4;
4360                                   } else {
4361                                        memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4362                                   }
4363                                   // restore extra stuff
4364                                   int dummy;
4365                                   matrix_->generalExpanded(this, 6, dummy);
4366                                   double pivotTolerance = factorization_->pivotTolerance();
4367                                   if(pivotTolerance < 0.2)
4368                                        factorization_->pivotTolerance(0.2);
4369                                   else if(progress_.timesFlagged() > 2)
4370                                        factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4371                                   if (alphaAccuracy_ != -1.0)
4372                                        alphaAccuracy_ = -2.0;
4373                                   if (internalFactorize(1)) {
4374                                        CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4375                                        CoinMemcpyN(savedSolution_ + numberColumns_ ,
4376                                                    numberRows_, rowActivityWork_);
4377                                        CoinMemcpyN(savedSolution_ ,
4378                                                    numberColumns_, columnActivityWork_);
4379                                        // restore extra stuff
4380                                        int dummy;
4381                                        matrix_->generalExpanded(this, 6, dummy);
4382                                        // debug
4383                                        int returnCode = internalFactorize(1);
4384                                        while (returnCode) {
4385                                             // ouch
4386                                             // switch off dense
4387                                             int saveDense = factorization_->denseThreshold();
4388                                             factorization_->setDenseThreshold(0);
4389                                             // Go to safe
4390                                             factorization_->pivotTolerance(0.99);
4391                                             // make sure will do safe factorization
4392                                             pivotVariable_[0] = -1;
4393                                             returnCode = internalFactorize(2);
4394                                             factorization_->setDenseThreshold(saveDense);
4395                                        }
4396                                   }
4397                                   resetFakeBounds(0);
4398                                   type = 2; // so will restore weights
4399                                   // get primal and dual solutions
4400                                   gutsOfSolution(givenDuals, NULL);
4401                                   if (numberPivots < 2) {
4402                                        // need to reject something
4403                                        char x = isColumn(sequenceOut_) ? 'C' : 'R';
4404                                        handler_->message(CLP_SIMPLEX_FLAG, messages_)
4405                                                  << x << sequenceWithin(sequenceOut_)
4406                                                  << CoinMessageEol;
4407#ifdef COIN_DEVELOP
4408                                        printf("flag d\n");
4409#endif
4410                                        setFlagged(sequenceOut_);
4412                                        progress_.incrementTimesFlagged();
4413                                   }
4414                                   if (numberPivots < 10)
4416#ifdef COIN_DEVELOP
4417                                   printf("obj now %g\n", objectiveValue_);
4418#endif
4419                                   progress_.modifyObjective(objectiveValue_
4420                                                             - bestPossibleImprovement_);
4421#endif
4422                              }
4423                         }
4424                    } else {
4425                         // in fast dual give up
4426#ifdef COIN_DEVELOP
4427                         printf("In fast dual?\n");
4428#endif
4429                         problemStatus_ = 3;
4430                    }
4431               } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4432                    numberTimesOptimal_ = 0;
4433               }
4434          }
4435#endif
4436     }
4437     // Up tolerance if looks a bit odd
4438     if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4439          if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4440               int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4441               if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4442                    if (factorization_->pivotTolerance() < 0.9) {
4443                         // up tolerance
4444                         factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4445                         //printf("tol now %g\n",factorization_->pivotTolerance());
4446                         progress_.clearIterationNumbers();
4447                    }
4448               }
4449          }
4450     }
4451     // Check if looping
4452     int loop;
4453     if (!givenDuals && type != 2)
4454          loop = progress_.looping();
4455     else
4456          loop = -1;
4457     if (progress_.reallyBadTimes() > 10) {
4458          problemStatus_ = 10; // instead - try other algorithm
4459#if COIN_DEVELOP>2
4460          printf("returning at %d\n", __LINE__);
4461#endif
4462     }
4463     int situationChanged = 0;
4464     if (loop >= 0) {
4465          problemStatus_ = loop; //exit if in loop
4466          if (!problemStatus_) {
4467               // declaring victory
4468               numberPrimalInfeasibilities_ = 0;
4469               sumPrimalInfeasibilities_ = 0.0;
4470          } else {
4471               problemStatus_ = 10; // instead - try other algorithm
4472#if COIN_DEVELOP>2
4473               printf("returning at %d\n", __LINE__);
4474#endif
4475          }
4476          return;
4477     } else if (loop < -1) {
4478          // something may have changed
4479          gutsOfSolution(NULL, NULL);
4480          situationChanged = 1;
4481     }
4483     if((progressFlag_ & 2) != 0) {
4484          situationChanged = 2;
4485     }
4486     progressFlag_ &= (~3); //reset progress flag
4487     if ((progressFlag_ & 4) != 0) {
4488          // save copy of cost_
4489          int nTotal = numberRows_ + numberColumns_;
4490          memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4491     }
4492     /*if (!numberIterations_&&sumDualInfeasibilities_)
4493       printf("OBJ %g sumPinf %g sumDinf %g\n",
4494        objectiveValue(),sumPrimalInfeasibilities_,
4495        sumDualInfeasibilities_);*/
4496     // mark as having gone optimal if looks like it
4497     if (!numberPrimalInfeasibilities_&&
4498         !numberDualInfeasibilities_)
4499       progressFlag_ |= 8;
4500     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4501          handler_->message(CLP_SIMPLEX_STATUS, messages_)
4502                    << numberIterations_ << objectiveValue();
4503          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4504                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4505          handler_->printing(sumDualInfeasibilities_ > 0.0)
4506                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
4507          handler_->printing(numberDualInfeasibilitiesWithoutFree_
4508                             < numberDualInfeasibilities_)
4509                    << numberDualInfeasibilitiesWithoutFree_;
4510          handler_->message() << CoinMessageEol;
4511     }
4512#if 0
4513     printf("IT %d %g %g(%d) %g(%d)\n",
4514            numberIterations_, objectiveValue(),
4515            sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4516            sumDualInfeasibilities_, numberDualInfeasibilities_);
4517#endif
4518     double approximateObjective = objectiveValue_;
4519#ifdef CLP_REPORT_PROGRESS
4520     if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4521          char temp[20];
4522          sprintf(temp, "x_sol%d.out", ixxxxxx);
4523          FILE * fp = fopen(temp, "w");
4524          int nTotal = numberRows_ + numberColumns_;
4525          for (int i = 0; i < nTotal; i++)
4526               fprintf(fp, "%d %d %g %g %g %g %g\n",
4527                       i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4528          fclose(fp);
4529          if (ixxxxxx == ixxyyyy)
4530               exit(6);
4531     }
4532#endif
4533     realDualInfeasibilities = sumDualInfeasibilities_;
4534     double saveTolerance = dualTolerance_;
4535     // If we need to carry on cleaning variables
4536     if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4537          for (int iRow = 0; iRow < numberRows_; iRow++) {
4538               int iPivot = pivotVariable_[iRow];
4539               if (!flagged(iPivot) && pivoted(iPivot)) {
4540                    // carry on
4541                    numberPrimalInfeasibilities_ = -1;
4542                    sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4543                    sumPrimalInfeasibilities_ = 1.0;
4544                    break;
4545               }
4546          }
4547     }
4548     /* If we are primal feasible and any dual infeasibilities are on
4549        free variables then it is better to go to primal */
4550     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4551               numberDualInfeasibilities_)
4552          problemStatus_ = 10;
4553     // dual bound coming in
4554     //double saveDualBound = dualBound_;
4555     bool needCleanFake = false;
4556     while (problemStatus_ <= -3) {
4557          int cleanDuals = 0;
4558          if (situationChanged != 0)
4559               cleanDuals = 1;
4560          int numberChangedBounds = 0;
4561          int doOriginalTolerance = 0;
4562          if ( lastCleaned == numberIterations_)
4563               doOriginalTolerance = 1;
4564          // check optimal
4565          // give code benefit of doubt
4566          if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4567                    sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4568               // say optimal (with these bounds etc)
4569               numberDualInfeasibilities_ = 0;
4570               sumDualInfeasibilities_ = 0.0;
4571               numberPrimalInfeasibilities_ = 0;
4572               sumPrimalInfeasibilities_ = 0.0;
4573          }
4574          //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4575          if (dualFeasible() || problemStatus_ == -4) {
4576               progress_.modifyObjective(objectiveValue_
4577                                         - bestPossibleImprovement_);
4578#ifdef COIN_DEVELOP
4579               if (sumDualInfeasibilities_ || bestPossibleImprovement_)
4580                    printf("improve %g dualinf %g -> %g\n",
4581                           bestPossibleImprovement_, sumDualInfeasibilities_,
4582                           sumDualInfeasibilities_ * dualBound_);
4583#endif
4584               // see if cutoff reached
4585               double limit = 0.0;
4586               getDblParam(ClpDualObjectiveLimit, limit);
4587#if 0
4588               if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4589                         limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4590                    //looks infeasible on objective
4591                    if (perturbation_ == 101) {
4592                         cleanDuals = 1;
4593                         // Save costs
4594                         int numberTotal = numberRows_ + numberColumns_;
4595                         double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4596                         // make sure fake bounds are back
4597                         changeBounds(1, NULL, changeCost);
4598                         createRim4(false);
4599                         // make sure duals are current
4600                         computeDuals(givenDuals);
4601                         checkDualSolution();
4602                         if(objectiveValue()*optimizationDirection_ >
4603                                   limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4604                              perturbation_ = 102; // stop any perturbations
4605                              printf("cutoff test succeeded\n");
4606                         } else {
4607                              printf("cutoff test failed\n");
4608                              // put back
4609                              memcpy(cost_, saveCost, numberTotal * sizeof(double));
4610                              // make sure duals are current
4611                              computeDuals(givenDuals);
4612                              checkDualSolution();
4613                              progress_.modifyObjective(-COIN_DBL_MAX);
4614                              problemStatus_ = -1;
4615                         }
4616                         delete [] saveCost;
4617                    }
4618               }
4619#endif
4620               if (primalFeasible() && !givenDuals) {
4621                    // may be optimal - or may be bounds are wrong
4622                    handler_->message(CLP_DUAL_BOUNDS, messages_)
4623                              << dualBound_
4624                              << CoinMessageEol;
4625                    // save solution in case unbounded
4626                    double * saveColumnSolution = NULL;
4627                    double * saveRowSolution = NULL;
4628                    bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
4629                    if (!inCbc) {
4630                         saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4631                         saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4632                    }
4633                    numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4634                    if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4635                         //looks optimal - do we need to reset tolerance
4636                         if (perturbation_ == 101) {
4637                              perturbation_ = 102; // stop any perturbations
4638                              cleanDuals = 1;
4639                              // make sure fake bounds are back
4640                              //computeObjectiveValue();
4641                              changeBounds(1, NULL, changeCost);
4642                              //computeObjectiveValue();
4643                              createRim4(false);
4644                              // make sure duals are current
4645                              computeDuals(givenDuals);
4646                              checkDualSolution();
4647                              progress_.modifyObjective(-COIN_DBL_MAX);
4648#define DUAL_TRY_FASTER
4649#ifdef DUAL_TRY_FASTER
4650                              if (numberDualInfeasibilities_) {
4651#endif
4652                                   numberChanged_ = 1; // force something to happen
4653                                   lastCleaned = numberIterations_ - 1;
4654#ifdef DUAL_TRY_FASTER
4655                              } else {
4656                                   //double value = objectiveValue_;
4657                                   computeObjectiveValue(true);
4658                                   //printf("old %g new %g\n",value,objectiveValue_);
4659                                   //numberChanged_=1;
4660                              }
4661#endif
4662                         }
4663                         if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4664                                   (numberChanged_ || (specialOptions_ & 4096) == 0)) {
4665                              doOriginalTolerance = 2;
4666                              numberTimesOptimal_++;
4667                              changeMade_++; // say something changed
4668                              if (numberTimesOptimal_ == 1) {
4669                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4670                              } else {
4671                                   if (numberTimesOptimal_ == 2) {
4672                                        // better to have small tolerance even if slower
4673                                        factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4674                                   }
4675                                   dualTolerance_ = dblParam_[ClpDualTolerance];
4676                                   dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4677                              }
4678                              cleanDuals = 2; // If nothing changed optimal else primal
4679                         } else {
4680                              problemStatus_ = 0; // optimal
4681                              if (lastCleaned < numberIterations_ && numberChanged_) {
4682                                   handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4683                                             << CoinMessageEol;
4684                              }
4685                         }
4686                    } else {
4687                         cleanDuals = 1;
4688                         if (doOriginalTolerance == 1) {
4689                              // check unbounded
4690                              // find a variable with bad dj
4691                              int iSequence;
4692                              int iChosen = -1;
4693                              if (!inCbc) {
4694                                   double largest = 100.0 * primalTolerance_;
4695                                   for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4696                                             iSequence++) {
4697                                        double djValue = dj_[iSequence];
4698                                        double originalLo = originalLower(iSequence);
4699                                        double originalUp = originalUpper(iSequence);
4700                                        if (fabs(djValue) > fabs(largest)) {
4701                                             if (getStatus(iSequence) != basic) {
4702                                                  if (djValue > 0 && originalLo < -1.0e20) {
4703                                                       if (djValue > fabs(largest)) {
4704                                                            largest = djValue;
4705                                                            iChosen = iSequence;
4706                                                       }
4707                                                  } else if (djValue < 0 && originalUp > 1.0e20) {
4708                                                       if (-djValue > fabs(largest)) {
4709                                                            largest = djValue;
4710                                                            iChosen = iSequence;
4711                                                       }
4712                                                  }
4713                                             }
4714                                        }
4715                                   }
4716                              }
4717                              if (iChosen >= 0) {
4718                                   int iSave = sequenceIn_;
4719                                   sequenceIn_ = iChosen;
4720                                   unpack(rowArray_[1]);
4721                                   sequenceIn_ = iSave;
4722                                   // if dual infeasibilities then must be free vector so add in dual
4723                                   if (numberDualInfeasibilities_) {
4724                                        if (fabs(changeCost) > 1.0e-5)
4725                                             printf("Odd free/unbounded combo\n");
4726                                        changeCost += cost_[iChosen];
4727                                   }
4728                                   problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4729                                                                   changeCost);
4730                                   rowArray_[1]->clear();
4731                              } else {
4732                                   problemStatus_ = -3;
4733                              }
4734                              if (problemStatus_ == 2 && perturbation_ == 101) {
4735                                   perturbation_ = 102; // stop any perturbations
4736                                   cleanDuals = 1;
4737                                   createRim4(false);
4738                                   progress_.modifyObjective(-COIN_DBL_MAX);
4739                                   problemStatus_ = -1;
4740                              }
4741                              if (problemStatus_ == 2) {
4742                                   // it is unbounded - restore solution
4743                                   // but first add in changes to non-basic
4744                                   int iColumn;
4745                                   double * original = columnArray_[0]->denseVector();
4746                                   for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4747                                        if(getColumnStatus(iColumn) != basic)
4748                                             ray_[iColumn] +=
4749                                                  saveColumnSolution[iColumn] - original[iColumn];
4750                                        columnActivityWork_[iColumn] = original[iColumn];
4751                                   }
4752                                   CoinMemcpyN(saveRowSolution, numberRows_,
4753                                               rowActivityWork_);
4754                              }
4755                         } else {
4756                              doOriginalTolerance = 2;
4757                              rowArray_[0]->clear();
4758                         }
4759                    }
4760                    delete [] saveColumnSolution;
4761                    delete [] saveRowSolution;
4762               }
4763               if (problemStatus_ == -4 || problemStatus_ == -5) {
4764                    // may be infeasible - or may be bounds are wrong
4765                    numberChangedBounds = changeBounds(0, NULL, changeCost);
4766                    needCleanFake = true;
4767                    /* Should this be here as makes no difference to being feasible.
4768                       But seems to make a difference to run times. */
4769                    if (perturbation_ == 101 && 0) {
4770                         perturbation_ = 102; // stop any perturbations
4771                         cleanDuals = 1;
4772                         numberChangedBounds = 1;
4773                         // make sure fake bounds are back
4774                         changeBounds(1, NULL, changeCost);
4775                         needCleanFake = true;
4776                         createRim4(false);
4777                         progress_.modifyObjective(-COIN_DBL_MAX);
4778                    }
4779                    if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 ||
4780                              (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4781                              (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) {
4782                         problemStatus_ = 1; // infeasible
4783                         if (perturbation_ == 101) {
4784                              perturbation_ = 102; // stop any perturbations
4785                              //cleanDuals=1;
4786                              //numberChangedBounds=1;
4787                              //createRim4(false);
4788                         }
4789                    } else {
4790                         problemStatus_ = -1; //iterate
4791                         cleanDuals = 1;
4792                         if (numberChangedBounds <= 0)
4793                              doOriginalTolerance = 2;
4794                         // and delete ray which has been created
4795                         delete [] ray_;
4796                         ray_ = NULL;
4797                    }
4798
4799               }
4800          } else {
4801               cleanDuals = 1;
4802          }
4803          if (problemStatus_ < 0) {
4804               if (doOriginalTolerance == 2) {
4805                    // put back original tolerance
4806                    lastCleaned = numberIterations_;
4807                    numberChanged_ = 0; // Number of variables with changed costs
4808                    handler_->message(CLP_DUAL_ORIGINAL, messages_)
4809                              << CoinMessageEol;
4810                    perturbation_ = 102; // stop any perturbations
4811#if 0
4812                    double * xcost = new double[numberRows_+numberColumns_];
4813                    double * xlower = new double[numberRows_+numberColumns_];
4814                    double * xupper = new double[numberRows_+numberColumns_];
4815                    double * xdj = new double[numberRows_+numberColumns_];
4816                    double * xsolution = new double[numberRows_+numberColumns_];
4817                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4818                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4819                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4820                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4821                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4822#endif
4823                    createRim4(false);
4824                    progress_.modifyObjective(-COIN_DBL_MAX);
4825                    // make sure duals are current
4826                    computeDuals(givenDuals);
4827                    checkDualSolution();
4828#if 0
4829                    int i;
4830                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4831                         if (cost_[i] != xcost[i])
4832                              printf("** %d old cost %g new %g sol %g\n",
4833                                     i, xcost[i], cost_[i], solution_[i]);
4834                         if (lower_[i] != xlower[i])
4835                              printf("** %d old lower %g new %g sol %g\n",
4836                                     i, xlower[i], lower_[i], solution_[i]);
4837                         if (upper_[i] != xupper[i])
4838                              printf("** %d old upper %g new %g sol %g\n",
4839                                     i, xupper[i], upper_[i], solution_[i]);
4840                         if (dj_[i] != xdj[i])
4841                              printf("** %d old dj %g new %g sol %g\n",
4842                                     i, xdj[i], dj_[i], solution_[i]);
4843                         if (solution_[i] != xsolution[i])
4844                              printf("** %d old solution %g new %g sol %g\n",
4845                                     i, xsolution[i], solution_[i], solution_[i]);
4846                    }
4847                    //delete [] xcost;
4848                    //delete [] xupper;
4849                    //delete [] xlower;
4850                    //delete [] xdj;
4851                    //delete [] xsolution;
4852#endif
4853                    // put back bounds as they were if was optimal
4854                    if (doOriginalTolerance == 2 && cleanDuals != 2) {
4855                         changeMade_++; // say something changed
4856                         /* We may have already changed some bounds in this function
4857                            so save numberFake_ and add in.
4858
4859                            Worst that can happen is that we waste a bit of time  - but it must be finite.
4860                         */
4861                         //int saveNumberFake = numberFake_;
4862                         //resetFakeBounds(-1);
4863                         changeBounds(3, NULL, changeCost);
4864                         needCleanFake = true;
4865                         //numberFake_ += saveNumberFake;
4866                         //resetFakeBounds(-1);
4867                         cleanDuals = 2;
4868                         //cleanDuals=1;
4869                    }
4870#if 0
4871                    //int i;
4872                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4873                         if (cost_[i] != xcost[i])
4874                              printf("** %d old cost %g new %g sol %g\n",
4875                                     i, xcost[i], cost_[i], solution_[i]);
4876                         if (lower_[i] != xlower[i])
4877                              printf("** %d old lower %g new %g sol %g\n",
4878                                     i, xlower[i], lower_[i], solution_[i]);
4879                         if (upper_[i] != xupper[i])
4880                              printf("** %d old upper %g new %g sol %g\n",
4881                                     i, xupper[i], upper_[i], solution_[i]);
4882                         if (dj_[i] != xdj[i])
4883                              printf("** %d old dj %g new %g sol %g\n",
4884                                     i, xdj[i], dj_[i], solution_[i]);
4885                         if (solution_[i] != xsolution[i])
4886                              printf("** %d old solution %g new %g sol %g\n",
4887                                     i, xsolution[i], solution_[i], solution_[i]);
4888                    }
4889                    delete [] xcost;
4890                    delete [] xupper;
4891                    delete [] xlower;
4892                    delete [] xdj;
4893                    delete [] xsolution;
4894#endif
4895               }
4896               if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4897                    // make sure dual feasible
4898                    // look at all rows and columns
4899                    rowArray_[0]->clear();
4900                    columnArray_[0]->clear();
4901                    double objectiveChange = 0.0;
4902                    double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4903                    if (!numberIterations_) {
4904                      int nTotal = numberRows_ + numberColumns_;
4905                      if (arraysNotCreated) {
4906                        // create save arrays
4907                        delete [] saveStatus_;
4908                        delete [] savedSolution_;
4909                        saveStatus_ = new unsigned char [nTotal];
4910                        savedSolution_ = new double [nTotal];
4911                        arraysNotCreated = false;
4912                      }
4913                      // save arrays
4914                      CoinMemcpyN(status_, nTotal, saveStatus_);
4915                      CoinMemcpyN(rowActivityWork_,
4916                                  numberRows_, savedSolution_ + numberColumns_);
4917                      CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4918                    }
4919#if 0
4920                    double * xcost = new double[numberRows_+numberColumns_];
4921                    double * xlower = new double[numberRows_+numberColumns_];
4922                    double * xupper = new double[numberRows_+numberColumns_];
4923                    double * xdj = new double[numberRows_+numberColumns_];
4924                    double * xsolution = new double[numberRows_+numberColumns_];
4925                    CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4926                    CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4927                    CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4928                    CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4929                    CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4930#endif
4931                    if (givenDuals)
4932                         dualTolerance_ = 1.0e50;
4933                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4934                                      0.0, objectiveChange, true);
4935                    dualTolerance_ = saveTolerance;
4936#if 0
4937                    int i;
4938                    for (i = 0; i < numberRows_ + numberColumns_; i++) {
4939                         if (cost_[i] != xcost[i])
4940                              printf("** %d old cost %g new %g sol %g\n",
4941                                     i, xcost[i], cost_[i], solution_[i]);
4942                         if (lower_[i] != xlower[i])
4943                              printf("** %d old lower %g new %g sol %g\n",
4944                                     i, xlower[i], lower_[i], solution_[i]);
4945                         if (upper_[i] != xupper[i])
4946                              printf("** %d old upper %g new %g sol %g\n",
4947                                     i, xupper[i], upper_[i], solution_[i]);
4948                         if (dj_[i] != xdj[i])
4949                              printf("** %d old dj %g new %g sol %g\n",
4950                                     i, xdj[i], dj_[i], solution_[i]);
4951                         if (solution_[i] != xsolution[i])
4952                              printf("** %d old solution %g new %g sol %g\n",
4953                                     i, xsolution[i], solution_[i], solution_[i]);
4954                    }
4955                    delete [] xcost;
4956                    delete [] xupper;
4957                    delete [] xlower;
4958                    delete [] xdj;
4959                    delete [] xsolution;
4960#endif
4961                    // for now - recompute all
4962                    gutsOfSolution(NULL, NULL);
4963                    if (givenDuals)
4964                         dualTolerance_ = 1.0e50;
4965                    updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4966                                      0.0, objectiveChange, true);
4967                    dualTolerance_ = saveTolerance;
4968                    if (!numberIterations_ && sumPrimalInfeasibilities_ >
4969                        1.0e5*(savePrimalInfeasibilities+1.0e3) &&
4970                        (moreSpecialOptions_ & 256) == 0) {
4971                      // Use primal
4972                      int nTotal = numberRows_ + numberColumns_;
4973                      CoinMemcpyN(saveStatus_, nTotal, status_);
4974                      CoinMemcpyN(savedSolution_ + numberColumns_ ,
4975                                  numberRows_, rowActivityWork_);
4976                      CoinMemcpyN(savedSolution_ ,
4977                                  numberColumns_, columnActivityWork_);
4978                      //numberDualInfeasibilities_=1;
4979                      problemStatus_ = 10;
4980                      situationChanged = 0;
4981                      //numberPrimalInfeasibilities_=1;
4982                    }
4983                    //assert(numberDualInfeasibilitiesWithoutFree_==0);
4984                    if (numberDualInfeasibilities_) {
4985                         if (numberPrimalInfeasibilities_ || numberPivots) {
4986                              problemStatus_ = -1; // carry on as normal
4987                         } else {
4988                              problemStatus_ = 10; // try primal
4989#if COIN_DEVELOP>1
4990                              printf("returning at %d\n", __LINE__);
4991#endif
4992                         }
4993                    } else if (situationChanged == 2) {
4994                         problemStatus_ = -1; // carry on as normal
4995                         // need to reset bounds
4996                         changeBounds(3, NULL, changeCost);
4997                    }
4998                    situationChanged = 0;
4999               } else {
5000                    // iterate
5001                    if (cleanDuals != 2) {
5002                         problemStatus_ = -1;
5003                    } else {
5004                         problemStatus_ = 10; // try primal
5005#if COIN_DEVELOP>2
5006                         printf("returning at %d\n", __LINE__);
5007#endif
5008                    }
5009               }
5010          }
5011     }
5012     // unflag all variables (we may want to wait a bit?)
5013     if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5014          int iRow;
5015          int numberFlagged = 0;
5016          for (iRow = 0; iRow < numberRows_; iRow++) {
5017               int iPivot = pivotVariable_[iRow];
5018               if (flagged(iPivot)) {
5019                    numberFlagged++;
5020                    clearFlagged(iPivot);
5021               }
5022          }
5023#ifdef COIN_DEVELOP
5024          if (numberFlagged) {
5025               printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5026                      problemStatus_, numberPrimalInfeasibilities_,
5027                      numberTimesOptimal_);
5028          }
5029#endif
5030          unflagVariables = numberFlagged > 0;
5031          if (numberFlagged && !numberPivots) {
5032               /* looks like trouble as we have not done any iterations.
5033               Try changing pivot tolerance then give it a few goes and give up */
5034               if (factorization_->pivotTolerance() < 0.9) {
5035                    factorization_->pivotTolerance(0.99);
5036                    problemStatus_ = -1;
5037               } else if (numberTimesOptimal_ < 3) {
5038                    numberTimesOptimal_++;
5039                    problemStatus_ = -1;
5040               } else {
5041                    unflagVariables = false;
5042                    //secondaryStatus_ = 1; // and say probably infeasible
5043                    if ((moreSpecialOptions_ & 256) == 0) {
5044                         // try primal
5045                         problemStatus_ = 10;
5046                    } else {
5047                         // almost certainly infeasible
5048                         problemStatus_ = 1;
5049                    }
5050#if COIN_DEVELOP>1
5051                    printf("returning at %d\n", __LINE__);
5052#endif
5053               }
5054          }
5055     }
5056     if (problemStatus_ < 0) {
5057          if (needCleanFake) {
5058               double dummyChangeCost = 0.0;
5059               changeBounds(3, NULL, dummyChangeCost);
5060          }
5061#if 0
5062          if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5063                    CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5064          } else {
5065               lastObjectiveValue_ = objectiveValue_;
5066          }
5067#endif
5068          if (type == 0 || type == 1) {
5069               if (!type && arraysNotCreated) {
5070                    // create save arrays
5071                    delete [] saveStatus_;
5072                    delete [] savedSolution_;
5073                    saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5074                    savedSolution_ = new double [numberRows_+numberColumns_];
5075               }
5076               // save arrays
5077               CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5078               CoinMemcpyN(rowActivityWork_,
5079                           numberRows_, savedSolution_ + numberColumns_);
5080               CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5081               // save extra stuff
5082               int dummy;
5083               matrix_->generalExpanded(this, 5, dummy);
5084          }
5085          if (weightsSaved) {
5086               // restore weights (if saved) - also recompute infeasibility list
5087               if (!reallyBadProblems && (largestPrimalError_ < 100.0 || numberPivots > 10)) {
5088                    if (tentativeStatus > -3)
5089                         dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5090                    else
5091                         dualRowPivot_->saveWeights(this, 3);
5092               } else {
5093                    // reset weights or scale back
5094                    dualRowPivot_->saveWeights(this, 6);
5095               }
5096          } else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5097               dualRowPivot_->saveWeights(this, 3);
5098          }
5099     }
5100     // see if cutoff reached
5101     double limit = 0.0;
5102     getDblParam(ClpDualObjectiveLimit, limit);
5103#if 0
5104     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5105               limit + 100.0) {
5106          printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5107                 limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5108     }
5109#endif
5110     if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5111               limit && !numberAtFakeBound()) {
5112          bool looksInfeasible = !numberDualInfeasibilities_;
5113          if (objectiveValue()*optimizationDirection_ > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5114                    sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5115               looksInfeasible = true;
5116          if (looksInfeasible) {
5117               // Even if not perturbed internal costs may have changed
5118               // be careful
5119               if (true || numberIterations_) {
5120                    if(computeInternalObjectiveValue() > limit) {
5121                         problemStatus_ = 1;
5122                         secondaryStatus_ = 1; // and say was on cutoff
5123                    }
5124               } else {
5125                    problemStatus_ = 1;
5126                    secondaryStatus_ = 1; // and say was on cutoff
5127               }
5128          }
5129     }
5130     // If we are in trouble and in branch and bound give up
5131     if ((specialOptions_ & 1024) != 0) {
5133          if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5135          } else if (largestPrimalError_ > 1.0e-2
5136                     && objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5138          }
5140               if (factorization_->pivotTolerance() < 0.9) {
5141                    // up tolerance
5142                    factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5143               } else if (numberIterations_ > 10000) {
5144                    if (handler_->logLevel() > 2)
5146                    problemStatus_ = 1;
5147                    secondaryStatus_ = 1; // and say was on cutoff
5148               } else if (largestPrimalError_ > 1.0e5) {
5149                    {
5150                         int iBigB = -1;
5151                         double bigB = 0.0;
5152                         int iBigN = -1;
5153                         double bigN = 0.0;
5154                         for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5155                              double value = fabs(solution_[i]);
5156                              if (getStatus(i) == basic) {
5157                                   if (value > bigB) {
5158                                        bigB = value;
5159                                        iBigB = i;
5160                                   }
5161                              } else {
5162                                   if (value > bigN) {
5163                                        bigN = value;
5164                                        iBigN = i;
5165                                   }
5166                              }
5167                         }
5168#ifdef CLP_INVESTIGATE
5169                         if (bigB > 1.0e8 || bigN > 1.0e8) {
5170                              if (handler_->logLevel() > 0)
5171                                   printf("it %d - basic %d %g, nonbasic %d %g\n",
5172                                          numberIterations_, iBigB, bigB, iBigN, bigN);
5173                         }
5174#endif
5175                    }
5176#if COIN_DEVELOP!=2
5177                    if (handler_->logLevel() > 2)
5178#endif
5180                    allSlackBasis(true);
5181                    problemStatus_ = 10;
5182               }
5183          }
5184     }
5185     if (problemStatus_ < 0 && !changeMade_) {
5186          problemStatus_ = 4; // unknown
5187     }
5188     lastGoodIteration_ = numberIterations_;
5189     if (numberIterations_ > lastBadIteration_ + 100)
5190          moreSpecialOptions_ &= ~16; // clear check accuracy flag
5191     if (problemStatus_ < 0) {
5192          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
5193          if (sumDualInfeasibilities_)
5194               numberDualInfeasibilities_ = 1;
5195     }
5196#ifdef CLP_REPORT_PROGRESS
5197     if (ixxxxxx > ixxyyyy - 3) {
5198          printf("objectiveValue_ %g\n", objectiveValue_);
5199          handler_->setLogLevel(63);
5200          int nTotal = numberColumns_ + numberRows_;
5201          double newObj = 0.0;
5202          for (int i = 0; i < nTotal; i++) {
5203               if (solution_[i])
5204                    newObj += solution_[i] * cost_[i];
5205          }
5206          printf("xxx obj %g\n", newObj);
5207          // for now - recompute all
5208          gutsOfSolution(NULL, NULL);
5209          newObj = 0.0;
5210          for (int i = 0; i < nTotal; i++) {
5211               if (solution_[i])
5212                    newObj += solution_[i] * cost_[i];
5213          }
5214          printf("yyy obj %g %g\n", newObj, objectiveValue_);
5215          progress_.modifyObjective(objectiveValue_
5216                                    - bestPossibleImprovement_);
5217     }
5218#endif
5219#if 1
5220     double thisObj = progress_.lastObjective(0);
5221     double lastObj = progress_.lastObjective(1);
5222     if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4
5223               && givenDuals == NULL && firstFree_ < 0) {
5224          int maxFactor = factorization_->maximumPivots();
5225          if (maxFactor > 10) {
5226               if (forceFactorization_ < 0)
5227                    forceFactorization_ = maxFactor;
5228               forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
5229               //printf("Reducing factorization frequency\n");
5230          }
5231     }
5232#endif
5233     // Allow matrices to be sorted etc
5234     int fake = -999; // signal sort
5235     matrix_->correctSequence(this, fake, fake);
5236     if (alphaAccuracy_ > 0.0)
5237          alphaAccuracy_ = 1.0;
5238     // If we are stopping - use plausible objective
5239     // Maybe only in fast dual
5240     if (problemStatus_ > 2)
5241          objectiveValue_ = approximateObjective;
5242     if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
5243         fabs(objectiveValue_) > 1.0e10 )
5244       problemStatus_ = 10; // infeasible - but has looked feasible
5245}
5246/* While updateDualsInDual sees what effect is of flip
5247   this does actual flipping.
5248   If change >0.0 then value in array >0.0 => from lower to upper
5249*/
5250void
5251ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
5252                           CoinIndexedVector * columnArray)
5253{
5254     int number;
5255     int * which;
5256
5257     int iSection;