source: branches/devel-1/ClpSimplexDual.cpp @ 19

Last change on this file since 19 was 19, checked in by ladanyi, 17 years ago

reordering Clp

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5/* Notes on implementation of dual simplex algorithm.
6
7   When dual feasible:
8
9   If primal feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to leave basis (normally going to nearest bound) (B).  We
11   now need to find an incoming variable which will leave problem
12   dual feasible so we get the row of the tableau corresponding to
13   the basic variable (with the correct sign depending if basic variable
14   above or below feasibility region - as that affects whether reduced
15   cost on outgoing variable has to be positive or negative).
16
17   We now perform a ratio test to determine which incoming variable will
18   preserve dual feasibility (C).  If no variable found then problem
19   is infeasible (in primal sense).  If there is a variable, we then
20   perform pivot and repeat.  Trivial?
21
22   -------------------------------------------
23
24   A) How do we get dual feasible?  If all variables have bounds then
25   it is trivial to get feasible by putting non-basic variables to
26   correct bounds.  OSL did not have a phase 1/phase 2 approach but
27   instead effectively put fake bounds on variables and this is the
28   approach here, although I had hoped to make it cleaner.
29
30   If there is a weight of X on getting dual feasible:
31     Non-basic variables with negative reduced costs are put to
32     lesser of their upper bound and their lower bound + X.
33     Similarly, mutatis mutandis, for positive reduced costs.
34
35   Free variables should normally be in basis, otherwise I have
36   coding which may be able to come out (and may not be correct).
37
38   In OSL, this weight was changed heuristically, here at present
39   it is only increased if problem looks finished.  If problem is
40   feasible I check for unboundedness.  If not unbounded we
41   could play with going into primal.  As long as weights increase
42   any algorithm would be finite.
43   
44   B) Which outgoing variable to choose is a virtual base class.
45   For difficult problems steepest edge is preferred while for
46   very easy (large) problems we will need partial scan.
47
48   C) Sounds easy, but this is hardest part of algorithm.
49      1) Instead of stopping at first choice, we may be able
50      to flip that variable to other bound and if objective
51      still improving choose again.  These mini iterations can
52      increase speed by orders of magnitude but we may need to
53      go to more of a bucket choice of variable rather than looking
54      at them one by one (for speed).
55      2) Accuracy.  Reduced costs may be of wrong sign but less than
56      tolerance.  Pivoting on these makes objective go backwards.
57      OSL modified cost so a zero move was made, Gill et al
58      (in primal analogue) modified so a strictly positive move was
59      made.  It is not quite as neat in dual but that is what we
60      try and do.  The two problems are that re-factorizations can
61      change reduced costs above and below tolerances and that when
62      finished we need to reset costs and try again.
63      3) Degeneracy.  Gill et al helps but may not be enough.  We
64      may need more.  Also it can improve speed a lot if we perturb
65      the costs significantly. 
66
67  References:
68     Forrest and Goldfarb, Steepest-edge simplex algorithms for
69       linear programming - Mathematical Programming 1992
70     Forrest and Tomlin, Implementing the simplex method for
71       the Optimization Subroutine Library - IBM Systems Journal 1992
72     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
73       Procedure for Linear and Nonlinear Programming SOL report 1988
74
75
76  TODO:
77 
78  a) Better recovery procedures.  At present I never check on forward
79     progress.  There is checkpoint/restart with reducing
80     re-factorization frequency, but this is only on singular
81     factorizations.
82  b) Fast methods for large easy problems (and also the option for
83     the code to automatically choose which method).
84  c) We need to be able to stop in various ways for OSI - this
85     is fairly easy.
86
87 */
88
89#if defined(_MSC_VER)
90// Turn off compiler warning about long names
91#  pragma warning(disable:4786)
92#endif
93
94#include <math.h>
95
96#include "CoinHelperFunctions.hpp"
97#include "ClpSimplexDual.hpp"
98#include "ClpFactorization.hpp"
99#include "CoinPackedMatrix.hpp"
100#include "CoinIndexedVector.hpp"
101#include "CoinWarmStartBasis.hpp"
102#include "ClpDualRowDantzig.hpp"
103#include "ClpMessage.hpp"
104#include <cfloat>
105#include <cassert>
106#include <string>
107#include <stdio.h>
108#include <iostream>
109// This returns a non const array filled with input from scalar
110// or actual array
111template <class T> inline T*
112copyOfArray( const T * array, const int size, T value)
113{
114  T * arrayNew = new T[size];
115  if (array) 
116    CoinDisjointCopyN(array,size,arrayNew);
117  else
118    CoinFillN ( arrayNew, size,value);
119  return arrayNew;
120}
121
122// This returns a non const array filled with actual array (or NULL)
123template <class T> inline T*
124copyOfArray( const T * array, const int size)
125{
126  if (array) {
127    T * arrayNew = new T[size];
128    CoinDisjointCopyN(array,size,arrayNew);
129    return arrayNew;
130  } else {
131    return NULL;
132  }
133}
134// dual
135int ClpSimplexDual::dual ( )
136{
137
138  /* *** Method
139     This is a vanilla version of dual simplex.
140
141     It tries to be a single phase approach with a weight of 1.0 being
142     given to getting optimal and a weight of dualBound_ being
143     given to getting dual feasible.  In this version I have used the
144     idea that this weight can be thought of as a fake bound.  If the
145     distance between the lower and upper bounds on a variable is less
146     than the feasibility weight then we are always better off flipping
147     to other bound to make dual feasible.  If the distance is greater
148     then we make up a fake bound dualBound_ away from one bound.
149     If we end up optimal or primal infeasible, we check to see if
150     bounds okay.  If so we have finished, if not we increase dualBound_
151     and continue (after checking if unbounded). I am undecided about
152     free variables - there is coding but I am not sure about it.  At
153     present I put them in basis anyway.
154
155     The code is designed to take advantage of sparsity so arrays are
156     seldom zeroed out from scratch or gone over in their entirety.
157     The only exception is a full scan to find outgoing variable.  This
158     will be changed to keep an updated list of infeasibilities (or squares
159     if steepest edge).  Also on easy problems we don't need full scan - just
160     pick first reasonable.
161
162     One problem is how to tackle degeneracy and accuracy.  At present
163     I am using the modification of costs which I put in OSL and which was
164     extended by Gill et al.  I am still not sure of the exact details.
165
166     The flow of dual is three while loops as follows:
167
168     while (not finished) {
169
170       while (not clean solution) {
171
172          Factorize and/or clean up solution by flipping variables so
173          dual feasible.  If looks finished check fake dual bounds.
174          Repeat until status is iterating (-1) or finished (0,1,2)
175
176       }
177
178       while (status==-1) {
179
180         Iterate until no pivot in or out or time to re-factorize.
181
182         Flow is:
183
184         choose pivot row (outgoing variable).  if none then
185         we are primal feasible so looks as if done but we need to
186         break and check bounds etc.
187
188         Get pivot row in tableau
189
190         Choose incoming column.  If we don't find one then we look
191         primal infeasible so break and check bounds etc.  (Also the
192         pivot tolerance is larger after any iterations so that may be
193         reason)
194
195         If we do find incoming column, we may have to adjust costs to
196         keep going forwards (anti-degeneracy).  Check pivot will be stable
197         and if unstable throw away iteration (we will need to implement
198         flagging of basic variables sometime) and break to re-factorize.
199         If minor error re-factorize after iteration.
200
201         Update everything (this may involve flipping variables to stay
202         dual feasible.
203
204       }
205
206     }
207
208     At present we never check we are going forwards.  I overdid that in
209     OSL so will try and make a last resort.
210
211     Needs partial scan pivot out option.
212     Needs dantzig, uninitialized and full steepest edge options (can still
213     use partial scan)
214
215     May need other anti-degeneracy measures, especially if we try and use
216     loose tolerances as a way to solve in fewer iterations.
217
218     I like idea of dynamic scaling.  This gives opportunity to decouple
219     different implications of scaling for accuracy, iteration count and
220     feasibility tolerance.
221
222  */
223
224
225  // sanity check
226  assert (numberRows_==matrix_->getNumRows());
227  assert (numberColumns_==matrix_->getNumCols());
228  // for moment all arrays must exist
229  assert(columnLower_);
230  assert(columnUpper_);
231  assert(rowLower_);
232  assert(rowUpper_);
233
234#ifdef CLP_DEBUG
235  int debugIteration=-1;
236#endif
237
238  algorithm_ = -1;
239  dualTolerance_=dblParam_[ClpDualTolerance];
240  primalTolerance_=dblParam_[ClpPrimalTolerance];
241
242  // put in standard form (and make row copy)
243  // create modifiable copies of model rim and do optional scaling
244  createRim(7+8+16,true);
245
246  // save dual bound
247  double saveDualBound_ = dualBound_;
248
249  int iRow,iColumn;
250  // Do initial factorization
251  // and set certain stuff
252  // We can either set increasing rows so ...IsBasic gives pivot row
253  // or we can just increment iBasic one by one
254  // for now let ...iBasic give pivot row
255  factorization_->increasingRows(2);
256  // row activities have negative sign
257  factorization_->slackValue(-1.0);
258  factorization_->zeroTolerance(1.0e-13);
259  // save if sparse factorization wanted
260  int saveSparse = factorization_->sparseThreshold();
261
262  int factorizationStatus = internalFactorize(0);
263  if (factorizationStatus<0)
264    return 1; // some error
265  else if (factorizationStatus)
266    handler_->message(CLP_SINGULARITIES,messages_)
267    <<factorizationStatus
268    <<CoinMessageEol;
269
270  // If user asked for perturbation - do it
271  int savePerturbation = perturbation_;
272
273  if (perturbation_<100) 
274    perturb();
275
276  double objectiveChange;
277  // for dual we will change bounds using dualBound_
278  // for this we need clean basis so it is after factorize
279  gutsOfSolution(rowActivityWork_,columnActivityWork_);
280
281  numberFake_ =0; // Number of variables at fake bounds
282  changeBounds(true,NULL,objectiveChange);
283
284  problemStatus_ = -1;
285  numberIterations_=0;
286
287  int lastCleaned=0; // last time objective or bounds cleaned up
288
289  // number of times we have declared optimality
290  numberTimesOptimal_=0;
291
292  // This says whether to restore things etc
293  int factorType=0;
294  /*
295    Status of problem:
296    0 - optimal
297    1 - infeasible
298    2 - unbounded
299    -1 - iterating
300    -2 - factorization wanted
301    -3 - redo checking without factorization
302    -4 - looks infeasible
303  */
304  while (problemStatus_<0) {
305    // clear
306    for (iRow=0;iRow<4;iRow++) {
307      rowArray_[iRow]->clear();
308    }   
309   
310    for (iColumn=0;iColumn<2;iColumn++) {
311      columnArray_[iColumn]->clear();
312    }   
313
314    // give matrix (and model costs and bounds a chance to be
315    // refreshed (normally null)
316    matrix_->refresh(this);
317    // If getting nowhere - why not give it a kick
318#if 0
319    // does not seem to work too well - do some more work
320    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
321      perturb();
322#endif
323    // may factorize, checks if problem finished
324    statusOfProblemInDual(lastCleaned,factorType);
325
326    // Say good factorization
327    factorType=1;
328    if (saveSparse) {
329      // use default at present
330      factorization_->sparseThreshold(0);
331      factorization_->goSparse();
332    }
333
334    // Do iterations
335    whileIterating();
336  }
337
338  assert(!numberFake_||problemStatus_); // all bounds should be okay
339  // at present we are leaving factorization around
340  // maybe we should empty it
341  deleteRim();
342  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
343    <<objectiveValue()
344    <<CoinMessageEol;
345  // Restore any saved stuff
346  perturbation_ = savePerturbation;
347  factorization_->sparseThreshold(saveSparse);
348  dualBound_ = saveDualBound_;
349  return problemStatus_;
350}
351/* Reasons to come out:
352   -1 iterations etc
353   -2 inaccuracy
354   -3 slight inaccuracy (and done iterations)
355   +0 looks optimal (might be unbounded - but we will investigate)
356   +1 looks infeasible
357   +3 max iterations
358 */
359int
360ClpSimplexDual::whileIterating()
361{
362  // status stays at -1 while iterating, >=0 finished, -2 to invert
363  // status -3 to go to top without an invert
364  int returnCode = -1;
365  while (problemStatus_==-1) {
366#ifdef CLP_DEBUG
367    {
368      int i;
369      for (i=0;i<4;i++) {
370        rowArray_[i]->checkClear();
371      }   
372      for (i=0;i<2;i++) {
373        columnArray_[i]->checkClear();
374      }   
375    }     
376#endif
377#if CLP_DEBUG>2
378    // very expensive
379    if (numberIterations_>0&&numberIterations_<-801) {
380      handler_->setLogLevel(63);
381      double saveValue = objectiveValue_;
382      double * saveRow1 = new double[numberRows_];
383      double * saveRow2 = new double[numberRows_];
384      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
385      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
386      double * saveColumn1 = new double[numberColumns_];
387      double * saveColumn2 = new double[numberColumns_];
388      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
389      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
390      gutsOfSolution(rowActivityWork_,columnActivityWork_);
391      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
392             numberIterations_,
393             saveValue,objectiveValue_,sumDualInfeasibilities_);
394      if (saveValue>objectiveValue_+1.0e-2)
395        printf("**bad**\n");
396      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
397      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
398      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
399      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
400      delete [] saveRow1;
401      delete [] saveRow2;
402      delete [] saveColumn1;
403      delete [] saveColumn2;
404      objectiveValue_=saveValue;
405    }
406#endif
407#ifdef CLP_DEBUG
408    {
409      int iSequence, number=numberRows_+numberColumns_;
410      for (iSequence=0;iSequence<number;iSequence++) {
411        double lowerValue=lower_[iSequence];
412        double upperValue=upper_[iSequence];
413        double value=solution_[iSequence];
414        if(getStatus(iSequence)!=ClpSimplex::basic) {
415          assert(lowerValue>-1.0e20);
416          assert(upperValue<1.0e20);
417        }
418        switch(getStatus(iSequence)) {
419         
420        case ClpSimplex::basic:
421          break;
422        case ClpSimplex::isFree:
423        case ClpSimplex::superBasic:
424          break;
425        case ClpSimplex::atUpperBound:
426          assert (fabs(value-upperValue)<=primalTolerance_) ;
427          break;
428        case ClpSimplex::atLowerBound:
429          assert (fabs(value-lowerValue)<=primalTolerance_) ;
430          break;
431        }
432      }
433    }
434    if(numberIterations_==debugIteration) {
435      printf("dodgy iteration coming up\n");
436    }
437#endif
438    // choose row to go out
439    dualRow();
440    if (pivotRow_>=0) {
441      // we found a pivot row
442      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
443        <<pivotRow_
444        <<CoinMessageEol;
445      // check accuracy of weights
446      dualRowPivot_->checkAccuracy();
447      // get sign for finding row of tableau
448      rowArray_[0]->insert(pivotRow_,directionOut_);
449      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
450      // put row of tableau in rowArray[0] and columnArray[0]
451      matrix_->transposeTimes(this,-1.0,
452                              rowArray_[0],columnArray_[1],columnArray_[0]);
453      // rowArray has pi equivalent
454      // do ratio test
455      dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
456                 rowArray_[3]);
457      if (sequenceIn_>=0) {
458        // normal iteration
459        // update the incoming column
460        unpack(rowArray_[1]);
461        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
462        // and update dual weights (can do in parallel - with extra array)
463        dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
464                                     rowArray_[1]);
465        // see if update stable
466        double btranAlpha = -alpha_*directionOut_; // for check
467        alpha_=(*rowArray_[1])[pivotRow_];
468#ifdef CLP_DEBUG
469        if ((handler_->logLevel()&32))
470          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
471#endif
472        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
473            fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
474          handler_->message(CLP_DUAL_CHECK,messages_)
475            <<btranAlpha
476            <<alpha_
477            <<CoinMessageEol;
478          dualRowPivot_->unrollWeights();
479          if (factorization_->pivots()) {
480            problemStatus_=-2; // factorize now
481            rowArray_[0]->clear();
482            rowArray_[1]->clear();
483            columnArray_[0]->clear();
484            returnCode=-2;
485            break;
486          } else {
487            // take on more relaxed criterion
488            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
489                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
490              // need to reject something
491              char x = isColumn(sequenceOut_) ? 'C' :'R';
492              handler_->message(CLP_SIMPLEX_FLAG,messages_)
493                <<x<<sequenceWithin(sequenceOut_)
494                <<CoinMessageEol;
495              setFlagged(sequenceOut_);
496              lastBadIteration_ = numberIterations_; // say be more cautious
497              rowArray_[0]->clear();
498              rowArray_[1]->clear();
499              columnArray_[0]->clear();
500              continue;
501            }
502          }
503        }
504        // update duals BEFORE replaceColumn so can do updateColumn
505        double objectiveChange=0.0;
506        // do duals first as variables may flip bounds
507        // rowArray_[0] and columnArray_[0] may have flips
508        // so use rowArray_[3] for work array from here on
509        int nswapped = 
510          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
511                            objectiveChange);
512        // which will change basic solution
513        if (nswapped) {
514#ifdef CLP_DEBUG
515          if ((handler_->logLevel()&16))
516            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
517                   dualOut_,valueOut_,lowerOut_,upperOut_);
518          double oldOut=dualOut_;
519#endif
520          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
521          dualRowPivot_->updatePrimalSolution(rowArray_[2],
522                                              1.0,objectiveChange);
523         
524          // recompute dualOut_
525          valueOut_ = solution_[sequenceOut_];
526          if (directionOut_<0) {
527            dualOut_ = valueOut_ - upperOut_;
528          } else {
529            dualOut_ = lowerOut_ - valueOut_;
530          }
531#ifdef CLP_DEBUG
532          if ((handler_->logLevel()&16))
533            printf("%g\n",dualOut_);
534          assert(dualOut_<=oldOut);
535#endif
536          if(dualOut_<0.0&&factorization_->pivots()) {
537            // going backwards - factorize
538            dualRowPivot_->unrollWeights();
539            problemStatus_=-2; // factorize now
540            returnCode=-2;
541            break;
542          }
543        }
544        // amount primal will move
545        double movement = -dualOut_*directionOut_/alpha_;
546        // so objective should increase by fabs(dj)*movement
547        // but we already have objective change - so check will be good
548        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
549#ifdef CLP_DEBUG
550          if (handler_->logLevel()&32)
551            printf("movement %g, swap change %g, rest %g  * %g\n",
552                   objectiveChange+fabs(movement*dualIn_),
553                   objectiveChange,movement,dualIn_);
554#endif
555          if(factorization_->pivots()>5) {
556            // going backwards - factorize
557            dualRowPivot_->unrollWeights();
558            problemStatus_=-2; // factorize now
559            returnCode=-2;
560            break;
561          }
562        }
563        // if stable replace in basis
564        int updateStatus = factorization_->replaceColumn(rowArray_[2],
565                                                         pivotRow_,
566                                                         alpha_);
567        if (updateStatus==1) {
568          // slight error
569          if (factorization_->pivots()>5) {
570            problemStatus_=-2; // factorize now
571            returnCode=-3;
572          }
573        } else if (updateStatus==2) {
574          // major error
575          dualRowPivot_->unrollWeights();
576          // later we may need to unwind more e.g. fake bounds
577          if (factorization_->pivots()) {
578            problemStatus_=-2; // factorize now
579            returnCode=-2;
580            break;
581          } else {
582            // need to reject something
583            char x = isColumn(sequenceOut_) ? 'C' :'R';
584            handler_->message(CLP_SIMPLEX_FLAG,messages_)
585              <<x<<sequenceWithin(sequenceOut_)
586              <<CoinMessageEol;
587            setFlagged(sequenceOut_);
588            lastBadIteration_ = numberIterations_; // say be more cautious
589            rowArray_[0]->clear();
590            rowArray_[1]->clear();
591            columnArray_[0]->clear();
592            continue;
593          }
594        } else if (updateStatus==3) {
595          // out of memory
596          // increase space if not many iterations
597          if (factorization_->pivots()<
598              0.5*factorization_->maximumPivots()&&
599              factorization_->pivots()<200)
600            factorization_->areaFactor(
601                                       factorization_->areaFactor() * 1.1);
602          problemStatus_=-2; // factorize now
603        } 
604        // update primal solution
605        if (theta_<0.0) {
606#ifdef CLP_DEBUG
607          if (handler_->logLevel()&32)
608            printf("negative theta %g\n",theta_);
609#endif
610          theta_=0.0;
611        }
612        // do actual flips
613        flipBounds(rowArray_[0],columnArray_[0],theta_);
614        dualRowPivot_->updatePrimalSolution(rowArray_[1],
615                                            movement,
616                                            objectiveChange);
617#ifdef CLP_DEBUG
618        double oldobj=objectiveValue_;
619#endif
620        int whatNext=housekeeping(objectiveChange);
621        // and set bounds correctly
622        originalBound(sequenceIn_); 
623        changeBound(sequenceOut_);
624#ifdef CLP_DEBUG
625        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
626          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
627#endif
628        if (whatNext==1) {
629          problemStatus_ =-2; // refactorize
630        } else if (whatNext==2) {
631          // maximum iterations or equivalent
632          problemStatus_= 3;
633          returnCode=3;
634          break;
635        }
636      } else {
637        // no incoming column is valid
638#ifdef CLP_DEBUG
639        if (handler_->logLevel()&32)
640          printf("** no column pivot\n");
641#endif
642        if (factorization_->pivots()<5) {
643          problemStatus_=-4; //say looks infeasible
644          // create ray anyway
645          delete [] ray_;
646          ray_ = new double [ numberRows_];
647          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
648        }
649        rowArray_[0]->clear();
650        columnArray_[0]->clear();
651        returnCode=1;
652        break;
653      }
654    } else {
655      // no pivot row
656#ifdef CLP_DEBUG
657      if (handler_->logLevel()&32)
658        printf("** no row pivot\n");
659#endif
660      if (!factorization_->pivots()) {
661        // may have crept through - so may be optimal
662        //problemStatus_=-5; //say looks unbounded
663        problemStatus_=0;
664        // check any flagged variables
665        int iRow;
666        for (iRow=0;iRow<numberRows_;iRow++) {
667          int iPivot=pivotVariable_[iRow];
668          if (flagged(iPivot))
669            break;
670        }
671        if (iRow<numberRows_) {
672#ifdef CLP_DEBUG
673          std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
674#endif
675          problemStatus_=-4; //say looks infeasible
676          // create ray anyway
677          delete [] ray_;
678          ray_ = new double [ numberRows_];
679          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
680        }
681      }
682      returnCode=0;
683      break;
684    }
685  }
686  return returnCode;
687}
688/* The duals are updated by the given arrays.
689   Returns number of infeasibilities.
690   rowArray and columnarray will have flipped
691   The output vector has movement (row length array) */
692int
693ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
694                                  CoinIndexedVector * columnArray,
695                                  CoinIndexedVector * outputArray,
696                                  double theta,
697                                  double & objectiveChange)
698{
699 
700  outputArray->clear();
701 
702  double * work;
703  int number;
704  int * which;
705 
706  int numberInfeasibilities=0;
707  int numberRowInfeasibilities=0;
708 
709  // see whether we will be doing full recompute
710  bool fullRecompute= (rowArray->getNumElements()==numberRows_&&
711                       columnArray->getNumElements()==numberColumns_);
712  int numberAtFake=0;
713 
714  // use a tighter tolerance except for all being okay
715  double tolerance = dualTolerance_;
716 
717  double changeObj=0.0;
718 
719  int iSection;
720 
721  for (iSection=0;iSection<2;iSection++) {
722    int i;
723    double * solution = solutionRegion(iSection);
724    double * reducedCost = djRegion(iSection);
725    double * lower = lowerRegion(iSection);
726    double * upper = upperRegion(iSection);
727    double * cost = costRegion(iSection);
728    int addSequence;
729    if (!iSection) {
730      addSequence = numberColumns_;
731      work = rowArray->denseVector();
732      number = rowArray->getNumElements();
733      which = rowArray->getIndices();
734    } else {
735      // set number of infeasibilities in row array
736      addSequence=0;
737      numberRowInfeasibilities=numberInfeasibilities;
738      rowArray->setNumElements(numberInfeasibilities);
739      numberInfeasibilities=0;
740      work = columnArray->denseVector();
741      number = columnArray->getNumElements();
742      which = columnArray->getIndices();
743    }
744   
745    for (i=0;i<number;i++) {
746      int iSequence = which[i];
747      double alphaI = work[iSequence];
748      double value = reducedCost[iSequence]-theta*alphaI;
749      work[iSequence]=0.0;
750      reducedCost[iSequence]=value;
751     
752      if (!fixed(iSequence+addSequence)) {
753        double movement=0.0;
754        FakeBound bound = getFakeBound(iSequence+addSequence);
755        Status status = getStatus(iSequence+addSequence);
756
757        switch(status) {
758         
759        case ClpSimplex::basic:
760        case ClpSimplex::superBasic:
761          break;
762        case ClpSimplex::isFree:
763          if (fabs(value)>tolerance) { 
764#ifdef CLP_DEBUG
765            if (handler_->logLevel()&32)
766              printf("%d %d, free has dj of %g, alpha %g\n",
767                     iSection,iSequence,value,alphaI);
768#endif
769          }
770          break;
771        case ClpSimplex::atUpperBound:
772          if (value>tolerance) {
773            // to lower bound (if swap)
774            // put back alpha
775            which[numberInfeasibilities++]=iSequence;
776            work[iSequence]=alphaI;
777            movement = lower[iSequence]-upper[iSequence];
778#ifdef CLP_DEBUG
779            if ((handler_->logLevel()&32))
780              printf("%d %d, new dj %g, alpha %g, movement %g\n",
781                   iSection,iSequence,value,alphaI,movement);
782#endif
783            changeObj += movement*cost[iSequence];
784            if (bound==ClpSimplexDual::bothFake||
785                bound==ClpSimplexDual::lowerFake) 
786              numberAtFake++;
787          } else if (fullRecompute) {
788            // at correct bound
789            if (bound==ClpSimplexDual::bothFake||
790                bound==ClpSimplexDual::upperFake) {
791              // but flip if dj would allow
792              if (bound==ClpSimplexDual::upperFake&&
793                  value>=-tolerance) {
794                movement = lower[iSequence]-upper[iSequence];
795                setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
796                solution[iSequence] = lower[iSequence];
797                changeObj += movement*cost[iSequence];
798              } else {
799                numberAtFake++;
800              }
801            }
802          }
803          break;
804        case ClpSimplex::atLowerBound:
805          if (value<-tolerance) {
806            // to upper bound
807            // put back alpha
808            which[numberInfeasibilities++]=iSequence;
809            work[iSequence]=alphaI;
810            movement = upper[iSequence] - lower[iSequence];
811#ifdef CLP_DEBUG
812            if ((handler_->logLevel()&32))
813              printf("%d %d, new dj %g, alpha %g, movement %g\n",
814                   iSection,iSequence,value,alphaI,movement);
815#endif
816            changeObj += movement*cost[iSequence];
817            if (bound==ClpSimplexDual::bothFake||
818                bound==ClpSimplexDual::upperFake) 
819              numberAtFake++;
820          } else if (fullRecompute) {
821            // at correct bound
822            if (bound==ClpSimplexDual::bothFake||
823                bound==ClpSimplexDual::lowerFake) {
824              // but flip if dj would allow
825              if (bound==ClpSimplexDual::lowerFake&&
826                  value<=tolerance) {
827                movement = upper[iSequence] - lower[iSequence];
828                setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
829                solution[iSequence] = upper[iSequence];
830                changeObj += movement*cost[iSequence];
831              } else {
832                numberAtFake++;
833              }
834            }
835          }
836          break;
837        }
838        if (!fullRecompute) {
839          if (movement) {
840            if (!iSection) {
841              // row (sign ?)
842              outputArray->quickAdd(iSequence,-movement);
843            } else {
844              matrix_->add(this,outputArray,iSequence,movement);
845            }
846          }
847        }
848      }
849    }
850  }
851#ifdef CLP_DEBUG
852  if (fullRecompute&&numberAtFake&&(handler_->logLevel()&16)!=0) 
853    printf("%d fake after full update\n",numberAtFake);
854#endif
855  outputArray->stopQuickAdd();
856  // set number of infeasibilities
857  columnArray->setNumElements(numberInfeasibilities);
858  numberInfeasibilities += numberRowInfeasibilities;
859  if (fullRecompute) {
860    // do actual flips
861    flipBounds(rowArray,columnArray,theta);
862    numberFake_ = numberAtFake;
863  }
864  objectiveChange += changeObj;
865  return numberInfeasibilities;
866}
867/*
868   Chooses dual pivot row
869   Would be faster with separate region to scan
870   and will have this (with square of infeasibility) when steepest
871   For easy problems we can just choose one of the first rows we look at
872*/
873void
874ClpSimplexDual::dualRow()
875{
876  // get pivot row using whichever method it is
877  pivotRow_=dualRowPivot_->pivotRow();
878  if (pivotRow_>=0) {
879    int iPivot=pivotVariable_[pivotRow_];
880    sequenceOut_ = iPivot;
881    if (iPivot>=numberColumns_) {
882      // slack
883      iPivot-=numberColumns_;
884      valueOut_=rowActivityWork_[iPivot];
885      lowerOut_=rowLowerWork_[iPivot];
886      upperOut_=rowUpperWork_[iPivot];
887    } else {
888      // column
889      valueOut_=columnActivityWork_[iPivot];
890      lowerOut_=columnLowerWork_[iPivot];
891      upperOut_=columnUpperWork_[iPivot];
892    }
893    // if we have problems we could try other way and hope we get a
894    // zero pivot?
895    if (valueOut_>upperOut_) {
896      directionOut_ = -1;
897      dualOut_ = valueOut_ - upperOut_;
898    } else {
899      directionOut_ = 1;
900      dualOut_ = lowerOut_ - valueOut_;
901    }
902#ifdef CLP_DEBUG
903    assert(dualOut_>=0.0);
904#endif
905  }
906  return ;
907}
908// Checks if any fake bounds active - if so returns number and modifies
909// dualBound_ and everything.
910// Free variables will be left as free
911// Returns number of bounds changed if >=0
912// Returns -1 if not initialize and no effect
913// Fills in changeVector which can be used to see if unbounded
914// and cost of change vector
915int
916ClpSimplexDual::changeBounds(bool initialize,
917                                 CoinIndexedVector * outputArray,
918                                 double & changeCost)
919{
920  if (!initialize) {
921    int numberInfeasibilities;
922    double newBound;
923    newBound = 5.0*dualBound_;
924    numberInfeasibilities=0;
925    changeCost=0.0;
926    // put back original bounds and then check
927    createRim(3);
928    int iSequence;
929    // bounds will get bigger - just look at ones at bounds
930    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
931      double lowerValue=lower_[iSequence];
932      double upperValue=upper_[iSequence];
933      double value=solution_[iSequence];
934      setFakeBound(iSequence,ClpSimplexDual::noFake);
935      switch(getStatus(iSequence)) {
936
937      case ClpSimplex::basic:
938        break;
939      case ClpSimplex::isFree:
940      case ClpSimplex::superBasic:
941        break;
942      case ClpSimplex::atUpperBound:
943        if (fabs(value-upperValue)>primalTolerance_) 
944          numberInfeasibilities++;
945        break;
946      case ClpSimplex::atLowerBound:
947        if (fabs(value-lowerValue)>primalTolerance_) 
948          numberInfeasibilities++;
949        break;
950      }
951    }
952    if (numberInfeasibilities) {
953      int iSequence;
954      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
955        double lowerValue=lower_[iSequence];
956        double upperValue=upper_[iSequence];
957        double newLowerValue;
958        double newUpperValue;
959        Status status = getStatus(iSequence);
960        if (status==ClpSimplex::atUpperBound||
961            status==ClpSimplex::atLowerBound) {
962          double value = solution_[iSequence];
963          if (value-lowerValue<=upperValue-value) {
964            newLowerValue = max(lowerValue,value-0.666667*newBound);
965            newUpperValue = min(upperValue,newLowerValue+newBound);
966          } else {
967            newUpperValue = min(upperValue,value+0.666667*newBound);
968            newLowerValue = max(lowerValue,newUpperValue-newBound);
969          }
970          lower_[iSequence]=newLowerValue;
971          upper_[iSequence]=newUpperValue;
972          if (newLowerValue > lowerValue) {
973            if (newUpperValue < upperValue) 
974              setFakeBound(iSequence,ClpSimplexDual::bothFake);
975            else
976              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
977          } else {
978            if (newUpperValue < upperValue) 
979              setFakeBound(iSequence,ClpSimplexDual::upperFake);
980          }
981          if (status==ClpSimplex::atUpperBound)
982            solution_[iSequence] = newUpperValue;
983          else
984            solution_[iSequence] = newLowerValue;
985          double movement = solution_[iSequence] - value;
986          if (movement&&outputArray) {
987            if (iSequence>=numberColumns_) {
988              outputArray->quickAdd(iSequence,-movement);
989              changeCost += movement*cost_[iSequence];
990            } else {
991              matrix_->add(this,outputArray,iSequence,movement);
992              changeCost += movement*cost_[iSequence];
993            }
994          }
995        }
996      }
997      dualBound_ = newBound;
998      if (outputArray)
999        outputArray->stopQuickAdd();
1000    } else {
1001      numberInfeasibilities=-1;
1002    }
1003    return numberInfeasibilities;
1004  } else {
1005    int iSequence;
1006     
1007    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1008      Status status = getStatus(iSequence);
1009      if (status==ClpSimplex::atUpperBound||
1010          status==ClpSimplex::atLowerBound) {
1011        double lowerValue=lower_[iSequence];
1012        double upperValue=upper_[iSequence];
1013        double value = solution_[iSequence];
1014        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1015          if (lowerValue-value>-0.5*dualBound_||
1016              upperValue-value<0.5*dualBound_) {
1017            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1018              if (upperValue > lowerValue + dualBound_) {
1019                upper_[iSequence]=lowerValue+dualBound_;
1020                setFakeBound(iSequence,ClpSimplexDual::upperFake);
1021              }
1022            } else {
1023              if (lowerValue < upperValue - dualBound_) {
1024                lower_[iSequence]=upperValue-dualBound_;
1025                setFakeBound(iSequence,ClpSimplexDual::lowerFake);
1026              }
1027            }
1028          } else {
1029            lower_[iSequence]=-0.5*dualBound_;
1030            upper_[iSequence]= 0.5*dualBound_;
1031            setFakeBound(iSequence,ClpSimplexDual::bothFake);
1032          }
1033        }
1034      }
1035    }
1036    return 1;
1037  }
1038}
1039/*
1040   Row array has row part of pivot row (as duals so sign may be switched)
1041   Column array has column part.
1042   This chooses pivot column.
1043   Spare array will be needed when we start getting clever.
1044   We will check for basic so spare array will never overflow.
1045   If necessary will modify costs
1046*/
1047void
1048ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
1049                       CoinIndexedVector * columnArray,
1050                       CoinIndexedVector * spareArray,
1051                       CoinIndexedVector * spareArray2)
1052{
1053  double * work;
1054  int number;
1055  int * which;
1056  double * reducedCost;
1057
1058  int iSection;
1059
1060  sequenceIn_=-1;
1061  int numberPossiblySwapped=0;
1062  int numberRemaining=0;
1063 
1064  double totalThru=0.0; // for when variables flip
1065  double acceptablePivot=1.0e-7;
1066  if (factorization_->pivots())
1067    acceptablePivot=1.0e-5; // if we have iterated be more strict
1068  double bestEverPivot=acceptablePivot;
1069  int lastSequence = -1;
1070  double lastPivot=0.0;
1071  double upperTheta;
1072  double newTolerance = dualTolerance_;
1073  // will we need to increase tolerance
1074  bool thisIncrease=false;
1075  // If we think we need to modify costs (not if something from broad sweep)
1076  bool modifyCosts=false;
1077  // Increase in objective due to swapping bounds (may be negative)
1078  double increaseInObjective=0.0;
1079
1080  // use spareArrays to put ones looked at in
1081  // we are going to flip flop between
1082  int iFlip = 0;
1083  // Possible list of pivots
1084  int interesting[2];
1085  // where possible swapped ones are
1086  int swapped[2];
1087  // for zeroing out arrays after
1088  int marker[2][2];
1089  // pivot elements
1090  double * array[2], * spare, * spare2;
1091  // indices
1092  int * indices[2], * index, * index2;
1093  spareArray->clear();
1094  spareArray2->clear();
1095  array[0] = spareArray->denseVector();
1096  indices[0] = spareArray->getIndices();
1097  spare = array[0];
1098  index = indices[0];
1099  array[1] = spareArray2->denseVector();
1100  indices[1] = spareArray2->getIndices();
1101  int i;
1102  double * lower;
1103  double * upper;
1104
1105  // initialize lists
1106  for (i=0;i<2;i++) {
1107    interesting[i]=0;
1108    swapped[i]=numberColumns_;
1109    marker[i][0]=0;
1110    marker[i][1]=numberColumns_;
1111  }
1112
1113  /*
1114    First we get a list of possible pivots.  We can also see if the
1115    problem looks infeasible or whether we want to pivot in free variable.
1116    This may make objective go backwards but can only happen a finite
1117    number of times and I do want free variables basic.
1118
1119    Then we flip back and forth.  At the start of each iteration
1120    interesting[iFlip] should have possible candidates and swapped[iFlip]
1121    will have pivots if we decide to take a previous pivot.
1122    At end of each iteration interesting[1-iFlip] should have
1123    candidates if we go through this theta and swapped[1-iFlip]
1124    pivots if we don't go through.
1125
1126    At first we increase theta and see what happens.  We start
1127    theta at a reasonable guess.  If in right area then we do bit by bit.
1128
1129   */
1130
1131  // do first pass to get possibles
1132  // We can also see if infeasible or pivoting on free
1133  double tentativeTheta = 1.0e22;
1134  upperTheta = 1.0e31;
1135  double freePivot = acceptablePivot;
1136  for (iSection=0;iSection<2;iSection++) {
1137
1138    int addSequence;
1139
1140    if (!iSection) {
1141      lower = rowLowerWork_;
1142      upper = rowUpperWork_;
1143      work = rowArray->denseVector();
1144      number = rowArray->getNumElements();
1145      which = rowArray->getIndices();
1146      reducedCost = rowReducedCost_;
1147      addSequence = numberColumns_;
1148    } else {
1149      lower = columnLowerWork_;
1150      upper = columnUpperWork_;
1151      work = columnArray->denseVector();
1152      number = columnArray->getNumElements();
1153      which = columnArray->getIndices();
1154      reducedCost = reducedCostWork_;
1155      addSequence = 0;
1156    }
1157   
1158    for (i=0;i<number;i++) {
1159      int iSequence = which[i];
1160      double alpha = work[iSequence];
1161      if (fixed(iSequence+addSequence)||!alpha) 
1162        continue; // skip fixed ones or (zeroed out)
1163      double oldValue = reducedCost[iSequence];
1164      double value = oldValue-tentativeTheta*alpha;
1165      int keep = 0;
1166
1167      switch(getStatus(iSequence+addSequence)) {
1168         
1169      case ClpSimplex::basic:
1170        break;
1171      case ClpSimplex::isFree:
1172      case ClpSimplex::superBasic:
1173        if (oldValue>dualTolerance_) {
1174          if (value<-newTolerance) 
1175            keep = 2;
1176        } else if (oldValue<-dualTolerance_) {
1177          if (value>newTolerance) 
1178            keep = 2;
1179        } else {
1180          if (alpha>=acceptablePivot) 
1181            keep = 2;
1182          else if (-alpha>=acceptablePivot) 
1183            keep = 2;
1184        }
1185        break;
1186      case ClpSimplex::atUpperBound:
1187        assert (oldValue<=dualTolerance_*1.0001);
1188        if (value>newTolerance) {
1189          keep = 1;
1190          value = oldValue-upperTheta*alpha;
1191          if (value>newTolerance && -alpha>=acceptablePivot) 
1192            upperTheta = (oldValue-newTolerance)/alpha;
1193        }
1194        break;
1195      case ClpSimplex::atLowerBound:
1196        assert (oldValue>=-dualTolerance_*1.0001);
1197        if (value<-newTolerance) {
1198          keep = 1;
1199          value = oldValue-upperTheta*alpha;
1200          if (value<-newTolerance && alpha>=acceptablePivot) 
1201            upperTheta = (oldValue+newTolerance)/alpha;
1202        }
1203        break;
1204      }
1205      if (keep) {
1206        if (keep==2) {
1207          // free - choose largest
1208          if (fabs(alpha)>freePivot) {
1209            freePivot=fabs(alpha);
1210            sequenceIn_ = iSequence + addSequence;
1211            theta_=oldValue/alpha;
1212          }
1213        } else {
1214          // add to list
1215          spare[numberRemaining]=alpha;
1216          index[numberRemaining++]=iSequence+addSequence;
1217        }
1218      }
1219    }
1220  }
1221  interesting[0]=numberRemaining;
1222  marker[0][0] = numberRemaining;
1223
1224  if (!numberRemaining)
1225    return; // Looks infeasible
1226
1227  if (sequenceIn_>=0) {
1228    // free variable - always choose
1229  } else {
1230
1231    theta_=1.0e50;
1232    // now flip flop between spare arrays until reasonable theta
1233    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1234   
1235    // loops increasing tentative theta until can't go through
1236   
1237    while (tentativeTheta < 1.0e22) {
1238      double thruThis = 0.0;
1239     
1240      double bestPivot=acceptablePivot;
1241      int bestSequence=-1;
1242     
1243      numberPossiblySwapped = numberColumns_;
1244      numberRemaining = 0;
1245
1246      upperTheta = 1.0e50;
1247
1248      spare = array[iFlip];
1249      index = indices[iFlip];
1250      spare2 = array[1-iFlip];
1251      index2 = indices[1-iFlip];
1252
1253      // try 3 different ways
1254      // 1 bias increase by ones with slightly wrong djs
1255      // 2 bias by all
1256      // 3 bias by all - tolerance (doesn't seem very good)
1257#define TRYBIAS 1
1258
1259
1260      double increaseInThis=0.0; //objective increase in this loop
1261     
1262      for (i=0;i<interesting[iFlip];i++) {
1263        int iSequence = index[i];
1264        double alpha = spare[i];
1265        double oldValue = dj_[iSequence];
1266        double value = oldValue-tentativeTheta*alpha;
1267
1268        if (alpha < 0.0) {
1269          //at upper bound
1270          if (value>newTolerance) {
1271            double range = upper_[iSequence] - lower_[iSequence];
1272            thruThis -= range*alpha;
1273#if TRYBIAS==1
1274            if (oldValue>0.0)
1275              increaseInThis -= oldValue*range;
1276#elif TRYBIAS==2
1277            increaseInThis -= oldValue*range;
1278#else
1279            increaseInThis -= (oldValue+dualTolerance_)*range;
1280#endif
1281            // goes on swapped list (also means candidates if too many)
1282            spare2[--numberPossiblySwapped]=alpha;
1283            index2[numberPossiblySwapped]=iSequence;
1284            if (fabs(alpha)>bestPivot) {
1285              bestPivot=fabs(alpha);
1286              bestSequence=numberPossiblySwapped;
1287            }
1288          } else {
1289            value = oldValue-upperTheta*alpha;
1290            if (value>newTolerance && -alpha>=acceptablePivot) 
1291              upperTheta = (oldValue-newTolerance)/alpha;
1292            spare2[numberRemaining]=alpha;
1293            index2[numberRemaining++]=iSequence;
1294          }
1295        } else {
1296          // at lower bound
1297          if (value<-newTolerance) {
1298            double range = upper_[iSequence] - lower_[iSequence];
1299            thruThis += range*alpha;
1300            //?? is this correct - and should we look at good ones
1301#if TRYBIAS==1
1302            if (oldValue<0.0)
1303              increaseInThis += oldValue*range;
1304#elif TRYBIAS==2
1305            increaseInThis += oldValue*range;
1306#else
1307            increaseInThis += (oldValue-dualTolerance_)*range;
1308#endif
1309            // goes on swapped list (also means candidates if too many)
1310            spare2[--numberPossiblySwapped]=alpha;
1311            index2[numberPossiblySwapped]=iSequence;
1312            if (fabs(alpha)>bestPivot) {
1313              bestPivot=fabs(alpha);
1314              bestSequence=numberPossiblySwapped;
1315            }
1316          } else {
1317            value = oldValue-upperTheta*alpha;
1318            if (value<-newTolerance && alpha>=acceptablePivot) 
1319              upperTheta = (oldValue+newTolerance)/alpha;
1320            spare2[numberRemaining]=alpha;
1321            index2[numberRemaining++]=iSequence;
1322          }
1323        }
1324      }
1325      swapped[1-iFlip]=numberPossiblySwapped;
1326      interesting[1-iFlip]=numberRemaining;
1327      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1328      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1329     
1330      if (totalThru+thruThis>=fabs(dualOut_)||
1331          increaseInObjective+increaseInThis<0.0) {
1332        // We should be pivoting in this batch
1333        // so compress down to this lot
1334        numberRemaining=0;
1335        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1336          spare[numberRemaining]=spare2[i];
1337          index[numberRemaining++]=index2[i];
1338        }
1339        interesting[iFlip]=numberRemaining;
1340        int iTry;
1341#define MAXTRY 100
1342        // first get ratio with tolerance
1343        for (iTry=0;iTry<MAXTRY;iTry++) {
1344         
1345          upperTheta=1.0e50;
1346          numberPossiblySwapped = numberColumns_;
1347          numberRemaining = 0;
1348
1349          increaseInThis=0.0; //objective increase in this loop
1350
1351          thruThis=0.0;
1352
1353          spare = array[iFlip];
1354          index = indices[iFlip];
1355          spare2 = array[1-iFlip];
1356          index2 = indices[1-iFlip];
1357     
1358          for (i=0;i<interesting[iFlip];i++) {
1359            int iSequence=index[i];
1360            double alpha=spare[i];
1361            double oldValue = dj_[iSequence];
1362            double value = oldValue-upperTheta*alpha;
1363           
1364            if (alpha < 0.0) {
1365              //at upper bound
1366              if (value>newTolerance) {
1367                if (-alpha>=acceptablePivot) {
1368                  upperTheta = (oldValue-newTolerance)/alpha;
1369                }
1370              }
1371            } else {
1372              // at lower bound
1373              if (value<-newTolerance) {
1374                if (alpha>=acceptablePivot) {
1375                  upperTheta = (oldValue+newTolerance)/alpha;
1376                }
1377              }
1378            }
1379          }
1380          bestPivot=acceptablePivot;
1381          sequenceIn_=-1;
1382          // now choose largest and sum all ones which will go through
1383#define MINIMUMTHETA 1.0e-12
1384          for (i=0;i<interesting[iFlip];i++) {
1385            int iSequence=index[i];
1386            double alpha=spare[i];
1387            double value = dj_[iSequence]-upperTheta*alpha;
1388            double badDj=0.0;
1389           
1390            bool addToSwapped=false;
1391           
1392            if (alpha < 0.0) {
1393              //at upper bound
1394              if (value>=0.0) { 
1395                addToSwapped=true;
1396#if TRYBIAS==1
1397                badDj = -max(dj_[iSequence],0.0);
1398#elif TRYBIAS==2
1399                badDj = -dj_[iSequence];
1400#else
1401                badDj = -dj_[iSequence]-dualTolerance_;
1402#endif
1403              }
1404            } else {
1405              // at lower bound
1406              if (value<=0.0) {
1407                addToSwapped=true;
1408#if TRYBIAS==1
1409                badDj = min(dj_[iSequence],0.0);
1410#elif TRYBIAS==2
1411                badDj = dj_[iSequence];
1412#else
1413                badDj = dj_[iSequence]-dualTolerance_;
1414#endif
1415              }
1416            }
1417            if (!addToSwapped) {
1418              // add to list of remaining
1419              spare2[numberRemaining]=alpha;
1420              index2[numberRemaining++]=iSequence;
1421            } else {
1422              // add to list of swapped
1423              spare2[--numberPossiblySwapped]=alpha;
1424              index2[numberPossiblySwapped]=iSequence;
1425              // select if largest pivot
1426              if (fabs(alpha)>bestPivot) {
1427                sequenceIn_ = numberPossiblySwapped;
1428                bestPivot =  fabs(alpha);
1429                theta_ = dj_[iSequence]/alpha;
1430              }
1431              double range = upper[iSequence] - lower[iSequence];
1432              thruThis += range*fabs(alpha);
1433              increaseInThis += badDj*range;
1434            }
1435          }
1436          swapped[1-iFlip]=numberPossiblySwapped;
1437          interesting[1-iFlip]=numberRemaining;
1438          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1439          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1440          // If we stop now this will be increase in objective (I think)
1441          double increase = (fabs(dualOut_)-totalThru)*theta_;
1442          increase += increaseInObjective;
1443          if (theta_<0.0)
1444            thruThis += fabs(dualOut_); // force using this one
1445          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
1446            // back
1447            // We may need to be more careful - we could do by
1448            // switch so we always do fine grained?
1449            bestPivot=0.0;
1450          } else {
1451            // add in
1452            totalThru += thruThis;
1453            increaseInObjective += increaseInThis;
1454          }
1455          if (bestPivot<0.1*bestEverPivot&&
1456              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1457            // back to previous one
1458            sequenceIn_=lastSequence;
1459            // swap regions
1460            iFlip = 1-iFlip;
1461            break;
1462          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
1463            if (lastPivot>acceptablePivot) {
1464              // back to previous one
1465              sequenceIn_=lastSequence;
1466              // swap regions
1467              iFlip = 1-iFlip;
1468            } else {
1469              // can only get here if all pivots too small
1470            }
1471            break;
1472          } else if (totalThru>=fabs(dualOut_)) {
1473            modifyCosts=true; // fine grain - we can modify costs
1474            break; // no point trying another loop
1475          } else {
1476            lastSequence=sequenceIn_;
1477            if (bestPivot>bestEverPivot)
1478              bestEverPivot=bestPivot;
1479            iFlip = 1 -iFlip;
1480            modifyCosts=true; // fine grain - we can modify costs
1481        }
1482        }
1483        if (iTry==MAXTRY)
1484          iFlip = 1-iFlip; // flip back
1485        break;
1486      } else {
1487        // skip this lot
1488        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
1489          bestEverPivot=bestPivot;
1490          lastSequence=bestSequence;
1491        } else {
1492          // keep old swapped
1493          memcpy(array[1-iFlip]+swapped[iFlip],
1494                 array[iFlip]+swapped[iFlip],
1495                 (numberColumns_-swapped[iFlip])*sizeof(double));
1496          memcpy(indices[1-iFlip]+swapped[iFlip],
1497                 indices[iFlip]+swapped[iFlip],
1498                 (numberColumns_-swapped[iFlip])*sizeof(int));
1499          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
1500          swapped[1-iFlip]=swapped[iFlip];
1501        }
1502        increaseInObjective += increaseInThis;
1503        iFlip = 1 - iFlip; // swap regions
1504        tentativeTheta = 2.0*upperTheta;
1505        totalThru += thruThis;
1506      }
1507    }
1508   
1509    // can get here without sequenceIn_ set but with lastSequence
1510    if (sequenceIn_<0&&lastSequence>=0) {
1511      // back to previous one
1512      sequenceIn_=lastSequence;
1513      // swap regions
1514      iFlip = 1-iFlip;
1515    }
1516   
1517    if (sequenceIn_>=0) {
1518      // at this stage sequenceIn_ is just pointer into index array
1519      // flip just so we can use iFlip
1520      iFlip = 1 -iFlip;
1521      spare = array[iFlip];
1522      index = indices[iFlip];
1523      double oldValue;
1524      double alpha = spare[sequenceIn_];
1525      sequenceIn_ = indices[iFlip][sequenceIn_];
1526      oldValue = dj_[sequenceIn_];
1527      theta_ = oldValue/alpha;
1528      if (theta_<MINIMUMTHETA) {
1529        // can't pivot to zero
1530        if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
1531          theta_=MINIMUMTHETA;
1532        } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
1533          theta_=MINIMUMTHETA;
1534          thisIncrease=true;
1535        } else {
1536          theta_=(oldValue+newTolerance)/alpha;
1537          assert(theta_>=0.0);
1538          thisIncrease=true;
1539        } 
1540      }
1541      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
1542      if (modifyCosts) {
1543        int i;
1544        double * workRow = rowArray->denseVector();
1545        double * workColumn = columnArray->denseVector();
1546        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
1547          int iSequence=index[i];
1548          double alpha;
1549          if (iSequence>=numberColumns_)
1550            alpha=workRow[iSequence-numberColumns_];
1551          else
1552            alpha=workColumn[iSequence];
1553          double value = dj_[iSequence]-theta_*alpha;
1554           
1555          // can't be free here
1556         
1557          if (alpha < 0.0) {
1558            //at upper bound
1559            if (value>dualTolerance_) {
1560              thisIncrease=true;
1561#define MODIFYCOST 2
1562#if MODIFYCOST
1563              // modify cost to hit new tolerance
1564              double modification = alpha*theta_-dj_[iSequence]
1565                +newTolerance;
1566              //modification = min(modification,dualTolerance_);
1567              //assert (fabs(modification)<1.0e-7);
1568              dj_[iSequence] += modification;
1569              cost_[iSequence] +=  modification;
1570#endif
1571            }
1572          } else {
1573            // at lower bound
1574            if (-value>dualTolerance_) {
1575              thisIncrease=true;
1576#if MODIFYCOST
1577              // modify cost to hit new tolerance
1578              double modification = alpha*theta_-dj_[iSequence]
1579                -newTolerance;
1580              //modification = max(modification,-dualTolerance_);
1581              //assert (fabs(modification)<1.0e-7);
1582              dj_[iSequence] += modification;
1583              cost_[iSequence] +=  modification;
1584#endif
1585            }
1586          }
1587        }
1588      }
1589    }
1590  }
1591
1592  if (sequenceIn_>=0) {
1593    if (sequenceIn_>=numberColumns_) {
1594      //slack
1595      alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
1596    } else {
1597      // column
1598      alpha_ = columnArray->denseVector()[sequenceIn_];
1599    }
1600    lowerIn_ = lower_[sequenceIn_];
1601    upperIn_ = upper_[sequenceIn_];
1602    valueIn_ = solution_[sequenceIn_];
1603    dualIn_ = dj_[sequenceIn_];
1604
1605    if (numberTimesOptimal_) {
1606      // can we adjust cost back closer to original
1607      //*** add coding
1608    }
1609#if MODIFYCOST>1   
1610    // modify cost to hit zero exactly
1611    // so (dualIn_+modification)==theta_*alpha_
1612    double modification = theta_*alpha_-dualIn_;
1613    dualIn_ += modification;
1614    dj_[sequenceIn_]=dualIn_;
1615    cost_[sequenceIn_] += modification;
1616    //assert (fabs(modification)<1.0e-6);
1617#ifdef CLP_DEBUG
1618    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
1619      printf("exact %d new cost %g, change %g\n",sequenceIn_,
1620             cost_[sequenceIn_],modification);
1621#endif
1622#endif
1623   
1624    if (alpha_<0.0) {
1625      // as if from upper bound
1626      directionIn_=-1;
1627      upperIn_=valueIn_;
1628    } else {
1629      // as if from lower bound
1630      directionIn_=1;
1631      lowerIn_=valueIn_;
1632    }
1633  }
1634  if (thisIncrease) {
1635    newTolerance = dualTolerance_+1.0e-4*dblParam_[ClpDualTolerance];
1636  }
1637
1638  // clear arrays
1639
1640  for (i=0;i<2;i++) {
1641    memset(array[i],0,marker[i][0]*sizeof(double));
1642    memset(array[i]+marker[i][1],0,
1643           (numberColumns_-marker[i][1])*sizeof(double));
1644  }
1645}
1646/* Checks if tentative optimal actually means unbounded
1647   Returns -3 if not, 2 if is unbounded */
1648int 
1649ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
1650                                   CoinIndexedVector * spare,
1651                                   double changeCost)
1652{
1653  int status=2; // say unbounded
1654  factorization_->updateColumn(spare,ray);
1655  // get reduced cost
1656  int i;
1657  int number=ray->getNumElements();
1658  int * index = ray->getIndices();
1659  double * array = ray->denseVector();
1660  for (i=0;i<number;i++) {
1661    int iRow=index[i];
1662    int iPivot=pivotVariable_[iRow];
1663    changeCost -= cost(iPivot)*array[iRow];
1664  }
1665  double way;
1666  if (changeCost>0.0) {
1667    //try going down
1668    way=1.0;
1669  } else if (changeCost<0.0) {
1670    //try going up
1671    way=-1.0;
1672  } else {
1673#ifdef CLP_DEBUG
1674    printf("can't decide on up or down\n");
1675#endif
1676    way=0.0;
1677    status=-3;
1678  }
1679  double movement=1.0e10*way; // some largish number
1680  double zeroTolerance = 1.0e-14*dualBound_;
1681  for (i=0;i<number;i++) {
1682    int iRow=index[i];
1683    int iPivot=pivotVariable_[iRow];
1684    double arrayValue = array[iRow];
1685    if (fabs(arrayValue)<zeroTolerance)
1686      arrayValue=0.0;
1687    double newValue=solution(iPivot)+movement*arrayValue;
1688    if (newValue>upper(iPivot)+primalTolerance_||
1689        newValue<lower(iPivot)-primalTolerance_)
1690      status=-3; // not unbounded
1691  }
1692  if (status==2) {
1693    // create ray
1694    delete [] ray_;
1695    ray_ = new double [numberColumns_];
1696    CoinFillN(ray_,numberColumns_,0.0);
1697    for (i=0;i<number;i++) {
1698      int iRow=index[i];
1699      int iPivot=pivotVariable_[iRow];
1700      double arrayValue = array[iRow];
1701      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
1702        ray_[iPivot] = way* array[iRow];
1703    }
1704  }
1705  ray->clear();
1706  return status;
1707}
1708/* Checks if finished.  Updates status */
1709void 
1710ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type)
1711{
1712  if (type==2) {
1713    // trouble - restore solution
1714    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1715    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1716           numberRows_*sizeof(double));
1717    memcpy(columnActivityWork_,savedSolution_ ,
1718           numberColumns_*sizeof(double));
1719    forceFactorization_=1; // a bit drastic but ..
1720    changeMade_++; // say something changed
1721  }
1722  int tentativeStatus = problemStatus_;
1723  double changeCost;
1724
1725  if (problemStatus_>-3) {
1726    // factorize
1727    // later on we will need to recover from singularities
1728    // also we could skip if first time
1729    // save dual weights
1730    dualRowPivot_->saveWeights(this,1);
1731    if (type) {
1732      // is factorization okay?
1733      if (internalFactorize(1)) {
1734        // no - restore previous basis
1735        assert (type==1);
1736        changeMade_++; // say something changed
1737        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
1738        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
1739               numberRows_*sizeof(double));
1740        memcpy(columnActivityWork_,savedSolution_ ,
1741               numberColumns_*sizeof(double));
1742        // get correct bounds on all variables
1743        double dummyChangeCost=0.0;
1744        changeBounds(true,rowArray_[2],dummyChangeCost);
1745        // throw away change
1746        rowArray_[2]->clear();
1747        forceFactorization_=1; // a bit drastic but ..
1748        type = 2;
1749        assert (internalFactorize(1)==0);
1750      }
1751    }
1752    problemStatus_=-3;
1753  }
1754  // at this stage status is -3 or -4 if looks infeasible
1755  // get primal and dual solutions
1756  gutsOfSolution(rowActivityWork_,columnActivityWork_);
1757#ifdef CLP_DEBUG
1758  if (!rowScale_&&(handler_->logLevel()&32)) {
1759    double * objectiveSimplex
1760      = copyOfArray(objective_,numberColumns_,0.0);
1761    double * rowObjectiveSimplex
1762      = copyOfArray(rowObjective_,numberRows_,0.0);
1763    int i;
1764    double largest;
1765    largest=0.0;
1766    for (i=0;i<numberRows_;i++) {
1767      rowObjectiveSimplex[i] *= optimizationDirection_;
1768      double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
1769      if (difference>largest)
1770        largest=difference;
1771    }
1772    for (i=0;i<numberColumns_;i++) {
1773      objectiveSimplex[i] *= optimizationDirection_;
1774      double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
1775      if (difference>largest)
1776        largest=difference;
1777    }
1778    if ((handler_->logLevel()&16))
1779      printf("difference in obj %g\n",largest);
1780    delete [] objectiveSimplex;
1781    delete [] rowObjectiveSimplex;
1782  }
1783#endif
1784  handler_->message(CLP_SIMPLEX_STATUS,messages_)
1785    <<numberIterations_<<objectiveValue();
1786  handler_->printing(sumPrimalInfeasibilities_>0.0)
1787    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
1788  handler_->printing(sumDualInfeasibilities_>0.0)
1789    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
1790  handler_->printing(numberDualInfeasibilitiesWithoutFree_
1791                     <numberDualInfeasibilities_)
1792                       <<numberDualInfeasibilities_-
1793    numberDualInfeasibilitiesWithoutFree_;
1794  handler_->message()<<CoinMessageEol;
1795  while (problemStatus_<=-3) {
1796    bool cleanDuals=false;
1797    int numberChangedBounds=0;
1798    int doOriginalTolerance=0;
1799    if ( lastCleaned==numberIterations_)
1800      doOriginalTolerance=1;
1801    // check optimal
1802    if (dualFeasible()||problemStatus_==-4) {
1803      if (primalFeasible()) {
1804        // may be optimal - or may be bounds are wrong
1805        handler_->message(CLP_DUAL_BOUNDS,messages_)
1806          <<dualBound_
1807          <<CoinMessageEol;
1808        // save solution in case unbounded
1809        CoinDisjointCopyN(columnActivityWork_,numberColumns_,
1810                          columnArray_[0]->denseVector());
1811        CoinDisjointCopyN(rowActivityWork_,numberRows_,
1812                          rowArray_[2]->denseVector());
1813        numberChangedBounds=changeBounds(false,rowArray_[0],changeCost);
1814        if (numberChangedBounds<=0) {
1815          //looks optimal - do we need to reset tolerance
1816          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
1817            doOriginalTolerance=2;
1818            numberTimesOptimal_++;
1819            changeMade_++; // say something changed
1820            if (numberTimesOptimal_==1) {
1821              dualTolerance_ = min(dualTolerance_,1.0e-8);
1822              // better to have small tolerance even if slower
1823              factorization_->zeroTolerance(1.0e-15);
1824            }
1825          } else {
1826            problemStatus_=0; // optimal
1827            if (lastCleaned<numberIterations_) {
1828              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1829                <<CoinMessageEol;
1830            }
1831          }
1832        } else {
1833          cleanDuals=true;
1834          if (doOriginalTolerance==1) {
1835            // check unbounded
1836            problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
1837                                            changeCost);
1838            if (problemStatus_==2) {
1839              // it is unbounded - restore solution
1840              // but first add in changes to non-basic
1841              int iColumn;
1842              double * original = columnArray_[0]->denseVector();
1843              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1844                if(getColumnStatus(iColumn)!= ClpSimplex::basic)
1845                  ray_[iColumn] += 
1846                    columnActivityWork_[iColumn]-original[iColumn];
1847                columnActivityWork_[iColumn] = original[iColumn];
1848              }
1849              CoinDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
1850                                rowActivityWork_);
1851            }
1852          } else {
1853            doOriginalTolerance=2;
1854            rowArray_[0]->clear();
1855          }
1856        }
1857        CoinFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
1858        CoinFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
1859      }
1860      if (problemStatus_==-4) {
1861        // may be infeasible - or may be bounds are wrong
1862        handler_->message(CLP_DUAL_CHECKB,messages_)
1863          <<dualBound_
1864          <<CoinMessageEol;
1865        numberChangedBounds=changeBounds(false,NULL,changeCost);
1866        if (numberChangedBounds<=0||dualBound_>1.0e20||
1867            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
1868          problemStatus_=1; // infeasible
1869        } else {
1870          problemStatus_=-1; //iterate
1871          cleanDuals=true;
1872          doOriginalTolerance=2;
1873          // and delete ray which has been created
1874          delete [] ray_;
1875          ray_ = NULL;
1876        }
1877      }
1878    } else {
1879      cleanDuals=true;
1880    }
1881    if (problemStatus_<0) {
1882      if (doOriginalTolerance==2) {
1883        // put back original tolerance
1884        lastCleaned=numberIterations_;
1885        handler_->message(CLP_DUAL_ORIGINAL,messages_)
1886          <<CoinMessageEol;
1887
1888        perturbation_=102; // stop any perturbations
1889        createRim(4);
1890        // make sure duals are current
1891        computeDuals();
1892        // put back bounds as they were if was optimal
1893        if (doOriginalTolerance==2) {
1894          changeMade_++; // say something changed
1895          changeBounds(true,NULL,changeCost);
1896          cleanDuals=true;
1897        }
1898      }
1899      if (cleanDuals) {
1900        // make sure dual feasible
1901        // look at all rows and columns
1902        rowArray_[0]->clear();
1903        CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
1904        rowArray_[0]->setNumElements(numberRows_);
1905        columnArray_[0]->clear();
1906        CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
1907        columnArray_[0]->setNumElements(numberColumns_);
1908        double objectiveChange=0.0;
1909        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
1910                          0.0,objectiveChange);
1911        // for now - recompute all
1912        gutsOfSolution(rowActivityWork_,columnActivityWork_);
1913        assert(numberDualInfeasibilitiesWithoutFree_==0);
1914        if (numberDualInfeasibilities_) {
1915          // bad free variables
1916          if (primalFeasible()) {
1917            std::cerr<<"Free variable problem?"<<std::endl;
1918            abort(); // what now
1919          }
1920          problemStatus_=-1; // carry on as normal
1921        }
1922      } else {
1923        // iterate
1924        problemStatus_=-1;
1925      }
1926    }
1927  }
1928  if (type==0||type==1) {
1929    if (!type) {
1930      // create save arrays
1931      delete [] saveStatus_;
1932      delete [] savedSolution_;
1933      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1934      savedSolution_ = new double [numberRows_+numberColumns_];
1935    }
1936    // save arrays
1937    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
1938    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
1939           numberRows_*sizeof(double));
1940    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
1941  }
1942
1943  // restore weights (if saved) - also recompute infeasibility list
1944  if (tentativeStatus>-3) 
1945    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
1946  else
1947    dualRowPivot_->saveWeights(this,3);
1948  // unflag all variables (we may want to wait a bit?)
1949  int iRow;
1950  for (iRow=0;iRow<numberRows_;iRow++) {
1951    int iPivot=pivotVariable_[iRow];
1952    clearFlagged(iPivot);
1953  }
1954  if (problemStatus_<0&&!changeMade_) {
1955    problemStatus_=4; // unknown
1956  }
1957
1958}
1959/* While updateDualsInDual sees what effect is of flip
1960   this does actual flipping.
1961   If change >0.0 then value in array >0.0 => from lower to upper
1962*/
1963void 
1964ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
1965                  CoinIndexedVector * columnArray,
1966                  double change)
1967{
1968  double * work;
1969  int number;
1970  int * which;
1971 
1972  int iSection;
1973
1974  for (iSection=0;iSection<2;iSection++) {
1975    int i;
1976    double * solution = solutionRegion(iSection);
1977    double * lower = lowerRegion(iSection);
1978    double * upper = upperRegion(iSection);
1979    int addSequence;
1980    if (!iSection) {
1981      work = rowArray->denseVector();
1982      number = rowArray->getNumElements();
1983      which = rowArray->getIndices();
1984      addSequence = numberColumns_;
1985    } else {
1986      work = columnArray->denseVector();
1987      number = columnArray->getNumElements();
1988      which = columnArray->getIndices();
1989      addSequence = 0;
1990    }
1991   
1992    for (i=0;i<number;i++) {
1993      int iSequence = which[i];
1994#ifndef NDEBUG
1995      double value = work[iSequence]*change;
1996#endif
1997      work[iSequence]=0.0;
1998      Status status = getStatus(iSequence+addSequence);
1999
2000      switch(status) {
2001
2002      case ClpSimplex::basic:
2003      case ClpSimplex::isFree:
2004      case ClpSimplex::superBasic:
2005        break;
2006      case ClpSimplex::atUpperBound:
2007        assert (value<=0.0);
2008        // to lower bound
2009        setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
2010        solution[iSequence] = lower[iSequence];
2011        break;
2012      case ClpSimplex::atLowerBound:
2013        assert (value>=0.0);
2014        // to upper bound
2015        setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
2016        solution[iSequence] = upper[iSequence];
2017        break;
2018      }
2019    }
2020  }
2021  rowArray->setNumElements(0);
2022  columnArray->setNumElements(0);
2023}
2024// Restores bound to original bound
2025void 
2026ClpSimplexDual::originalBound( int iSequence)
2027{
2028  if (iSequence>=numberColumns_) {
2029    // rows
2030    int iRow = iSequence-numberColumns_;
2031    rowLowerWork_[iRow]=rowLower_[iRow];
2032    rowUpperWork_[iRow]=rowUpper_[iRow];
2033    if (rowScale_) {
2034      if (rowLowerWork_[iRow]>-1.0e50)
2035        rowLowerWork_[iRow] *= rowScale_[iRow];
2036      if (rowUpperWork_[iRow]<1.0e50)
2037        rowUpperWork_[iRow] *= rowScale_[iRow];
2038    }
2039  } else {
2040    // columns
2041    columnLowerWork_[iSequence]=columnLower_[iSequence];
2042    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2043    if (rowScale_) {
2044      double multiplier = 1.0/columnScale_[iSequence];
2045      if (columnLowerWork_[iSequence]>-1.0e50)
2046        columnLowerWork_[iSequence] *= multiplier;
2047      if (columnUpperWork_[iSequence]<1.0e50)
2048        columnUpperWork_[iSequence] *= multiplier;
2049    }
2050  }
2051  setFakeBound(iSequence,ClpSimplex::noFake);
2052}
2053/* As changeBounds but just changes new bounds for a single variable.
2054   Returns true if change */
2055bool 
2056ClpSimplexDual::changeBound( int iSequence)
2057{
2058  // old values
2059  double oldLower=lower_[iSequence];
2060  double oldUpper=upper_[iSequence];
2061  double value=solution_[iSequence];
2062  bool modified=false;
2063  originalBound(iSequence);
2064  // original values
2065  double lowerValue=lower_[iSequence];
2066  double upperValue=upper_[iSequence];
2067  // back to altered values
2068  lower_[iSequence] = oldLower;
2069  upper_[iSequence] = oldUpper;
2070  if (value==oldLower) {
2071    if (upperValue > oldLower + dualBound_) {
2072      upper_[iSequence]=oldLower+dualBound_;
2073      setFakeBound(iSequence,ClpSimplex::upperFake);
2074      modified=true;
2075    }
2076  } else if (value==oldUpper) {
2077    if (lowerValue < oldUpper - dualBound_) {
2078      lower_[iSequence]=oldUpper-dualBound_;
2079      setFakeBound(iSequence,ClpSimplex::lowerFake);
2080      modified=true;
2081    }
2082  } else {
2083    assert(value==oldLower||value==oldUpper);
2084  }
2085  return modified;
2086}
2087// Perturbs problem
2088void 
2089ClpSimplexDual::perturb()
2090{
2091  if (perturbation_>100)
2092    return; //perturbed already
2093  int iRow,iColumn;
2094  // dual perturbation
2095  double perturbation=1.0e-20;
2096  // maximum fraction of cost to perturb
2097  double maximumFraction = 1.0e-4;
2098  if (perturbation_==100) {
2099    perturbation = 1.0e-4;
2100    for (iRow=0;iRow<numberRows_;iRow++) {
2101      double value = fabs(rowActivityWork_[iRow]*rowObjectiveWork_[iRow]);
2102      perturbation = max(perturbation,value);
2103    }
2104    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
2105      double value = 
2106        fabs(columnActivityWork_[iColumn]*objectiveWork_[iColumn]);
2107      perturbation = max(perturbation,value);
2108    }
2109    perturbation *= 1.0e-8;
2110  } else if (perturbation_<100) {
2111    perturbation = pow(10.0,perturbation_);
2112    // user is in charge
2113    maximumFraction = 1.0e100;
2114  }
2115  // modify costs
2116  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2117    <<perturbation
2118    <<CoinMessageEol;
2119  for (iRow=0;iRow<numberRows_;iRow++) {
2120    double value = perturbation;
2121    double currentValue = rowObjectiveWork_[iRow];
2122    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2123    if (rowLowerWork_[iRow]>-largeValue_) {
2124      if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
2125        value *= drand48();
2126      else
2127        value *= -drand48();
2128    } else if (rowUpperWork_[iRow]<largeValue_) {
2129      value *= -drand48();
2130    } else {
2131      value=0.0;
2132    }
2133    rowObjectiveWork_[iRow] += value;
2134  }
2135  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2136    double value = perturbation;
2137    double currentValue = objectiveWork_[iColumn];
2138    value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
2139    if (columnLowerWork_[iColumn]>-largeValue_) {
2140      if (fabs(columnLowerWork_[iColumn])<
2141          fabs(columnUpperWork_[iColumn])) 
2142        value *= drand48();
2143      else
2144        value *= -drand48();
2145    } else if (columnUpperWork_[iColumn]<largeValue_) {
2146      value *= -drand48();
2147    } else {
2148      value=0.0;
2149    }
2150    objectiveWork_[iColumn] += value;
2151  }
2152  // say perturbed
2153  perturbation_=102;
2154
2155}
2156/* For strong branching.  On input lower and upper are new bounds
2157   while on output they are change in objective function values
2158   (>1.0e50 infeasible).
2159   Return code is 0 if nothing interesting, -1 if infeasible both
2160   ways and +1 if infeasible one way (check values to see which one(s))
2161*/
2162int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
2163                                    double * newLower, double * newUpper,
2164                                    bool stopOnFirstInfeasible,
2165                                    bool alwaysFinish)
2166{
2167  int i;
2168  int returnCode=0;
2169  double saveObjective = objectiveValue_;
2170#if 1
2171  algorithm_ = -1;
2172
2173  //scaling(false);
2174
2175  // put in standard form (and make row copy)
2176  // create modifiable copies of model rim and do optional scaling
2177  createRim(7+8+16,true);
2178
2179  // change newLower and newUpper if scaled
2180
2181  // Do initial factorization
2182  // and set certain stuff
2183  // We can either set increasing rows so ...IsBasic gives pivot row
2184  // or we can just increment iBasic one by one
2185  // for now let ...iBasic give pivot row
2186  factorization_->increasingRows(2);
2187  // row activities have negative sign
2188  factorization_->slackValue(-1.0);
2189  factorization_->zeroTolerance(1.0e-13);
2190  // save if sparse factorization wanted
2191  int saveSparse = factorization_->sparseThreshold();
2192
2193  int factorizationStatus = internalFactorize(0);
2194  if (factorizationStatus<0)
2195    return 1; // some error
2196  else if (factorizationStatus)
2197    handler_->message(CLP_SINGULARITIES,messages_)
2198      <<factorizationStatus
2199      <<CoinMessageEol;
2200  if (saveSparse) {
2201    // use default at present
2202    factorization_->sparseThreshold(0);
2203    factorization_->goSparse();
2204  }
2205 
2206  // save stuff
2207  ClpFactorization saveFactorization(*factorization_);
2208  // save basis and solution
2209  double * saveSolution = new double[numberRows_+numberColumns_];
2210  memcpy(saveSolution,solution_,
2211         (numberRows_+numberColumns_)*sizeof(double));
2212  unsigned char * saveStatus = 
2213    new unsigned char [numberRows_+numberColumns_];
2214  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2215  // save bounds as createRim makes clean copies
2216  double * saveLower = new double[numberRows_+numberColumns_];
2217  memcpy(saveLower,lower_,
2218         (numberRows_+numberColumns_)*sizeof(double));
2219  double * saveUpper = new double[numberRows_+numberColumns_];
2220  memcpy(saveUpper,upper_,
2221         (numberRows_+numberColumns_)*sizeof(double));
2222  int * savePivot = new int [numberRows_];
2223  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
2224  // need to save/restore weights.
2225
2226  for (i=0;i<numberVariables;i++) {
2227    int iColumn = variables[i];
2228    double objectiveChange;
2229    double saveBound;
2230   
2231    // try down
2232   
2233    saveBound = columnUpper_[iColumn];
2234    // external view - in case really getting optimal
2235    columnUpper_[iColumn] = newUpper[i];
2236    if (scalingFlag_<=0) 
2237      upper_[iColumn] = newUpper[i];
2238    else 
2239      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
2240    // Start of fast iterations
2241    int status = fastDual(alwaysFinish);
2242
2243    // restore
2244    memcpy(solution_,saveSolution,
2245           (numberRows_+numberColumns_)*sizeof(double));
2246    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2247    memcpy(lower_,saveLower,
2248           (numberRows_+numberColumns_)*sizeof(double));
2249    memcpy(upper_,saveUpper,
2250           (numberRows_+numberColumns_)*sizeof(double));
2251    columnUpper_[iColumn] = saveBound;
2252    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2253    delete factorization_;
2254    factorization_ = new ClpFactorization(saveFactorization);
2255
2256    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2257      objectiveChange = objectiveValue_-saveObjective;
2258    } else {
2259      objectiveChange = 1.0e100;
2260    }
2261    newUpper[i]=objectiveChange;
2262#ifdef CLP_DEBUG
2263    printf("down on %d costs %g\n",iColumn,objectiveChange);
2264#endif
2265         
2266    // try up
2267   
2268    saveBound = columnLower_[iColumn];
2269    // external view - in case really getting optimal
2270    columnLower_[iColumn] = newLower[i];
2271    if (scalingFlag_<=0) 
2272      lower_[iColumn] = newLower[i];
2273    else 
2274      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
2275    // Start of fast iterations
2276    status = fastDual(alwaysFinish);
2277
2278    // restore
2279    memcpy(solution_,saveSolution,
2280           (numberRows_+numberColumns_)*sizeof(double));
2281    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2282    memcpy(lower_,saveLower,
2283           (numberRows_+numberColumns_)*sizeof(double));
2284    memcpy(upper_,saveUpper,
2285           (numberRows_+numberColumns_)*sizeof(double));
2286    columnLower_[iColumn] = saveBound;
2287    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
2288    delete factorization_;
2289    factorization_ = new ClpFactorization(saveFactorization);
2290
2291    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
2292      objectiveChange = objectiveValue_-saveObjective;
2293    } else {
2294      objectiveChange = 1.0e100;
2295    }
2296    newLower[i]=objectiveChange;
2297#ifdef CLP_DEBUG
2298    printf("up on %d costs %g\n",iColumn,objectiveChange);
2299#endif
2300         
2301    /* Possibilities are:
2302       Both sides feasible - store
2303       Neither side feasible - set objective high and exit
2304       One side feasible - change bounds and resolve
2305    */
2306    if (newUpper[i]<1.0e100) {
2307      if(newLower[i]<1.0e100) {
2308        // feasible - no action
2309      } else {
2310        // up feasible, down infeasible
2311        returnCode=1;
2312        if (stopOnFirstInfeasible)
2313          break;
2314      }
2315    } else {
2316      if(newLower[i]<1.0e100) {
2317        // down feasible, up infeasible
2318        returnCode=1;
2319        if (stopOnFirstInfeasible)
2320          break;
2321      } else {
2322        // neither side feasible
2323        returnCode=-1;
2324        break;
2325      }
2326    }
2327  }
2328  delete [] saveSolution;
2329  delete [] saveLower;
2330  delete [] saveUpper;
2331  delete [] saveStatus;
2332  delete [] savePivot;
2333
2334  // at present we are leaving factorization around
2335  // maybe we should empty it
2336  deleteRim();
2337  factorization_->sparseThreshold(saveSparse);
2338#else
2339  // save basis and solution
2340  double * rowSolution = new double[numberRows_];
2341  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
2342  double * columnSolution = new double[numberColumns_];
2343  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
2344  unsigned char * saveStatus = 
2345    new unsigned char [numberRows_+numberColumns_];
2346  if (!status_) {
2347    // odd but anyway setup all slack basis
2348    status_ = new unsigned char [numberColumns_+numberRows_];
2349    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2350  }
2351  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
2352  for (i=0;i<numberVariables;i++) {
2353    int iColumn = variables[i];
2354    double objectiveChange;
2355   
2356    // try down
2357   
2358    double saveUpper = columnUpper_[iColumn];
2359    columnUpper_[iColumn] = newUpper[i];
2360    dual();
2361
2362    // restore
2363    columnUpper_[iColumn] = saveUpper;
2364    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2365    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2366    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2367
2368    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2369      objectiveChange = objectiveValue_-saveObjective;
2370    } else {
2371      objectiveChange = 1.0e100;
2372    }
2373    newUpper[i]=objectiveChange;
2374#ifdef CLP_DEBUG
2375    printf("down on %d costs %g\n",iColumn,objectiveChange);
2376#endif
2377         
2378    // try up
2379   
2380    double saveLower = columnLower_[iColumn];
2381    columnLower_[iColumn] = newLower[i];
2382    dual();
2383
2384    // restore
2385    columnLower_[iColumn] = saveLower;
2386    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
2387    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
2388    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
2389
2390    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
2391      objectiveChange = objectiveValue_-saveObjective;
2392    } else {
2393      objectiveChange = 1.0e100;
2394    }
2395    newLower[i]=objectiveChange;
2396#ifdef CLP_DEBUG
2397    printf("up on %d costs %g\n",iColumn,objectiveChange);
2398#endif
2399         
2400    /* Possibilities are:
2401       Both sides feasible - store
2402       Neither side feasible - set objective high and exit
2403       One side feasible - change bounds and resolve
2404    */
2405    if (newUpper[i]<1.0e100) {
2406      if(newLower[i]<1.0e100) {
2407        // feasible - no action
2408      } else {
2409        // up feasible, down infeasible
2410        returnCode=1;
2411        if (stopOnFirstInfeasible)
2412          break;
2413      }
2414    } else {
2415      if(newLower[i]<1.0e100) {
2416        // down feasible, up infeasible
2417        returnCode=1;
2418        if (stopOnFirstInfeasible)
2419          break;
2420      } else {
2421        // neither side feasible
2422        returnCode=-1;
2423        break;
2424      }
2425    }
2426  }
2427  delete [] rowSolution;
2428  delete [] columnSolution;
2429  delete [] saveStatus;
2430#endif
2431  objectiveValue_ = saveObjective;
2432  return returnCode;
2433}
2434// treat no pivot as finished (unless interesting)
2435int ClpSimplexDual::fastDual(bool alwaysFinish)
2436{
2437  algorithm_ = -1;
2438  dualTolerance_=dblParam_[ClpDualTolerance];
2439  primalTolerance_=dblParam_[ClpPrimalTolerance];
2440
2441  // save dual bound
2442  double saveDualBound = dualBound_;
2443
2444  double objectiveChange;
2445  // for dual we will change bounds using dualBound_
2446  // for this we need clean basis so it is after factorize
2447  gutsOfSolution(rowActivityWork_,columnActivityWork_);
2448  numberFake_ =0; // Number of variables at fake bounds
2449  changeBounds(true,NULL,objectiveChange);
2450
2451  problemStatus_ = -1;
2452  numberIterations_=0;
2453
2454  int lastCleaned=0; // last time objective or bounds cleaned up
2455
2456  // number of times we have declared optimality
2457  numberTimesOptimal_=0;
2458
2459  // This says whether to restore things etc
2460  int factorType=0;
2461  /*
2462    Status of problem:
2463    0 - optimal
2464    1 - infeasible
2465    2 - unbounded
2466    -1 - iterating
2467    -2 - factorization wanted
2468    -3 - redo checking without factorization
2469    -4 - looks infeasible
2470
2471    BUT also from whileIterating return code is:
2472
2473   -1 iterations etc
2474   -2 inaccuracy
2475   -3 slight inaccuracy (and done iterations)
2476   +0 looks optimal (might be unbounded - but we will investigate)
2477   +1 looks infeasible
2478   +3 max iterations
2479
2480  */
2481
2482  int returnCode = 0;
2483
2484  while (problemStatus_<0) {
2485    int iRow,iColumn;
2486    // clear
2487    for (iRow=0;iRow<4;iRow++) {
2488      rowArray_[iRow]->clear();
2489    }   
2490   
2491    for (iColumn=0;iColumn<2;iColumn++) {
2492      columnArray_[iColumn]->clear();
2493    }   
2494
2495    // give matrix (and model costs and bounds a chance to be
2496    // refreshed (normally null)
2497    matrix_->refresh(this);
2498    // may factorize, checks if problem finished
2499    // should be able to speed this up on first time
2500    statusOfProblemInDual(lastCleaned,factorType);
2501
2502    // Say good factorization
2503    factorType=1;
2504
2505    // Do iterations
2506    if (problemStatus_<0) {
2507#if 1
2508      returnCode = whileIterating();
2509      if (!alwaysFinish&&returnCode<1) {
2510        double limit = 0.0;
2511        getDblParam(ClpDualObjectiveLimit, limit);
2512        if(objectiveValue()*optimizationDirection_<limit|| 
2513           numberAtFakeBound()) {
2514          // can't say anything interesting - might as well return
2515#ifdef CLP_DEBUG
2516          printf("returning from fastDual after %d iterations with code %d\n",
2517                 numberIterations_,returnCode);
2518#endif
2519          returnCode=1;
2520          break;
2521        }
2522      }
2523      returnCode=0;
2524#else
2525      whileIterating();
2526#endif
2527    }
2528  }
2529
2530  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
2531  dualBound_ = saveDualBound;
2532  return returnCode;
2533}
2534/* Checks number of variables at fake bounds.  This is used by fastDual
2535   so can exit gracefully before end */
2536int 
2537ClpSimplexDual::numberAtFakeBound()
2538{
2539  int iSequence;
2540  int numberFake=0;
2541 
2542  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2543    FakeBound bound = getFakeBound(iSequence);
2544    switch(getStatus(iSequence)) {
2545
2546    case ClpSimplex::basic:
2547      break;
2548    case ClpSimplex::isFree:
2549    case ClpSimplex::superBasic:
2550      break;
2551    case ClpSimplex::atUpperBound:
2552      if (bound==ClpSimplex::upperFake||bound==ClpSimplex::bothFake)
2553        numberFake++;
2554      break;
2555    case ClpSimplex::atLowerBound:
2556      if (bound==ClpSimplex::lowerFake||bound==ClpSimplex::bothFake)
2557        numberFake++;
2558      break;
2559    }
2560  }
2561  return numberFake;
2562}
Note: See TracBrowser for help on using the repository browser.