# source:branches/devel-1/ClpSimplexDual.cpp@29

Last change on this file since 29 was 29, checked in by forrest, 19 years ago

Presolve (no changes to Makefile)

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