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

Last change on this file since 21 was 21, checked in by forrest, 18 years ago

Check matrix more carefully

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