source: stable/1.16/Clp/src/ClpSimplexDual.cpp @ 2265

Last change on this file since 2265 was 2265, checked in by forrest, 3 years ago

fix nan dual infeasible

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