source: branches/pre/ClpPdco.cpp @ 212

Last change on this file since 212 was 212, checked in by forrest, 17 years ago

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 118.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5/* Notes on implementation of dual simplex algorithm.
6
7   When dual feasible:
8
9   If primal feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to leave basis (normally going to nearest bound) (B).  We
11   now need to find an incoming variable which will leave problem
12   dual feasible so we get the row of the tableau corresponding to
13   the basic variable (with the correct sign depending if basic variable
14   above or below feasibility region - as that affects whether reduced
15   cost on outgoing variable has to be positive or negative).
16
17   We now perform a ratio test to determine which incoming variable will
18   preserve dual feasibility (C).  If no variable found then problem
19   is infeasible (in primal sense).  If there is a variable, we then
20   perform pivot and repeat.  Trivial?
21
22   -------------------------------------------
23
24   A) How do we get dual feasible?  If all variables have bounds then
25   it is trivial to get feasible by putting non-basic variables to
26   correct bounds.  OSL did not have a phase 1/phase 2 approach but
27   instead effectively put fake bounds on variables and this is the
28   approach here, although I had hoped to make it cleaner.
29
30   If there is a weight of X on getting dual feasible:
31     Non-basic variables with negative reduced costs are put to
32     lesser of their upper bound and their lower bound + X.
33     Similarly, mutatis mutandis, for positive reduced costs.
34
35   Free variables should normally be in basis, otherwise I have
36   coding which may be able to come out (and may not be correct).
37
38   In OSL, this weight was changed heuristically, here at present
39   it is only increased if problem looks finished.  If problem is
40   feasible I check for unboundedness.  If not unbounded we
41   could play with going into primal.  As long as weights increase
42   any algorithm would be finite.
43   
44   B) Which outgoing variable to choose is a virtual base class.
45   For difficult problems steepest edge is preferred while for
46   very easy (large) problems we will need partial scan.
47
48   C) Sounds easy, but this is hardest part of algorithm.
49      1) Instead of stopping at first choice, we may be able
50      to flip that variable to other bound and if objective
51      still improving choose again.  These mini iterations can
52      increase speed by orders of magnitude but we may need to
53      go to more of a bucket choice of variable rather than looking
54      at them one by one (for speed).
55      2) Accuracy.  Reduced costs may be of wrong sign but less than
56      tolerance.  Pivoting on these makes objective go backwards.
57      OSL modified cost so a zero move was made, Gill et al
58      (in primal analogue) modified so a strictly positive move was
59      made.  It is not quite as neat in dual but that is what we
60      try and do.  The two problems are that re-factorizations can
61      change reduced costs above and below tolerances and that when
62      finished we need to reset costs and try again.
63      3) Degeneracy.  Gill et al helps but may not be enough.  We
64      may need more.  Also it can improve speed a lot if we perturb
65      the costs significantly. 
66
67  References:
68     Forrest and Goldfarb, Steepest-edge simplex algorithms for
69       linear programming - Mathematical Programming 1992
70     Forrest and Tomlin, Implementing the simplex method for
71       the Optimization Subroutine Library - IBM Systems Journal 1992
72     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
73       Procedure for Linear and Nonlinear Programming SOL report 1988
74
75
76  TODO:
77 
78  a) Better recovery procedures.  At present I never check on forward
79     progress.  There is checkpoint/restart with reducing
80     re-factorization frequency, but this is only on singular
81     factorizations.
82  b) Fast methods for large easy problems (and also the option for
83     the code to automatically choose which method).
84  c) We need to be able to stop in various ways for OSI - this
85     is fairly easy.
86
87 */
88
89
90#include "CoinPragma.hpp"
91
92#include <math.h>
93
94#include "CoinHelperFunctions.hpp"
95#include "ClpPdco.hpp"
96#include "ClpFactorization.hpp"
97#include "CoinPackedMatrix.hpp"
98#include "CoinIndexedVector.hpp"
99#include "CoinWarmStartBasis.hpp"
100#include "ClpDualRowDantzig.hpp"
101#include "ClpPlusMinusOneMatrix.hpp"
102#include "ClpMessage.hpp"
103#include <cfloat>
104#include <cassert>
105#include <string>
106#include <stdio.h>
107#include <iostream>
108//#define CLP_DEBUG 1
109#define CHECK_DJ
110// dual
111int ClpPdco::dual (int ifValuesPass )
112{
113
114  /* *** Method
115     This is a vanilla version of dual simplex.
116
117     It tries to be a single phase approach with a weight of 1.0 being
118     given to getting optimal and a weight of dualBound_ being
119     given to getting dual feasible.  In this version I have used the
120     idea that this weight can be thought of as a fake bound.  If the
121     distance between the lower and upper bounds on a variable is less
122     than the feasibility weight then we are always better off flipping
123     to other bound to make dual feasible.  If the distance is greater
124     then we make up a fake bound dualBound_ away from one bound.
125     If we end up optimal or primal infeasible, we check to see if
126     bounds okay.  If so we have finished, if not we increase dualBound_
127     and continue (after checking if unbounded). I am undecided about
128     free variables - there is coding but I am not sure about it.  At
129     present I put them in basis anyway.
130
131     The code is designed to take advantage of sparsity so arrays are
132     seldom zeroed out from scratch or gone over in their entirety.
133     The only exception is a full scan to find outgoing variable.  This
134     will be changed to keep an updated list of infeasibilities (or squares
135     if steepest edge).  Also on easy problems we don't need full scan - just
136     pick first reasonable.
137
138     One problem is how to tackle degeneracy and accuracy.  At present
139     I am using the modification of costs which I put in OSL and which was
140     extended by Gill et al.  I am still not sure of the exact details.
141
142     The flow of dual is three while loops as follows:
143
144     while (not finished) {
145
146       while (not clean solution) {
147
148          Factorize and/or clean up solution by flipping variables so
149          dual feasible.  If looks finished check fake dual bounds.
150          Repeat until status is iterating (-1) or finished (0,1,2)
151
152       }
153
154       while (status==-1) {
155
156         Iterate until no pivot in or out or time to re-factorize.
157
158         Flow is:
159
160         choose pivot row (outgoing variable).  if none then
161         we are primal feasible so looks as if done but we need to
162         break and check bounds etc.
163
164         Get pivot row in tableau
165
166         Choose incoming column.  If we don't find one then we look
167         primal infeasible so break and check bounds etc.  (Also the
168         pivot tolerance is larger after any iterations so that may be
169         reason)
170
171         If we do find incoming column, we may have to adjust costs to
172         keep going forwards (anti-degeneracy).  Check pivot will be stable
173         and if unstable throw away iteration (we will need to implement
174         flagging of basic variables sometime) and break to re-factorize.
175         If minor error re-factorize after iteration.
176
177         Update everything (this may involve flipping variables to stay
178         dual feasible.
179
180       }
181
182     }
183
184     At present we never check we are going forwards.  I overdid that in
185     OSL so will try and make a last resort.
186
187     Needs partial scan pivot out option.
188     Needs dantzig, uninitialized and full steepest edge options (can still
189     use partial scan)
190
191     May need other anti-degeneracy measures, especially if we try and use
192     loose tolerances as a way to solve in fewer iterations.
193
194     I like idea of dynamic scaling.  This gives opportunity to decouple
195     different implications of scaling for accuracy, iteration count and
196     feasibility tolerance.
197
198  */
199  algorithm_ = -1;
200
201  // save data
202  ClpDataSave data = saveData();
203
204  // Save so can see if doing after primal
205  int initialStatus=problemStatus_;
206
207  // If values pass then save given duals round check solution
208  double * saveDuals = NULL;
209  if (ifValuesPass) {
210    saveDuals = new double [numberRows_+numberColumns_];
211    memcpy(saveDuals,dual_,numberRows_*sizeof(double));
212  }
213 
214  // sanity check
215  // initialize - no values pass and algorithm_ is -1
216  // put in standard form (and make row copy)
217  // create modifiable copies of model rim and do optional scaling
218  // If problem looks okay
219  // Do initial factorization
220  // If user asked for perturbation - do it
221  if (!startup(0)) {
222    // looks okay
223   
224    // If values pass then scale pi
225    if (ifValuesPass) {
226      if (problemStatus_&&perturbation_<100) 
227        perturb();
228      int i;
229      if (scalingFlag_>0) {
230        for (i=0;i<numberRows_;i++) {
231          dual_[i] = saveDuals[i]/rowScale_[i];
232        }
233      } else {
234        memcpy(dual_,saveDuals,numberRows_*sizeof(double));
235      }
236      // now create my duals
237      for (i=0;i<numberRows_;i++) {
238        // slack
239        double value = dual_[i];
240        value += rowObjectiveWork_[i];
241        saveDuals[i+numberColumns_]=value;
242      }
243      ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
244      transposeTimes(-1.0,dual_,saveDuals);
245      memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
246      // set up possible ones
247      for (i=0;i<numberRows_+numberColumns_;i++)
248        clearPivoted(i);
249      int iRow;
250      for (iRow=0;iRow<numberRows_;iRow++) {
251        int iPivot=pivotVariable_[iRow];
252        if (fabs(saveDuals[iPivot])>dualTolerance_) {
253          if (getStatus(iPivot)!=isFree) 
254            setPivoted(iPivot);
255        }
256      }
257    } 
258
259    double objectiveChange;
260    numberFake_ =0; // Number of variables at fake bounds
261    changeBounds(true,NULL,objectiveChange);
262   
263    int lastCleaned=0; // last time objective or bounds cleaned up
264
265    if (!ifValuesPass) {
266      // Check optimal
267      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
268        problemStatus_=0;
269    }
270    if (problemStatus_<0&&perturbation_<100) {
271      perturb();
272      // Can't get here if values pass
273      gutsOfSolution(NULL,NULL);
274      if (handler_->logLevel()>2) {
275        handler_->message(CLP_SIMPLEX_STATUS,messages_)
276          <<numberIterations_<<objectiveValue();
277        handler_->printing(sumPrimalInfeasibilities_>0.0)
278          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
279        handler_->printing(sumDualInfeasibilities_>0.0)
280          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
281        handler_->printing(numberDualInfeasibilitiesWithoutFree_
282                           <numberDualInfeasibilities_)
283                             <<numberDualInfeasibilitiesWithoutFree_;
284        handler_->message()<<CoinMessageEol;
285      }
286    }
287       
288    // This says whether to restore things etc
289    int factorType=0;
290    /*
291      Status of problem:
292      0 - optimal
293      1 - infeasible
294      2 - unbounded
295      -1 - iterating
296      -2 - factorization wanted
297      -3 - redo checking without factorization
298      -4 - looks infeasible
299    */
300    while (problemStatus_<0) {
301      int iRow, iColumn;
302      // clear
303      for (iRow=0;iRow<4;iRow++) {
304        rowArray_[iRow]->clear();
305      }   
306     
307      for (iColumn=0;iColumn<2;iColumn++) {
308        columnArray_[iColumn]->clear();
309      }   
310     
311      // give matrix (and model costs and bounds a chance to be
312      // refreshed (normally null)
313      matrix_->refresh(this);
314      // If getting nowhere - why not give it a kick
315      // does not seem to work too well - do some more work
316      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
317          &&initialStatus!=10) {
318        perturb();
319        // Can't get here if values pass
320        gutsOfSolution(NULL,NULL);
321        if (handler_->logLevel()>2) {
322          handler_->message(CLP_SIMPLEX_STATUS,messages_)
323            <<numberIterations_<<objectiveValue();
324          handler_->printing(sumPrimalInfeasibilities_>0.0)
325            <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
326          handler_->printing(sumDualInfeasibilities_>0.0)
327            <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
328          handler_->printing(numberDualInfeasibilitiesWithoutFree_
329                             <numberDualInfeasibilities_)
330                               <<numberDualInfeasibilitiesWithoutFree_;
331          handler_->message()<<CoinMessageEol;
332        }
333      }
334      // may factorize, checks if problem finished
335      statusOfProblemInDual(lastCleaned,factorType,progress_,saveDuals);
336      // If values pass then do easy ones on first time
337      if (ifValuesPass&&
338          progress_->lastIterationNumber(0)<0) {
339        doEasyOnesInValuesPass(saveDuals);
340      }
341     
342      // Say good factorization
343      factorType=1;
344      if (data.sparseThreshold_) {
345        // use default at present
346        factorization_->sparseThreshold(0);
347        factorization_->goSparse();
348      }
349
350      // exit if victory declared
351      if (problemStatus_>=0)
352        break;
353     
354      // Do iterations
355      whileIterating(saveDuals);
356    }
357  }
358
359  assert (problemStatus_||!sumPrimalInfeasibilities_);
360
361  // clean up
362  finish();
363  delete [] saveDuals;
364
365  // Restore any saved stuff
366  restoreData(data);
367  return problemStatus_;
368}
369/* Reasons to come out:
370   -1 iterations etc
371   -2 inaccuracy
372   -3 slight inaccuracy (and done iterations)
373   +0 looks optimal (might be unbounded - but we will investigate)
374   +1 looks infeasible
375   +3 max iterations
376 */
377int
378ClpPdco::whileIterating(double * & givenDuals)
379{
380#ifdef CLP_DEBUG
381  int debugIteration=-1;
382#endif
383  {
384    int i;
385    for (i=0;i<4;i++) {
386      rowArray_[i]->clear();
387    }   
388    for (i=0;i<2;i++) {
389      columnArray_[i]->clear();
390    }   
391  }     
392  // if can't trust much and long way from optimal then relax
393  if (largestPrimalError_>10.0)
394    factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
395  else
396    factorization_->relaxAccuracyCheck(1.0);
397  // status stays at -1 while iterating, >=0 finished, -2 to invert
398  // status -3 to go to top without an invert
399  int returnCode = -1;
400
401  // compute average infeasibility for backward test
402  double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
403    ((double ) (numberPrimalInfeasibilities_+1));
404
405  // Get dubious weights
406  factorization_->getWeights(rowArray_[0]->getIndices());
407  CoinBigIndex * dubiousWeights = matrix_->dubiousWeights(this,rowArray_[0]->getIndices());
408  // If values pass then get list of candidates
409  int * candidateList = NULL;
410  int numberCandidates = 0;
411#ifdef CLP_DEBUG
412  bool wasInValuesPass= (givenDuals!=NULL);
413#endif
414  int candidate=-1;
415  if (givenDuals) {
416    candidateList = new int[numberRows_];
417    // move reduced costs across
418    memcpy(dj_,givenDuals,(numberRows_+numberColumns_)*sizeof(double));
419    int iRow;
420    for (iRow=0;iRow<numberRows_;iRow++) {
421      int iPivot=pivotVariable_[iRow];
422      if (flagged(iPivot))
423        continue;
424      if (fabs(dj_[iPivot])>dualTolerance_) {
425        if (pivoted(iPivot)) 
426          candidateList[numberCandidates++]=iRow;
427      } else {
428        clearPivoted(iPivot);
429      }
430    }
431    // and set first candidate
432    if (!numberCandidates) {
433      delete [] candidateList;
434      delete [] givenDuals;
435      givenDuals=NULL;
436      candidateList=NULL;
437      int iRow;
438      for (iRow=0;iRow<numberRows_;iRow++) {
439        int iPivot=pivotVariable_[iRow];
440        clearPivoted(iPivot);
441      }
442    }
443  }
444  while (problemStatus_==-1) {
445#ifdef CLP_DEBUG
446    if (givenDuals) {
447      double value5=0.0;
448      int i;
449      for (i=0;i<numberRows_+numberColumns_;i++) {
450        if (dj_[i]<-1.0e-6)
451          if (upper_[i]<1.0e20)
452            value5 += dj_[i]*upper_[i];
453          else
454            printf("bad dj %g on %d with large upper status %d\n",
455                   dj_[i],i,status_[i]&7);
456        else if (dj_[i] >1.0e-6)
457          if (lower_[i]>-1.0e20)
458            value5 += dj_[i]*lower_[i];
459          else
460            printf("bad dj %g on %d with large lower status %d\n",
461                   dj_[i],i,status_[i]&7);
462      }
463      printf("Values objective Value %g\n",value5);
464    }
465    if ((handler_->logLevel()&32)&&wasInValuesPass) {
466      double value5=0.0;
467      int i;
468      for (i=0;i<numberRows_+numberColumns_;i++) {
469        if (dj_[i]<-1.0e-6)
470          if (upper_[i]<1.0e20)
471            value5 += dj_[i]*upper_[i];
472        else if (dj_[i] >1.0e-6)
473          if (lower_[i]>-1.0e20)
474            value5 += dj_[i]*lower_[i];
475      }
476      printf("Values objective Value %g\n",value5);
477      {
478        int i;
479        for (i=0;i<numberRows_+numberColumns_;i++) {
480          int iSequence = i;
481          double oldValue;
482         
483          switch(getStatus(iSequence)) {
484           
485          case basic:
486          case ClpSimplex::isFixed:
487            break;
488          case isFree:
489          case superBasic:
490            abort();
491            break;
492          case atUpperBound:
493            oldValue = dj_[iSequence];
494            //assert (oldValue<=tolerance);
495            assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
496            break;
497          case atLowerBound:
498            oldValue = dj_[iSequence];
499            //assert (oldValue>=-tolerance);
500            assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
501            break;
502          }
503        }
504      }
505    }
506#endif
507#ifdef CLP_DEBUG
508    {
509      int i;
510      for (i=0;i<4;i++) {
511        rowArray_[i]->checkClear();
512      }   
513      for (i=0;i<2;i++) {
514        columnArray_[i]->checkClear();
515      }   
516    }     
517#endif
518#if CLP_DEBUG>2
519    // very expensive
520    if (numberIterations_>0&&numberIterations_<-801) {
521      handler_->setLogLevel(63);
522      double saveValue = objectiveValue_;
523      double * saveRow1 = new double[numberRows_];
524      double * saveRow2 = new double[numberRows_];
525      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
526      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
527      double * saveColumn1 = new double[numberColumns_];
528      double * saveColumn2 = new double[numberColumns_];
529      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
530      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
531      gutsOfSolution(NULL,NULL);
532      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
533             numberIterations_,
534             saveValue,objectiveValue_,sumDualInfeasibilities_);
535      if (saveValue>objectiveValue_+1.0e-2)
536        printf("**bad**\n");
537      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
538      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
539      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
540      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
541      delete [] saveRow1;
542      delete [] saveRow2;
543      delete [] saveColumn1;
544      delete [] saveColumn2;
545      objectiveValue_=saveValue;
546    }
547#endif
548#if 0
549    if (!factorization_->pivots()){
550      int iPivot;
551      double * array = rowArray_[3]->denseVector();
552      int i;
553      for (iPivot=0;iPivot<numberRows_;iPivot++) {
554        int iSequence = pivotVariable_[iPivot];
555        unpack(rowArray_[3],iSequence);
556        factorization_->updateColumn(rowArray_[2],rowArray_[3]);
557        assert (fabs(array[iPivot]-1.0)<1.0e-4);
558        array[iPivot]=0.0;
559        for (i=0;i<numberRows_;i++)
560          assert (fabs(array[i])<1.0e-4);
561        rowArray_[3]->clear();
562      }
563    }
564#endif
565#ifdef CLP_DEBUG
566    {
567      int iSequence, number=numberRows_+numberColumns_;
568      for (iSequence=0;iSequence<number;iSequence++) {
569        double lowerValue=lower_[iSequence];
570        double upperValue=upper_[iSequence];
571        double value=solution_[iSequence];
572        if(getStatus(iSequence)!=basic) {
573          assert(lowerValue>-1.0e20);
574          assert(upperValue<1.0e20);
575        }
576        switch(getStatus(iSequence)) {
577         
578        case basic:
579          break;
580        case isFree:
581        case superBasic:
582          break;
583        case atUpperBound:
584          assert (fabs(value-upperValue)<=primalTolerance_) ;
585          break;
586        case atLowerBound:
587        case ClpSimplex::isFixed:
588          assert (fabs(value-lowerValue)<=primalTolerance_) ;
589          break;
590        }
591      }
592    }
593    if(numberIterations_==debugIteration) {
594      printf("dodgy iteration coming up\n");
595    }
596#endif
597    // choose row to go out
598    // dualRow will go to virtual row pivot choice algorithm
599    // make sure values pass off if it should be
600    if (numberCandidates) 
601      candidate = candidateList[--numberCandidates];
602    else
603      candidate=-1;
604    dualRow(candidate);
605    if (pivotRow_>=0) {
606      // we found a pivot row
607      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
608        <<pivotRow_
609        <<CoinMessageEol;
610      // check accuracy of weights
611      dualRowPivot_->checkAccuracy();
612      // Get good size for pivot
613      double acceptablePivot=1.0e-7;
614      if (factorization_->pivots()>10)
615        acceptablePivot=1.0e-5; // if we have iterated be more strict
616      else if (factorization_->pivots()>5)
617        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
618      // get sign for finding row of tableau
619      if (candidate<0) {
620        // normal iteration
621        // create as packed
622        double direction=directionOut_;
623        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
624        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
625        // put row of tableau in rowArray[0] and columnArray[0]
626        matrix_->transposeTimes(this,-1.0,
627                              rowArray_[0],rowArray_[3],columnArray_[0]);
628        // do ratio test for normal iteration
629        dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
630                 rowArray_[3],acceptablePivot,dubiousWeights);
631      } else {
632        double direction=directionOut_;
633        rowArray_[0]->createPacked(1,&pivotRow_,&direction);
634        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
635        // put row of tableau in rowArray[0] and columnArray[0]
636        matrix_->transposeTimes(this,-1.0,
637                              rowArray_[0],rowArray_[3],columnArray_[0]);
638        acceptablePivot *= 10.0;
639        // do ratio test
640        checkPossibleValuesMove(rowArray_[0],columnArray_[0],
641                                            acceptablePivot,NULL);
642
643        // recompute true dualOut_
644        if (directionOut_<0) {
645          dualOut_ = valueOut_ - upperOut_;
646        } else {
647          dualOut_ = lowerOut_ - valueOut_;
648        }
649        // check what happened if was values pass
650        // may want to move part way i.e. movement
651        bool normalIteration = (sequenceIn_!=sequenceOut_);
652
653        clearPivoted(sequenceOut_);  // make sure won't be done again
654        // see if end of values pass
655        if (!numberCandidates) {
656          int iRow;
657          delete [] candidateList;
658          delete [] givenDuals;
659          candidate=-2; // -2 signals end
660          givenDuals=NULL;
661          handler_->message(CLP_END_VALUES_PASS,messages_)
662            <<numberIterations_;
663          candidateList=NULL;
664          for (iRow=0;iRow<numberRows_;iRow++) {
665            int iPivot=pivotVariable_[iRow];
666            //assert (fabs(dj_[iPivot]),1.0e-5);
667            clearPivoted(iPivot);
668          }
669        }
670        if (!normalIteration) {
671          //rowArray_[0]->cleanAndPackSafe(1.0e-60);
672          //columnArray_[0]->cleanAndPackSafe(1.0e-60);
673          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
674          if (candidate==-2)
675            problemStatus_=-2;
676          continue; // skip rest of iteration
677        } else {
678          // recompute dualOut_
679          if (directionOut_<0) {
680            dualOut_ = valueOut_ - upperOut_;
681          } else {
682            dualOut_ = lowerOut_ - valueOut_;
683          }
684        }
685      }
686      if (sequenceIn_>=0) {
687        // normal iteration
688        // update the incoming column
689        double btranAlpha = -alpha_*directionOut_; // for check
690        unpackPacked(rowArray_[1]);
691        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
692        // and update dual weights (can do in parallel - with extra array)
693        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
694                                              rowArray_[2],
695                                              rowArray_[1]);
696        // see if update stable
697#ifdef CLP_DEBUG
698        if ((handler_->logLevel()&32))
699          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
700#endif
701        double checkValue=1.0e-7;
702        // if can't trust much and long way from optimal then relax
703        if (largestPrimalError_>10.0)
704          checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
705        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
706            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
707          handler_->message(CLP_DUAL_CHECK,messages_)
708            <<btranAlpha
709            <<alpha_
710            <<CoinMessageEol;
711          if (factorization_->pivots()) {
712            dualRowPivot_->unrollWeights();
713            problemStatus_=-2; // factorize now
714            rowArray_[0]->clear();
715            rowArray_[1]->clear();
716            columnArray_[0]->clear();
717            returnCode=-2;
718            break;
719          } else {
720            // take on more relaxed criterion
721            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
722                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
723              dualRowPivot_->unrollWeights();
724              // need to reject something
725              char x = isColumn(sequenceOut_) ? 'C' :'R';
726              handler_->message(CLP_SIMPLEX_FLAG,messages_)
727                <<x<<sequenceWithin(sequenceOut_)
728                <<CoinMessageEol;
729              setFlagged(sequenceOut_);
730              lastBadIteration_ = numberIterations_; // say be more cautious
731              rowArray_[0]->clear();
732              rowArray_[1]->clear();
733              columnArray_[0]->clear();
734              continue;
735            }
736          }
737        }
738        // update duals BEFORE replaceColumn so can do updateColumn
739        double objectiveChange=0.0;
740        // do duals first as variables may flip bounds
741        // rowArray_[0] and columnArray_[0] may have flips
742        // so use rowArray_[3] for work array from here on
743        int nswapped = 0;
744        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
745        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
746        if (candidate==-1) 
747          nswapped = updateDualsInDual(rowArray_[0],columnArray_[0],
748                                       rowArray_[2],theta_,
749                                       objectiveChange,false);
750        else 
751          updateDualsInValuesPass(rowArray_[0],columnArray_[0],theta_);
752
753        // which will change basic solution
754        if (nswapped) {
755          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
756          dualRowPivot_->updatePrimalSolution(rowArray_[2],
757                                              1.0,objectiveChange);
758         
759          // recompute dualOut_
760          valueOut_ = solution_[sequenceOut_];
761          if (directionOut_<0) {
762            dualOut_ = valueOut_ - upperOut_;
763          } else {
764            dualOut_ = lowerOut_ - valueOut_;
765          }
766          if(dualOut_<-max(1.0e-12*averagePrimalInfeasibility,1.0e-8)
767             &&factorization_->pivots()&&
768             getStatus(sequenceIn_)!=isFree) {
769            // going backwards - factorize
770            dualRowPivot_->unrollWeights();
771            problemStatus_=-2; // factorize now
772            returnCode=-2;
773            break;
774          }
775        }
776        // amount primal will move
777        double movement = -dualOut_*directionOut_/alpha_;
778        // so objective should increase by fabs(dj)*movement
779        // but we already have objective change - so check will be good
780        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
781#ifdef CLP_DEBUG
782          if (handler_->logLevel()&32)
783            printf("movement %g, swap change %g, rest %g  * %g\n",
784                   objectiveChange+fabs(movement*dualIn_),
785                   objectiveChange,movement,dualIn_);
786#endif
787          if(factorization_->pivots()>5) {
788            // going backwards - factorize
789            dualRowPivot_->unrollWeights();
790            problemStatus_=-2; // factorize now
791            returnCode=-2;
792            break;
793          }
794        }
795        assert(fabs(dualOut_)<1.0e50);
796        // if stable replace in basis
797        int updateStatus = factorization_->replaceColumn(rowArray_[2],
798                                                         pivotRow_,
799                                                         alpha_);
800        // if no pivots, bad update but reasonable alpha - take and invert
801        if (updateStatus==2&&
802                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
803          updateStatus=4;
804        if (updateStatus==1||updateStatus==4) {
805          // slight error
806          if (factorization_->pivots()>5||updateStatus==4) {
807            problemStatus_=-2; // factorize now
808            returnCode=-3;
809          }
810        } else if (updateStatus==2) {
811          // major error
812          dualRowPivot_->unrollWeights();
813          // later we may need to unwind more e.g. fake bounds
814          if (factorization_->pivots()) {
815            problemStatus_=-2; // factorize now
816            returnCode=-2;
817            break;
818          } else {
819            // need to reject something
820            char x = isColumn(sequenceOut_) ? 'C' :'R';
821            handler_->message(CLP_SIMPLEX_FLAG,messages_)
822              <<x<<sequenceWithin(sequenceOut_)
823              <<CoinMessageEol;
824            setFlagged(sequenceOut_);
825            lastBadIteration_ = numberIterations_; // say be more cautious
826            rowArray_[0]->clear();
827            rowArray_[1]->clear();
828            columnArray_[0]->clear();
829            // make sure dual feasible
830            // look at all rows and columns
831            double objectiveChange=0.0;
832            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
833                              0.0,objectiveChange,true);
834            continue;
835          }
836        } else if (updateStatus==3) {
837          // out of memory
838          // increase space if not many iterations
839          if (factorization_->pivots()<
840              0.5*factorization_->maximumPivots()&&
841              factorization_->pivots()<200)
842            factorization_->areaFactor(
843                                       factorization_->areaFactor() * 1.1);
844          problemStatus_=-2; // factorize now
845        } else if (updateStatus==5) {
846          problemStatus_=-2; // factorize now
847        }
848        // update primal solution
849        if (theta_<0.0&&candidate==-1) {
850#ifdef CLP_DEBUG
851          if (handler_->logLevel()&32)
852            printf("negative theta %g\n",theta_);
853#endif
854          theta_=0.0;
855        }
856        // do actual flips
857        flipBounds(rowArray_[0],columnArray_[0],theta_);
858        //rowArray_[1]->expand();
859        dualRowPivot_->updatePrimalSolution(rowArray_[1],
860                                            movement,
861                                            objectiveChange);
862#ifdef CLP_DEBUG
863        double oldobj=objectiveValue_;
864#endif
865        // modify dualout
866        dualOut_ /= alpha_;
867        dualOut_ *= -directionOut_;
868        //setStatus(sequenceIn_,basic);
869        dj_[sequenceIn_]=0.0;
870        double oldValue=valueIn_;
871        if (directionIn_==-1) {
872          // as if from upper bound
873          valueIn_ = upperIn_+dualOut_;
874        } else {
875          // as if from lower bound
876          valueIn_ = lowerIn_+dualOut_;
877        }
878        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
879        // outgoing
880        // set dj to zero unless values pass
881        if (directionOut_>0) {
882          valueOut_ = lowerOut_;
883          if (candidate==-1)
884            dj_[sequenceOut_] = theta_;
885        } else {
886          valueOut_ = upperOut_;
887          if (candidate==-1)
888            dj_[sequenceOut_] = -theta_;
889        }
890        solution_[sequenceOut_]=valueOut_;
891        int whatNext=housekeeping(objectiveChange);
892        // and set bounds correctly
893        originalBound(sequenceIn_); 
894        changeBound(sequenceOut_);
895#ifdef CLP_DEBUG
896        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
897          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
898#endif
899        if (whatNext==1||candidate==-2) {
900          problemStatus_ =-2; // refactorize
901        } else if (whatNext==2) {
902          // maximum iterations or equivalent
903          problemStatus_= 3;
904          returnCode=3;
905          break;
906        }
907      } else {
908        // no incoming column is valid
909#ifdef CLP_DEBUG
910        if (handler_->logLevel()&32)
911          printf("** no column pivot\n");
912#endif
913        if (factorization_->pivots()<5) {
914          problemStatus_=-4; //say looks infeasible
915          // create ray anyway
916          delete [] ray_;
917          ray_ = new double [ numberRows_];
918          rowArray_[0]->expand(); // in case packed
919          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
920        }
921        rowArray_[0]->clear();
922        columnArray_[0]->clear();
923        returnCode=1;
924        break;
925      }
926    } else {
927      // no pivot row
928#ifdef CLP_DEBUG
929      if (handler_->logLevel()&32)
930        printf("** no row pivot\n");
931#endif
932      if (!factorization_->pivots()) {
933        // may have crept through - so may be optimal
934        // check any flagged variables
935        int iRow;
936        for (iRow=0;iRow<numberRows_;iRow++) {
937          int iPivot=pivotVariable_[iRow];
938          if (flagged(iPivot))
939            break;
940        }
941        if (numberFake_||numberDualInfeasibilities_) {
942          // may be dual infeasible
943          problemStatus_=-5;
944        } else {
945          if (iRow<numberRows_) {
946#ifdef CLP_DEBUG
947            std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
948            printf("Probably infeasible - pivot was %g\n",alpha_);
949#endif
950            if (fabs(alpha_)<1.0e-4) {
951              problemStatus_=1;
952            } else {
953#ifdef CLP_DEBUG
954              abort();
955#endif
956            }
957          } else {
958            problemStatus_=0;
959          }
960        }
961      } else {
962        problemStatus_=-3;
963      }
964      returnCode=0;
965      break;
966    }
967  }
968  if (givenDuals) {
969    memcpy(givenDuals,dj_,(numberRows_+numberColumns_)*sizeof(double));
970    // get rid of any values pass array
971    delete [] candidateList;
972  }
973  delete [] dubiousWeights;
974  return returnCode;
975}
976/* The duals are updated by the given arrays.
977   Returns number of infeasibilities.
978   rowArray and columnarray will have flipped
979   The output vector has movement (row length array) */
980int
981ClpPdco::updateDualsInDual(CoinIndexedVector * rowArray,
982                                  CoinIndexedVector * columnArray,
983                                  CoinIndexedVector * outputArray,
984                                  double theta,
985                                  double & objectiveChange,
986                                  bool fullRecompute)
987{
988 
989  outputArray->clear();
990 
991 
992  int numberInfeasibilities=0;
993  int numberRowInfeasibilities=0;
994 
995  // use a tighter tolerance except for all being okay
996  double tolerance = dualTolerance_;
997 
998  double changeObj=0.0;
999
1000  // Coding is very similar but we can save a bit by splitting
1001  // Do rows
1002  if (!fullRecompute) {
1003    int i;
1004    double * reducedCost = djRegion(0);
1005    const double * lower = lowerRegion(0);
1006    const double * upper = upperRegion(0);
1007    const double * cost = costRegion(0);
1008    double * work;
1009    int number;
1010    int * which;
1011    assert(rowArray->packedMode());
1012    work = rowArray->denseVector();
1013    number = rowArray->getNumElements();
1014    which = rowArray->getIndices();
1015    for (i=0;i<number;i++) {
1016      int iSequence = which[i];
1017      double alphaI = work[i];
1018      work[i]=0.0;
1019     
1020      Status status = getStatus(iSequence+numberColumns_);
1021      // more likely to be at upper bound ?
1022      if (status==atUpperBound) {
1023
1024        double value = reducedCost[iSequence]-theta*alphaI;
1025        reducedCost[iSequence]=value;
1026
1027        if (value>tolerance) {
1028          // to lower bound (if swap)
1029          which[numberInfeasibilities++]=iSequence;
1030          double movement = lower[iSequence]-upper[iSequence];
1031#ifdef CLP_DEBUG
1032          if ((handler_->logLevel()&32))
1033            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1034                   0,iSequence,value,alphaI,movement);
1035#endif
1036          changeObj += movement*cost[iSequence];
1037          outputArray->quickAdd(iSequence,-movement);
1038        }
1039      } else if (status==atLowerBound) {
1040
1041        double value = reducedCost[iSequence]-theta*alphaI;
1042        reducedCost[iSequence]=value;
1043
1044        if (value<-tolerance) {
1045          // to upper bound
1046          which[numberInfeasibilities++]=iSequence;
1047          double movement = upper[iSequence] - lower[iSequence];
1048#ifdef CLP_DEBUG
1049          if ((handler_->logLevel()&32))
1050            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1051                   0,iSequence,value,alphaI,movement);
1052#endif
1053          changeObj += movement*cost[iSequence];
1054          outputArray->quickAdd(iSequence,-movement);
1055        }
1056      }
1057    }
1058  } else  {
1059    double * solution = solutionRegion(0);
1060    double * reducedCost = djRegion(0);
1061    const double * lower = lowerRegion(0);
1062    const double * upper = upperRegion(0);
1063    const double * cost = costRegion(0);
1064    int * which;
1065    which = rowArray->getIndices();
1066    int iSequence;
1067    for (iSequence=0;iSequence<numberRows_;iSequence++) {
1068      double value = reducedCost[iSequence];
1069     
1070      Status status = getStatus(iSequence+numberColumns_);
1071      // more likely to be at upper bound ?
1072      if (status==atUpperBound) {
1073        double movement=0.0;
1074
1075        if (value>tolerance) {
1076          // to lower bound (if swap)
1077          // put back alpha
1078          which[numberInfeasibilities++]=iSequence;
1079          movement = lower[iSequence]-upper[iSequence];
1080          changeObj += movement*cost[iSequence];
1081          outputArray->quickAdd(iSequence,-movement);
1082        } else if (value>-tolerance) {
1083          // at correct bound but may swap
1084          FakeBound bound = getFakeBound(iSequence+numberColumns_);
1085          if (bound==ClpPdco::upperFake) {
1086            movement = lower[iSequence]-upper[iSequence];
1087            setStatus(iSequence+numberColumns_,atLowerBound);
1088            solution[iSequence] = lower[iSequence];
1089            changeObj += movement*cost[iSequence];
1090            numberFake_--;
1091            setFakeBound(iSequence+numberColumns_,noFake);
1092          }
1093        }
1094      } else if (status==atLowerBound) {
1095        double movement=0.0;
1096
1097        if (value<-tolerance) {
1098          // to upper bound
1099          // put back alpha
1100          which[numberInfeasibilities++]=iSequence;
1101          movement = upper[iSequence] - lower[iSequence];
1102          changeObj += movement*cost[iSequence];
1103          outputArray->quickAdd(iSequence,-movement);
1104        } else if (value<tolerance) {
1105          // at correct bound but may swap
1106          FakeBound bound = getFakeBound(iSequence+numberColumns_);
1107          if (bound==ClpPdco::lowerFake) {
1108            movement = upper[iSequence]-lower[iSequence];
1109            setStatus(iSequence+numberColumns_,atUpperBound);
1110            solution[iSequence] = upper[iSequence];
1111            changeObj += movement*cost[iSequence];
1112            numberFake_--;
1113            setFakeBound(iSequence+numberColumns_,noFake);
1114          }
1115        }
1116      }
1117    }
1118  }
1119
1120  // Do columns
1121  if (!fullRecompute) {
1122    int i;
1123    double * reducedCost = djRegion(1);
1124    const double * lower = lowerRegion(1);
1125    const double * upper = upperRegion(1);
1126    const double * cost = costRegion(1);
1127    double * work;
1128    int number;
1129    int * which;
1130    // set number of infeasibilities in row array
1131    numberRowInfeasibilities=numberInfeasibilities;
1132    rowArray->setNumElements(numberInfeasibilities);
1133    numberInfeasibilities=0;
1134    work = columnArray->denseVector();
1135    number = columnArray->getNumElements();
1136    which = columnArray->getIndices();
1137   
1138    for (i=0;i<number;i++) {
1139      int iSequence = which[i];
1140      double alphaI = work[i];
1141      work[i]=0.0;
1142     
1143      Status status = getStatus(iSequence);
1144      if (status==atLowerBound) {
1145        double value = reducedCost[iSequence]-theta*alphaI;
1146        reducedCost[iSequence]=value;
1147        double movement=0.0;
1148
1149        if (value<-tolerance) {
1150          // to upper bound
1151          which[numberInfeasibilities++]=iSequence;
1152          movement = upper[iSequence] - lower[iSequence];
1153#ifdef CLP_DEBUG
1154          if ((handler_->logLevel()&32))
1155            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1156                   1,iSequence,value,alphaI,movement);
1157#endif
1158          changeObj += movement*cost[iSequence];
1159          matrix_->add(this,outputArray,iSequence,movement);
1160        }
1161      } else if (status==atUpperBound) {
1162        double value = reducedCost[iSequence]-theta*alphaI;
1163        reducedCost[iSequence]=value;
1164        double movement=0.0;
1165
1166        if (value>tolerance) {
1167          // to lower bound (if swap)
1168          which[numberInfeasibilities++]=iSequence;
1169          movement = lower[iSequence]-upper[iSequence];
1170#ifdef CLP_DEBUG
1171          if ((handler_->logLevel()&32))
1172            printf("%d %d, new dj %g, alpha %g, movement %g\n",
1173                   1,iSequence,value,alphaI,movement);
1174#endif
1175          changeObj += movement*cost[iSequence];
1176          matrix_->add(this,outputArray,iSequence,movement);
1177        }
1178      } else if (status==isFree) {
1179        double value = reducedCost[iSequence]-theta*alphaI;
1180        reducedCost[iSequence]=value;
1181      }
1182    }
1183  } else {
1184    double * solution = solutionRegion(1);
1185    double * reducedCost = djRegion(1);
1186    const double * lower = lowerRegion(1);
1187    const double * upper = upperRegion(1);
1188    const double * cost = costRegion(1);
1189    int * which;
1190    // set number of infeasibilities in row array
1191    numberRowInfeasibilities=numberInfeasibilities;
1192    rowArray->setNumElements(numberInfeasibilities);
1193    numberInfeasibilities=0;
1194    which = columnArray->getIndices();
1195    int iSequence;
1196    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
1197      double value = reducedCost[iSequence];
1198     
1199      Status status = getStatus(iSequence);
1200      if (status==atLowerBound) {
1201        double movement=0.0;
1202
1203        if (value<-tolerance) {
1204          // to upper bound
1205          // put back alpha
1206          which[numberInfeasibilities++]=iSequence;
1207          movement = upper[iSequence] - lower[iSequence];
1208          changeObj += movement*cost[iSequence];
1209          matrix_->add(this,outputArray,iSequence,movement);
1210        } else if (value<tolerance) {
1211          // at correct bound but may swap
1212          FakeBound bound = getFakeBound(iSequence);
1213          if (bound==ClpPdco::lowerFake) {
1214            movement = upper[iSequence]-lower[iSequence];
1215            setStatus(iSequence,atUpperBound);
1216            solution[iSequence] = upper[iSequence];
1217            changeObj += movement*cost[iSequence];
1218            numberFake_--;
1219            setFakeBound(iSequence,noFake);
1220          }
1221        }
1222      } else if (status==atUpperBound) {
1223        double movement=0.0;
1224
1225        if (value>tolerance) {
1226          // to lower bound (if swap)
1227          // put back alpha
1228          which[numberInfeasibilities++]=iSequence;
1229          movement = lower[iSequence]-upper[iSequence];
1230          changeObj += movement*cost[iSequence];
1231          matrix_->add(this,outputArray,iSequence,movement);
1232        } else if (value>-tolerance) {
1233          // at correct bound but may swap
1234          FakeBound bound = getFakeBound(iSequence);
1235          if (bound==ClpPdco::upperFake) {
1236            movement = lower[iSequence]-upper[iSequence];
1237            setStatus(iSequence,atLowerBound);
1238            solution[iSequence] = lower[iSequence];
1239            changeObj += movement*cost[iSequence];
1240            numberFake_--;
1241            setFakeBound(iSequence,noFake);
1242          }
1243        }
1244      }
1245    }
1246  }
1247#ifdef CLP_DEBUG
1248  if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0) 
1249    printf("%d fake after full update\n",numberFake_);
1250#endif
1251  // set number of infeasibilities
1252  columnArray->setNumElements(numberInfeasibilities);
1253  numberInfeasibilities += numberRowInfeasibilities;
1254  if (fullRecompute) {
1255    // do actual flips
1256    flipBounds(rowArray,columnArray,0.0);
1257  }
1258  objectiveChange += changeObj;
1259  return numberInfeasibilities;
1260}
1261void
1262ClpPdco::updateDualsInValuesPass(CoinIndexedVector * rowArray,
1263                                  CoinIndexedVector * columnArray,
1264                                        double theta)
1265{
1266 
1267  // use a tighter tolerance except for all being okay
1268  double tolerance = dualTolerance_;
1269 
1270  // Coding is very similar but we can save a bit by splitting
1271  // Do rows
1272  {
1273    int i;
1274    double * reducedCost = djRegion(0);
1275    double * work;
1276    int number;
1277    int * which;
1278    work = rowArray->denseVector();
1279    number = rowArray->getNumElements();
1280    which = rowArray->getIndices();
1281    for (i=0;i<number;i++) {
1282      int iSequence = which[i];
1283      double alphaI = work[i];
1284      double value = reducedCost[iSequence]-theta*alphaI;
1285      work[i]=0.0;
1286      reducedCost[iSequence]=value;
1287     
1288      Status status = getStatus(iSequence+numberColumns_);
1289      // more likely to be at upper bound ?
1290      if (status==atUpperBound) {
1291
1292        if (value>tolerance) 
1293          reducedCost[iSequence]=0.0;
1294      } else if (status==atLowerBound) {
1295
1296        if (value<-tolerance) {
1297          reducedCost[iSequence]=0.0;
1298        }
1299      }
1300    }
1301  }
1302  rowArray->setNumElements(0);
1303
1304  // Do columns
1305  {
1306    int i;
1307    double * reducedCost = djRegion(1);
1308    double * work;
1309    int number;
1310    int * which;
1311    work = columnArray->denseVector();
1312    number = columnArray->getNumElements();
1313    which = columnArray->getIndices();
1314   
1315    for (i=0;i<number;i++) {
1316      int iSequence = which[i];
1317      double alphaI = work[i];
1318      double value = reducedCost[iSequence]-theta*alphaI;
1319      work[i]=0.0;
1320      reducedCost[iSequence]=value;
1321     
1322      Status status = getStatus(iSequence);
1323      if (status==atLowerBound) {
1324        if (value<-tolerance) 
1325          reducedCost[iSequence]=0.0;
1326      } else if (status==atUpperBound) {
1327        if (value>tolerance) 
1328          reducedCost[iSequence]=0.0;
1329      }
1330    }
1331  }
1332  columnArray->setNumElements(0);
1333}
1334/*
1335   Chooses dual pivot row
1336   Would be faster with separate region to scan
1337   and will have this (with square of infeasibility) when steepest
1338   For easy problems we can just choose one of the first rows we look at
1339*/
1340void
1341ClpPdco::dualRow(int alreadyChosen)
1342{
1343  // get pivot row using whichever method it is
1344  int chosenRow=-1;
1345  if (alreadyChosen<0) {
1346    // first see if any free variables and put them in basis
1347    int nextFree = nextSuperBasic();
1348    //nextFree=-1; //off
1349    if (nextFree>=0) {
1350      // unpack vector and find a good pivot
1351      unpack(rowArray_[1],nextFree);
1352      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
1353     
1354      double * work=rowArray_[1]->denseVector();
1355      int number=rowArray_[1]->getNumElements();
1356      int * which=rowArray_[1]->getIndices();
1357      double bestFeasibleAlpha=0.0;
1358      int bestFeasibleRow=-1;
1359      double bestInfeasibleAlpha=0.0;
1360      int bestInfeasibleRow=-1;
1361      int i;
1362     
1363      for (i=0;i<number;i++) {
1364        int iRow = which[i];
1365        double alpha = fabs(work[iRow]);
1366        if (alpha>1.0e-3) {
1367          int iSequence=pivotVariable_[iRow];
1368          double value = solution_[iSequence];
1369          double lower = lower_[iSequence];
1370          double upper = upper_[iSequence];
1371          double infeasibility=0.0;
1372          if (value>upper)
1373            infeasibility = value-upper;
1374          else if (value<lower)
1375            infeasibility = lower-value;
1376          if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
1377            if (!flagged(iSequence)) {
1378              bestInfeasibleAlpha = infeasibility*alpha;
1379              bestInfeasibleRow=iRow;
1380            }
1381          }
1382          if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
1383            bestFeasibleAlpha = alpha;
1384            bestFeasibleRow=iRow;
1385          }
1386        }
1387      }
1388      if (bestInfeasibleRow>=0) 
1389        chosenRow = bestInfeasibleRow;
1390      else if (bestFeasibleAlpha>1.0e-2)
1391        chosenRow = bestFeasibleRow;
1392      if (chosenRow>=0)
1393        pivotRow_=chosenRow;
1394      rowArray_[1]->clear();
1395    } 
1396  } else {
1397    // in values pass
1398    chosenRow=alreadyChosen;
1399    pivotRow_=chosenRow;
1400  }
1401  if (chosenRow<0) 
1402    pivotRow_=dualRowPivot_->pivotRow();
1403
1404  if (pivotRow_>=0) {
1405    sequenceOut_ = pivotVariable_[pivotRow_];
1406    valueOut_ = solution_[sequenceOut_];
1407    lowerOut_ = lower_[sequenceOut_];
1408    upperOut_ = upper_[sequenceOut_];
1409    if (alreadyChosen<0) {
1410      // if we have problems we could try other way and hope we get a
1411      // zero pivot?
1412      if (valueOut_>upperOut_) {
1413        directionOut_ = -1;
1414        dualOut_ = valueOut_ - upperOut_;
1415      } else if (valueOut_<lowerOut_) {
1416        directionOut_ = 1;
1417        dualOut_ = lowerOut_ - valueOut_;
1418      } else {
1419        // odd (could be free) - it's feasible - go to nearest
1420        if (valueOut_-lowerOut_<upperOut_-valueOut_) {
1421          directionOut_ = 1;
1422          dualOut_ = lowerOut_ - valueOut_;
1423        } else {
1424          directionOut_ = -1;
1425          dualOut_ = valueOut_ - upperOut_;
1426        }
1427      }
1428#ifdef CLP_DEBUG
1429      assert(dualOut_>=0.0);
1430#endif
1431    } else {
1432      // in values pass so just use sign of dj
1433      // We don't want to go through any barriers so set dualOut low
1434      // free variables will never be here
1435      dualOut_ = 1.0e-6;
1436      if (dj_[sequenceOut_]>0.0) {
1437        // this will give a -1 in pivot row (as slacks are -1.0)
1438        directionOut_ = 1;
1439      } else {
1440        directionOut_ = -1;
1441      }
1442    }
1443  }
1444  return ;
1445}
1446// Checks if any fake bounds active - if so returns number and modifies
1447// dualBound_ and everything.
1448// Free variables will be left as free
1449// Returns number of bounds changed if >=0
1450// Returns -1 if not initialize and no effect
1451// Fills in changeVector which can be used to see if unbounded
1452// and cost of change vector
1453int
1454ClpPdco::changeBounds(bool initialize,
1455                                 CoinIndexedVector * outputArray,
1456                                 double & changeCost)
1457{ 
1458  numberFake_ = 0;
1459  if (!initialize) {
1460    int numberInfeasibilities;
1461    double newBound;
1462    newBound = 5.0*dualBound_;
1463    numberInfeasibilities=0;
1464    changeCost=0.0;
1465    // put back original bounds and then check
1466    createRim(3);
1467    int iSequence;
1468    // bounds will get bigger - just look at ones at bounds
1469    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1470      double lowerValue=lower_[iSequence];
1471      double upperValue=upper_[iSequence];
1472      double value=solution_[iSequence];
1473      setFakeBound(iSequence,ClpPdco::noFake);
1474      switch(getStatus(iSequence)) {
1475
1476      case basic:
1477      case ClpSimplex::isFixed:
1478        break;
1479      case isFree:
1480      case superBasic:
1481        break;
1482      case atUpperBound:
1483        if (fabs(value-upperValue)>primalTolerance_) 
1484          numberInfeasibilities++;
1485        break;
1486      case atLowerBound:
1487        if (fabs(value-lowerValue)>primalTolerance_) 
1488          numberInfeasibilities++;
1489        break;
1490      }
1491    }
1492    // If dual infeasible then carry on
1493    if (numberInfeasibilities) {
1494      handler_->message(CLP_DUAL_CHECKB,messages_)
1495        <<newBound
1496        <<CoinMessageEol;
1497      int iSequence;
1498      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1499        double lowerValue=lower_[iSequence];
1500        double upperValue=upper_[iSequence];
1501        double newLowerValue;
1502        double newUpperValue;
1503        Status status = getStatus(iSequence);
1504        if (status==atUpperBound||
1505            status==atLowerBound) {
1506          double value = solution_[iSequence];
1507          if (value-lowerValue<=upperValue-value) {
1508            newLowerValue = max(lowerValue,value-0.666667*newBound);
1509            newUpperValue = min(upperValue,newLowerValue+newBound);
1510          } else {
1511            newUpperValue = min(upperValue,value+0.666667*newBound);
1512            newLowerValue = max(lowerValue,newUpperValue-newBound);
1513          }
1514          lower_[iSequence]=newLowerValue;
1515          upper_[iSequence]=newUpperValue;
1516          if (newLowerValue > lowerValue) {
1517            if (newUpperValue < upperValue) {
1518              setFakeBound(iSequence,ClpPdco::bothFake);
1519              numberFake_++;
1520            } else {
1521              setFakeBound(iSequence,ClpPdco::lowerFake);
1522              numberFake_++;
1523            }
1524          } else {
1525            if (newUpperValue < upperValue) {
1526              setFakeBound(iSequence,ClpPdco::upperFake);
1527              numberFake_++;
1528            }
1529          }
1530          if (status==atUpperBound)
1531            solution_[iSequence] = newUpperValue;
1532          else
1533            solution_[iSequence] = newLowerValue;
1534          double movement = solution_[iSequence] - value;
1535          if (movement&&outputArray) {
1536            if (iSequence>=numberColumns_) {
1537              outputArray->quickAdd(iSequence,-movement);
1538              changeCost += movement*cost_[iSequence];
1539            } else {
1540              matrix_->add(this,outputArray,iSequence,movement);
1541              changeCost += movement*cost_[iSequence];
1542            }
1543          }
1544        }
1545      }
1546      dualBound_ = newBound;
1547    } else {
1548      numberInfeasibilities=-1;
1549    }
1550    return numberInfeasibilities;
1551  } else {
1552    int iSequence;
1553     
1554    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
1555      Status status = getStatus(iSequence);
1556      if (status==atUpperBound||
1557          status==atLowerBound) {
1558        double lowerValue=lower_[iSequence];
1559        double upperValue=upper_[iSequence];
1560        double value = solution_[iSequence];
1561        if (lowerValue>-largeValue_||upperValue<largeValue_) {
1562          if (lowerValue-value>-0.5*dualBound_||
1563              upperValue-value<0.5*dualBound_) {
1564            if (fabs(lowerValue-value)<=fabs(upperValue-value)) {
1565              if (upperValue > lowerValue + dualBound_) {
1566                upper_[iSequence]=lowerValue+dualBound_;
1567                setFakeBound(iSequence,ClpPdco::upperFake);
1568                numberFake_++;
1569              }
1570            } else {
1571              if (lowerValue < upperValue - dualBound_) {
1572                lower_[iSequence]=upperValue-dualBound_;
1573                setFakeBound(iSequence,ClpPdco::lowerFake);
1574                numberFake_++;
1575              }
1576            }
1577          } else {
1578            lower_[iSequence]=-0.5*dualBound_;
1579            upper_[iSequence]= 0.5*dualBound_;
1580            setFakeBound(iSequence,ClpPdco::bothFake);
1581            numberFake_++;
1582          }
1583        } else {
1584          // set non basic free variables to fake bounds
1585          // I don't think we should ever get here
1586          assert(!("should not be here"));
1587          lower_[iSequence]=-0.5*dualBound_;
1588          upper_[iSequence]= 0.5*dualBound_;
1589          setFakeBound(iSequence,ClpPdco::bothFake);
1590          numberFake_++;
1591          setStatus(iSequence,atUpperBound);
1592          solution_[iSequence]=0.5*dualBound_;
1593        }
1594      }
1595    }
1596
1597    return 1;
1598  }
1599}
1600/*
1601   Row array has row part of pivot row (as duals so sign may be switched)
1602   Column array has column part.
1603   This chooses pivot column.
1604   Spare array will be needed when we start getting clever.
1605   We will check for basic so spare array will never overflow.
1606   If necessary will modify costs
1607*/
1608void
1609ClpPdco::dualColumn(CoinIndexedVector * rowArray,
1610                           CoinIndexedVector * columnArray,
1611                           CoinIndexedVector * spareArray,
1612                           CoinIndexedVector * spareArray2,
1613                           double acceptablePivot,
1614                           CoinBigIndex * dubiousWeights)
1615{
1616  double * work;
1617  int number;
1618  int * which;
1619  double * reducedCost;
1620  int iSection;
1621
1622  sequenceIn_=-1;
1623  int numberPossiblySwapped=0;
1624  int numberRemaining=0;
1625 
1626  double totalThru=0.0; // for when variables flip
1627
1628  double bestEverPivot=acceptablePivot;
1629  int lastSequence = -1;
1630  double lastPivot=0.0;
1631  double upperTheta;
1632  double newTolerance = dualTolerance_;
1633  //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
1634  // will we need to increase tolerance
1635  bool thisIncrease=false;
1636  // If we think we need to modify costs (not if something from broad sweep)
1637  bool modifyCosts=false;
1638  // Increase in objective due to swapping bounds (may be negative)
1639  double increaseInObjective=0.0;
1640
1641  // use spareArrays to put ones looked at in
1642  // we are going to flip flop between
1643  int iFlip = 0;
1644  // Possible list of pivots
1645  int interesting[2];
1646  // where possible swapped ones are
1647  int swapped[2];
1648  // for zeroing out arrays after
1649  int marker[2][2];
1650  // pivot elements
1651  double * array[2], * spare, * spare2;
1652  // indices
1653  int * indices[2], * index, * index2;
1654  spareArray->clear();
1655  spareArray2->clear();
1656  array[0] = spareArray->denseVector();
1657  indices[0] = spareArray->getIndices();
1658  spare = array[0];
1659  index = indices[0];
1660  array[1] = spareArray2->denseVector();
1661  indices[1] = spareArray2->getIndices();
1662  int i;
1663  const double * lower;
1664  const double * upper;
1665
1666  // initialize lists
1667  for (i=0;i<2;i++) {
1668    interesting[i]=0;
1669    swapped[i]=numberColumns_;
1670    marker[i][0]=0;
1671    marker[i][1]=numberColumns_;
1672  }
1673
1674  /*
1675    First we get a list of possible pivots.  We can also see if the
1676    problem looks infeasible or whether we want to pivot in free variable.
1677    This may make objective go backwards but can only happen a finite
1678    number of times and I do want free variables basic.
1679
1680    Then we flip back and forth.  At the start of each iteration
1681    interesting[iFlip] should have possible candidates and swapped[iFlip]
1682    will have pivots if we decide to take a previous pivot.
1683    At end of each iteration interesting[1-iFlip] should have
1684    candidates if we go through this theta and swapped[1-iFlip]
1685    pivots if we don't go through.
1686
1687    At first we increase theta and see what happens.  We start
1688    theta at a reasonable guess.  If in right area then we do bit by bit.
1689
1690   */
1691
1692  // do first pass to get possibles
1693  // We can also see if infeasible or pivoting on free
1694  double tentativeTheta = 1.0e22;
1695  upperTheta = 1.0e31;
1696  double freePivot = acceptablePivot;
1697
1698  for (iSection=0;iSection<2;iSection++) {
1699
1700    int addSequence;
1701
1702    if (!iSection) {
1703      lower = rowLowerWork_;
1704      upper = rowUpperWork_;
1705      work = rowArray->denseVector();
1706      number = rowArray->getNumElements();
1707      which = rowArray->getIndices();
1708      reducedCost = rowReducedCost_;
1709      addSequence = numberColumns_;
1710    } else {
1711      lower = columnLowerWork_;
1712      upper = columnUpperWork_;
1713      work = columnArray->denseVector();
1714      number = columnArray->getNumElements();
1715      which = columnArray->getIndices();
1716      reducedCost = reducedCostWork_;
1717      addSequence = 0;
1718    }
1719   
1720    for (i=0;i<number;i++) {
1721      int iSequence = which[i];
1722      double alpha;
1723      double oldValue;
1724      double value;
1725      bool keep;
1726
1727      switch(getStatus(iSequence+addSequence)) {
1728         
1729      case basic:
1730      case ClpSimplex::isFixed:
1731        break;
1732      case isFree:
1733      case superBasic:
1734        alpha = work[i];
1735        oldValue = reducedCost[iSequence];
1736        if (oldValue>dualTolerance_) {
1737          keep = true;
1738        } else if (oldValue<-dualTolerance_) {
1739          keep = true;
1740        } else {
1741          if (fabs(alpha)>10.0*acceptablePivot) 
1742            keep = true;
1743          else
1744            keep=false;
1745        }
1746        if (keep) {
1747          // free - choose largest
1748          if (fabs(alpha)>freePivot) {
1749            freePivot=fabs(alpha);
1750            sequenceIn_ = iSequence + addSequence;
1751            theta_=oldValue/alpha;
1752            alpha_=alpha;
1753          }
1754        }
1755        break;
1756      case atUpperBound:
1757        alpha = work[i];
1758        oldValue = reducedCost[iSequence];
1759        value = oldValue-tentativeTheta*alpha;
1760        assert (oldValue<=dualTolerance_*1.0001);
1761        if (value>newTolerance) {
1762          value = oldValue-upperTheta*alpha;
1763          if (value>newTolerance && -alpha>=acceptablePivot) 
1764            upperTheta = (oldValue-newTolerance)/alpha;
1765          // add to list
1766          spare[numberRemaining]=alpha;
1767          index[numberRemaining++]=iSequence+addSequence;
1768        }
1769        break;
1770      case atLowerBound:
1771        alpha = work[i];
1772        oldValue = reducedCost[iSequence];
1773        value = oldValue-tentativeTheta*alpha;
1774        assert (oldValue>=-dualTolerance_*1.0001);
1775        if (value<-newTolerance) {
1776          value = oldValue-upperTheta*alpha;
1777          if (value<-newTolerance && alpha>=acceptablePivot) 
1778            upperTheta = (oldValue+newTolerance)/alpha;
1779          // add to list
1780          spare[numberRemaining]=alpha;
1781          index[numberRemaining++]=iSequence+addSequence;
1782        }
1783        break;
1784      }
1785    }
1786  }
1787  interesting[0]=numberRemaining;
1788  marker[0][0] = numberRemaining;
1789
1790  if (!numberRemaining&&sequenceIn_<0)
1791    return ; // Looks infeasible
1792
1793  if (sequenceIn_>=0) {
1794    // free variable - always choose
1795  } else {
1796
1797    theta_=1.0e50;
1798    // now flip flop between spare arrays until reasonable theta
1799    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1800   
1801    // loops increasing tentative theta until can't go through
1802   
1803    while (tentativeTheta < 1.0e22) {
1804      double thruThis = 0.0;
1805     
1806      double bestPivot=acceptablePivot;
1807      int bestSequence=-1;
1808     
1809      numberPossiblySwapped = numberColumns_;
1810      numberRemaining = 0;
1811
1812      upperTheta = 1.0e50;
1813
1814      spare = array[iFlip];
1815      index = indices[iFlip];
1816      spare2 = array[1-iFlip];
1817      index2 = indices[1-iFlip];
1818
1819      // try 3 different ways
1820      // 1 bias increase by ones with slightly wrong djs
1821      // 2 bias by all
1822      // 3 bias by all - tolerance (doesn't seem very good)
1823#define TRYBIAS 1
1824
1825
1826      double increaseInThis=0.0; //objective increase in this loop
1827     
1828      for (i=0;i<interesting[iFlip];i++) {
1829        int iSequence = index[i];
1830        double alpha = spare[i];
1831        double oldValue = dj_[iSequence];
1832        double value = oldValue-tentativeTheta*alpha;
1833
1834        if (alpha < 0.0) {
1835          //at upper bound
1836          if (value>newTolerance) {
1837            double range = upper_[iSequence] - lower_[iSequence];
1838            thruThis -= range*alpha;
1839#if TRYBIAS==1
1840            if (oldValue>0.0)
1841              increaseInThis -= oldValue*range;
1842#elif TRYBIAS==2
1843            increaseInThis -= oldValue*range;
1844#else
1845            increaseInThis -= (oldValue+dualTolerance_)*range;
1846#endif
1847            // goes on swapped list (also means candidates if too many)
1848            spare2[--numberPossiblySwapped]=alpha;
1849            index2[numberPossiblySwapped]=iSequence;
1850            if (fabs(alpha)>bestPivot) {
1851              bestPivot=fabs(alpha);
1852              bestSequence=numberPossiblySwapped;
1853            }
1854          } else {
1855            value = oldValue-upperTheta*alpha;
1856            if (value>newTolerance && -alpha>=acceptablePivot) 
1857              upperTheta = (oldValue-newTolerance)/alpha;
1858            spare2[numberRemaining]=alpha;
1859            index2[numberRemaining++]=iSequence;
1860          }
1861        } else {
1862          // at lower bound
1863          if (value<-newTolerance) {
1864            double range = upper_[iSequence] - lower_[iSequence];
1865            thruThis += range*alpha;
1866            //?? is this correct - and should we look at good ones
1867#if TRYBIAS==1
1868            if (oldValue<0.0)
1869              increaseInThis += oldValue*range;
1870#elif TRYBIAS==2
1871            increaseInThis += oldValue*range;
1872#else
1873            increaseInThis += (oldValue-dualTolerance_)*range;
1874#endif
1875            // goes on swapped list (also means candidates if too many)
1876            spare2[--numberPossiblySwapped]=alpha;
1877            index2[numberPossiblySwapped]=iSequence;
1878            if (fabs(alpha)>bestPivot) {
1879              bestPivot=fabs(alpha);
1880              bestSequence=numberPossiblySwapped;
1881            }
1882          } else {
1883            value = oldValue-upperTheta*alpha;
1884            if (value<-newTolerance && alpha>=acceptablePivot) 
1885              upperTheta = (oldValue+newTolerance)/alpha;
1886            spare2[numberRemaining]=alpha;
1887            index2[numberRemaining++]=iSequence;
1888          }
1889        }
1890      }
1891      swapped[1-iFlip]=numberPossiblySwapped;
1892      interesting[1-iFlip]=numberRemaining;
1893      marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
1894      marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
1895     
1896      if (totalThru+thruThis>=fabs(dualOut_)||
1897          increaseInObjective+increaseInThis<0.0) {
1898        // We should be pivoting in this batch
1899        // so compress down to this lot
1900        numberRemaining=0;
1901        for (i=numberColumns_-1;i>=swapped[1-iFlip];i--) {
1902          spare[numberRemaining]=spare2[i];
1903          index[numberRemaining++]=index2[i];
1904        }
1905        interesting[iFlip]=numberRemaining;
1906        int iTry;
1907#define MAXTRY 100
1908        // first get ratio with tolerance
1909        for (iTry=0;iTry<MAXTRY;iTry++) {
1910         
1911          upperTheta=1.0e50;
1912          numberPossiblySwapped = numberColumns_;
1913          numberRemaining = 0;
1914
1915          increaseInThis=0.0; //objective increase in this loop
1916
1917          thruThis=0.0;
1918
1919          spare = array[iFlip];
1920          index = indices[iFlip];
1921          spare2 = array[1-iFlip];
1922          index2 = indices[1-iFlip];
1923     
1924          for (i=0;i<interesting[iFlip];i++) {
1925            int iSequence=index[i];
1926            double alpha=spare[i];
1927            double oldValue = dj_[iSequence];
1928            double value = oldValue-upperTheta*alpha;
1929           
1930            if (alpha < 0.0) {
1931              //at upper bound
1932              if (value>newTolerance) {
1933                if (-alpha>=acceptablePivot) {
1934                  upperTheta = (oldValue-newTolerance)/alpha;
1935                  // recompute value and make sure works
1936                  value = oldValue-upperTheta*alpha;
1937                  if (value<0)
1938                    upperTheta *= 1.0 +1.0e-11; // must be large
1939                }
1940              }
1941            } else {
1942              // at lower bound
1943              if (value<-newTolerance) {
1944                if (alpha>=acceptablePivot) {
1945                  upperTheta = (oldValue+newTolerance)/alpha;
1946                  // recompute value and make sure works
1947                  value = oldValue-upperTheta*alpha;
1948                  if (value>0)
1949                    upperTheta *= 1.0 +1.0e-11; // must be large
1950                }
1951              }
1952            }
1953          }
1954          bestPivot=acceptablePivot;
1955          sequenceIn_=-1;
1956          double bestWeight=COIN_DBL_MAX;
1957          double largestPivot=acceptablePivot;
1958          // now choose largest and sum all ones which will go through
1959          //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
1960          for (i=0;i<interesting[iFlip];i++) {
1961            int iSequence=index[i];
1962            double alpha=spare[i];
1963            double value = dj_[iSequence]-upperTheta*alpha;
1964            double badDj=0.0;
1965           
1966            bool addToSwapped=false;
1967           
1968            if (alpha < 0.0) {
1969              //at upper bound
1970              if (value>=0.0) { 
1971                addToSwapped=true;
1972#if TRYBIAS==1
1973                badDj = -max(dj_[iSequence],0.0);
1974#elif TRYBIAS==2
1975                badDj = -dj_[iSequence];
1976#else
1977                badDj = -dj_[iSequence]-dualTolerance_;
1978#endif
1979              }
1980            } else {
1981              // at lower bound
1982              if (value<=0.0) {
1983                addToSwapped=true;
1984#if TRYBIAS==1
1985                badDj = min(dj_[iSequence],0.0);
1986#elif TRYBIAS==2
1987                badDj = dj_[iSequence];
1988#else
1989                badDj = dj_[iSequence]-dualTolerance_;
1990#endif
1991              }
1992            }
1993            if (!addToSwapped) {
1994              // add to list of remaining
1995              spare2[numberRemaining]=alpha;
1996              index2[numberRemaining++]=iSequence;
1997            } else {
1998              // add to list of swapped
1999              spare2[--numberPossiblySwapped]=alpha;
2000              index2[numberPossiblySwapped]=iSequence;
2001              // select if largest pivot
2002              bool take=false;
2003              double absAlpha = fabs(alpha);
2004              // User could do anything to break ties here
2005              double weight;
2006              if (dubiousWeights)
2007                weight=dubiousWeights[iSequence];
2008              else
2009                weight=1.0;
2010              weight += CoinDrand48()*1.0e-2;
2011              if (absAlpha>2.0*bestPivot) {
2012                take=true;
2013              } else if (absAlpha>0.5*largestPivot) {
2014                // could multiply absAlpha and weight
2015                if (weight*bestPivot<bestWeight*absAlpha)
2016                  take=true;
2017              }
2018              if (take) {
2019                sequenceIn_ = numberPossiblySwapped;
2020                bestPivot =  absAlpha;
2021                theta_ = dj_[iSequence]/alpha;
2022                largestPivot = max(largestPivot,bestPivot);
2023                bestWeight = weight;
2024                //printf(" taken seq %d alpha %g weight %d\n",
2025                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2026              } else {
2027                //printf(" not taken seq %d alpha %g weight %d\n",
2028                //   iSequence,absAlpha,dubiousWeights[iSequence]);
2029              }
2030              double range = upper[iSequence] - lower[iSequence];
2031              thruThis += range*fabs(alpha);
2032              increaseInThis += badDj*range;
2033            }
2034          }
2035          swapped[1-iFlip]=numberPossiblySwapped;
2036          interesting[1-iFlip]=numberRemaining;
2037          marker[1-iFlip][0]= max(marker[1-iFlip][0],numberRemaining);
2038          marker[1-iFlip][1]= min(marker[1-iFlip][1],numberPossiblySwapped);
2039          // If we stop now this will be increase in objective (I think)
2040          double increase = (fabs(dualOut_)-totalThru)*theta_;
2041          increase += increaseInObjective;
2042          if (theta_<0.0)
2043            thruThis += fabs(dualOut_); // force using this one
2044          if (increaseInObjective<0.0&&increase<0.0&&lastSequence>=0) {
2045            // back
2046            // We may need to be more careful - we could do by
2047            // switch so we always do fine grained?
2048            bestPivot=0.0;
2049          } else {
2050            // add in
2051            totalThru += thruThis;
2052            increaseInObjective += increaseInThis;
2053          }
2054          if (bestPivot<0.1*bestEverPivot&&
2055              bestEverPivot>1.0e-6&&
2056              (bestPivot<1.0e-3||totalThru*2.0>fabs(dualOut_))) {
2057            // back to previous one
2058            sequenceIn_=lastSequence;
2059            // swap regions
2060            iFlip = 1-iFlip;
2061            break;
2062          } else if (sequenceIn_==-1&&upperTheta>largeValue_) {
2063            if (lastPivot>acceptablePivot) {
2064              // back to previous one
2065              sequenceIn_=lastSequence;
2066              // swap regions
2067              iFlip = 1-iFlip;
2068            } else {
2069              // can only get here if all pivots too small
2070            }
2071            break;
2072          } else if (totalThru>=fabs(dualOut_)) {
2073            modifyCosts=true; // fine grain - we can modify costs
2074            break; // no point trying another loop
2075          } else {
2076            lastSequence=sequenceIn_;
2077            if (bestPivot>bestEverPivot)
2078              bestEverPivot=bestPivot;
2079            iFlip = 1 -iFlip;
2080            modifyCosts=true; // fine grain - we can modify costs
2081          }
2082        }
2083        if (iTry==MAXTRY)
2084          iFlip = 1-iFlip; // flip back
2085        break;
2086      } else {
2087        // skip this lot
2088        if (bestPivot>1.0e-3||bestPivot>bestEverPivot) {
2089          bestEverPivot=bestPivot;
2090          lastSequence=bestSequence;
2091        } else {
2092          // keep old swapped
2093          memcpy(array[1-iFlip]+swapped[iFlip],
2094                 array[iFlip]+swapped[iFlip],
2095                 (numberColumns_-swapped[iFlip])*sizeof(double));
2096          memcpy(indices[1-iFlip]+swapped[iFlip],
2097                 indices[iFlip]+swapped[iFlip],
2098                 (numberColumns_-swapped[iFlip])*sizeof(int));
2099          marker[1-iFlip][1] = min(marker[1-iFlip][1],swapped[iFlip]);
2100          swapped[1-iFlip]=swapped[iFlip];
2101        }
2102        increaseInObjective += increaseInThis;
2103        iFlip = 1 - iFlip; // swap regions
2104        tentativeTheta = 2.0*upperTheta;
2105        totalThru += thruThis;
2106      }
2107    }
2108   
2109    // can get here without sequenceIn_ set but with lastSequence
2110    if (sequenceIn_<0&&lastSequence>=0) {
2111      // back to previous one
2112      sequenceIn_=lastSequence;
2113      // swap regions
2114      iFlip = 1-iFlip;
2115    }
2116   
2117#define MINIMUMTHETA 1.0e-12
2118    // Movement should be minimum for anti-degeneracy - unless
2119    // fixed variable out
2120    double minimumTheta;
2121    if (upperOut_>lowerOut_)
2122      minimumTheta=MINIMUMTHETA;
2123    else
2124      minimumTheta=0.0;
2125    if (sequenceIn_>=0) {
2126      // at this stage sequenceIn_ is just pointer into index array
2127      // flip just so we can use iFlip
2128      iFlip = 1 -iFlip;
2129      spare = array[iFlip];
2130      index = indices[iFlip];
2131      double oldValue;
2132      alpha_ = spare[sequenceIn_];
2133      sequenceIn_ = indices[iFlip][sequenceIn_];
2134      oldValue = dj_[sequenceIn_];
2135      theta_ = oldValue/alpha_;
2136      if (theta_<minimumTheta) {
2137        // can't pivot to zero
2138        if (oldValue-minimumTheta*alpha_>=-dualTolerance_) {
2139          theta_=minimumTheta;
2140        } else if (oldValue-minimumTheta*alpha_>=-newTolerance) {
2141          theta_=minimumTheta;
2142          thisIncrease=true;
2143        } else {
2144          theta_=(oldValue+newTolerance)/alpha_;
2145          assert(theta_>=0.0);
2146          thisIncrease=true;
2147        } 
2148      }
2149      // may need to adjust costs so all dual feasible AND pivoted is exactly 0
2150      //int costOffset = numberRows_+numberColumns_;
2151      if (modifyCosts) {
2152        int i;
2153        for (i=numberColumns_-1;i>=swapped[iFlip];i--) {
2154          int iSequence=index[i];
2155          double alpha=spare[i];
2156          double value = dj_[iSequence]-theta_*alpha;
2157           
2158          // can't be free here
2159         
2160          if (alpha < 0.0) {
2161            //at upper bound
2162            if (value>dualTolerance_) {
2163              thisIncrease=true;
2164#define MODIFYCOST 2
2165#if MODIFYCOST
2166              // modify cost to hit new tolerance
2167              double modification = alpha*theta_-dj_[iSequence]
2168                +newTolerance;
2169              //modification = min(modification,dualTolerance_);
2170              //assert (fabs(modification)<1.0e-7);
2171              dj_[iSequence] += modification;
2172              cost_[iSequence] +=  modification;
2173              //cost_[iSequence+costOffset] += modification; // save change
2174#endif
2175            }
2176          } else {
2177            // at lower bound
2178            if (-value>dualTolerance_) {
2179              thisIncrease=true;
2180#if MODIFYCOST
2181              // modify cost to hit new tolerance
2182              double modification = alpha*theta_-dj_[iSequence]
2183                -newTolerance;
2184              //modification = max(modification,-dualTolerance_);
2185              //assert (fabs(modification)<1.0e-7);
2186              dj_[iSequence] += modification;
2187              cost_[iSequence] +=  modification;
2188              //cost_[iSequence+costOffset] += modification; // save change
2189#endif
2190            }
2191          }
2192        }
2193      }
2194    }
2195  }
2196
2197  if (sequenceIn_>=0) {
2198    lowerIn_ = lower_[sequenceIn_];
2199    upperIn_ = upper_[sequenceIn_];
2200    valueIn_ = solution_[sequenceIn_];
2201    dualIn_ = dj_[sequenceIn_];
2202
2203    if (numberTimesOptimal_) {
2204      // can we adjust cost back closer to original
2205      //*** add coding
2206    }
2207#if MODIFYCOST>1   
2208    // modify cost to hit zero exactly
2209    // so (dualIn_+modification)==theta_*alpha_
2210    double modification = theta_*alpha_-dualIn_;
2211    dualIn_ += modification;
2212    dj_[sequenceIn_]=dualIn_;
2213    cost_[sequenceIn_] += modification;
2214    //int costOffset = numberRows_+numberColumns_;
2215    //cost_[sequenceIn_+costOffset] += modification; // save change
2216    //assert (fabs(modification)<1.0e-6);
2217#ifdef CLP_DEBUG
2218    if ((handler_->logLevel()&32)&&fabs(modification)>1.0e-15)
2219      printf("exact %d new cost %g, change %g\n",sequenceIn_,
2220             cost_[sequenceIn_],modification);
2221#endif
2222#endif
2223   
2224    if (alpha_<0.0) {
2225      // as if from upper bound
2226      directionIn_=-1;
2227      upperIn_=valueIn_;
2228    } else {
2229      // as if from lower bound
2230      directionIn_=1;
2231      lowerIn_=valueIn_;
2232    }
2233  }
2234  //if (thisIncrease)
2235  //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
2236
2237  // clear arrays
2238
2239  for (i=0;i<2;i++) {
2240    memset(array[i],0,marker[i][0]*sizeof(double));
2241    memset(array[i]+marker[i][1],0,
2242           (numberColumns_-marker[i][1])*sizeof(double));
2243  }
2244}
2245/* Checks if tentative optimal actually means unbounded
2246   Returns -3 if not, 2 if is unbounded */
2247int 
2248ClpPdco::checkUnbounded(CoinIndexedVector * ray,
2249                                   CoinIndexedVector * spare,
2250                                   double changeCost)
2251{
2252  int status=2; // say unbounded
2253  factorization_->updateColumn(spare,ray);
2254  // get reduced cost
2255  int i;
2256  int number=ray->getNumElements();
2257  int * index = ray->getIndices();
2258  double * array = ray->denseVector();
2259  for (i=0;i<number;i++) {
2260    int iRow=index[i];
2261    int iPivot=pivotVariable_[iRow];
2262    changeCost -= cost(iPivot)*array[iRow];
2263  }
2264  double way;
2265  if (changeCost>0.0) {
2266    //try going down
2267    way=1.0;
2268  } else if (changeCost<0.0) {
2269    //try going up
2270    way=-1.0;
2271  } else {
2272#ifdef CLP_DEBUG
2273    printf("can't decide on up or down\n");
2274#endif
2275    way=0.0;
2276    status=-3;
2277  }
2278  double movement=1.0e10*way; // some largish number
2279  double zeroTolerance = 1.0e-14*dualBound_;
2280  for (i=0;i<number;i++) {
2281    int iRow=index[i];
2282    int iPivot=pivotVariable_[iRow];
2283    double arrayValue = array[iRow];
2284    if (fabs(arrayValue)<zeroTolerance)
2285      arrayValue=0.0;
2286    double newValue=solution(iPivot)+movement*arrayValue;
2287    if (newValue>upper(iPivot)+primalTolerance_||
2288        newValue<lower(iPivot)-primalTolerance_)
2289      status=-3; // not unbounded
2290  }
2291  if (status==2) {
2292    // create ray
2293    delete [] ray_;
2294    ray_ = new double [numberColumns_];
2295    ClpFillN(ray_,numberColumns_,0.0);
2296    for (i=0;i<number;i++) {
2297      int iRow=index[i];
2298      int iPivot=pivotVariable_[iRow];
2299      double arrayValue = array[iRow];
2300      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2301        ray_[iPivot] = way* array[iRow];
2302    }
2303  }
2304  ray->clear();
2305  return status;
2306}
2307/* Checks if finished.  Updates status */
2308void 
2309ClpPdco::statusOfProblemInDual(int & lastCleaned,int type,
2310                                      ClpSimplexProgress * progress,
2311                                      double * givenDuals)
2312{
2313  bool normalType=true;
2314  if (type==2) {
2315    // trouble - restore solution
2316    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2317    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2318           numberRows_*sizeof(double));
2319    memcpy(columnActivityWork_,savedSolution_ ,
2320           numberColumns_*sizeof(double));
2321    forceFactorization_=1; // a bit drastic but ..
2322    changeMade_++; // say something changed
2323  }
2324  int tentativeStatus = problemStatus_;
2325  double changeCost;
2326  bool unflagVariables = true;
2327  if (problemStatus_>-3||factorization_->pivots()) {
2328    // factorize
2329    // later on we will need to recover from singularities
2330    // also we could skip if first time
2331    // save dual weights
2332    dualRowPivot_->saveWeights(this,1);
2333    if (type) {
2334      // is factorization okay?
2335      if (internalFactorize(1)) {
2336        // no - restore previous basis
2337        unflagVariables = false;
2338        assert (type==1);
2339        changeMade_++; // say something changed
2340        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2341        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2342               numberRows_*sizeof(double));
2343        memcpy(columnActivityWork_,savedSolution_ ,
2344               numberColumns_*sizeof(double));
2345        // get correct bounds on all variables
2346        double dummyChangeCost=0.0;
2347        changeBounds(true,rowArray_[2],dummyChangeCost);
2348        // throw away change
2349        int i;
2350        for (i=0;i<4;i++) 
2351          rowArray_[i]->clear();
2352        // need to reject something
2353        char x = isColumn(sequenceOut_) ? 'C' :'R';
2354        handler_->message(CLP_SIMPLEX_FLAG,messages_)
2355          <<x<<sequenceWithin(sequenceOut_)
2356          <<CoinMessageEol;
2357        setFlagged(sequenceOut_);
2358       
2359        forceFactorization_=1; // a bit drastic but ..
2360        type = 2;
2361        //assert (internalFactorize(1)==0);
2362        if (internalFactorize(1)) {
2363          memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2364          memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2365                 numberRows_*sizeof(double));
2366          memcpy(columnActivityWork_,savedSolution_ ,
2367                 numberColumns_*sizeof(double));
2368          // debug
2369          assert (internalFactorize(1)==0);
2370        }
2371      }
2372    }
2373    if (problemStatus_!=-4||factorization_->pivots()>10)
2374      problemStatus_=-3;
2375  }
2376  // at this stage status is -3 or -4 if looks infeasible
2377  // get primal and dual solutions
2378  gutsOfSolution(givenDuals,NULL);
2379  // Double check infeasibility if no action
2380  if (progress->lastIterationNumber(0)==numberIterations_) {
2381    if (dualRowPivot_->looksOptimal()) {
2382      numberPrimalInfeasibilities_ = 0;
2383      sumPrimalInfeasibilities_ = 0.0;
2384    }
2385  }
2386  // Check if looping
2387  int loop;
2388  if (!givenDuals&&type!=2) 
2389    loop = progress->looping();
2390  else
2391    loop=-1;
2392  int situationChanged=0;
2393  if (loop>=0) {
2394    problemStatus_ = loop; //exit if in loop
2395    if (!problemStatus_) {
2396      // declaring victory
2397      numberPrimalInfeasibilities_ = 0;
2398      sumPrimalInfeasibilities_ = 0.0;
2399    }
2400    return;
2401  } else if (loop<-1) {
2402    // something may have changed
2403    gutsOfSolution(NULL,NULL);
2404    situationChanged=1;
2405  }
2406  // really for free variables in
2407  if((progressFlag_&2)!=0) {
2408    situationChanged=2;
2409  }
2410  progressFlag_ = 0; //reset progress flag
2411  handler_->message(CLP_SIMPLEX_STATUS,messages_)
2412    <<numberIterations_<<objectiveValue();
2413  handler_->printing(sumPrimalInfeasibilities_>0.0)
2414    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
2415  handler_->printing(sumDualInfeasibilities_>0.0)
2416    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
2417  handler_->printing(numberDualInfeasibilitiesWithoutFree_
2418                     <numberDualInfeasibilities_)
2419                       <<numberDualInfeasibilitiesWithoutFree_;
2420  handler_->message()<<CoinMessageEol;
2421
2422  // dual bound coming in
2423  //double saveDualBound = dualBound_;
2424  while (problemStatus_<=-3) {
2425    int cleanDuals=0;
2426    if (situationChanged!=0)
2427      cleanDuals=1;
2428    int numberChangedBounds=0;
2429    int doOriginalTolerance=0;
2430    if ( lastCleaned==numberIterations_)
2431      doOriginalTolerance=1;
2432    // check optimal
2433    // give code benefit of doubt
2434    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
2435        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2436      // say optimal (with these bounds etc)
2437      numberDualInfeasibilities_ = 0;
2438      sumDualInfeasibilities_ = 0.0;
2439      numberPrimalInfeasibilities_ = 0;
2440      sumPrimalInfeasibilities_ = 0.0;
2441    }
2442    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
2443    if (dualFeasible()||problemStatus_==-4) {
2444      progress->modifyObjective(objectiveValue_
2445                               -sumDualInfeasibilities_*dualBound_);
2446      if (primalFeasible()) {
2447        normalType=false;
2448        // may be optimal - or may be bounds are wrong
2449        handler_->message(CLP_DUAL_BOUNDS,messages_)
2450          <<dualBound_
2451          <<CoinMessageEol;
2452        // save solution in case unbounded
2453        ClpDisjointCopyN(columnActivityWork_,numberColumns_,
2454                          columnArray_[0]->denseVector());
2455        ClpDisjointCopyN(rowActivityWork_,numberRows_,
2456                          rowArray_[2]->denseVector());
2457        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
2458        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
2459          //looks optimal - do we need to reset tolerance
2460          if (perturbation_==101) {
2461            perturbation_=102; // stop any perturbations
2462            cleanDuals=1;
2463            createRim(4);
2464          }
2465          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
2466            doOriginalTolerance=2;
2467            numberTimesOptimal_++;
2468            changeMade_++; // say something changed
2469            if (numberTimesOptimal_==1) {
2470              dualTolerance_ = dblParam_[ClpDualTolerance];
2471              // better to have small tolerance even if slower
2472              factorization_->zeroTolerance(1.0e-15);
2473            } else {
2474              dualTolerance_ = dblParam_[ClpDualTolerance];
2475              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
2476            }
2477            cleanDuals=2; // If nothing changed optimal else primal
2478          } else {
2479            problemStatus_=0; // optimal
2480            if (lastCleaned<numberIterations_) {
2481              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
2482                <<CoinMessageEol;
2483            }
2484          }
2485        } else {
2486          cleanDuals=1;
2487          if (doOriginalTolerance==1) {
2488            // check unbounded
2489            // find a variable with bad dj
2490            int iSequence;
2491            int iChosen=-1;
2492            double largest = 100.0*primalTolerance_;
2493            for (iSequence=0;iSequence<numberRows_+numberColumns_;
2494                 iSequence++) {
2495              double djValue = dj_[iSequence];
2496              double originalLo = originalLower(iSequence);
2497              double originalUp = originalUpper(iSequence);
2498              if (fabs(djValue)>fabs(largest)) {
2499                if (getStatus(iSequence)!=basic) {
2500                  if (djValue>0&&originalLo<-1.0e20) {
2501                    if (djValue>fabs(largest)) {
2502                      largest=djValue;
2503                      iChosen=iSequence;
2504                    }
2505                  } else if (djValue<0&&originalUp>1.0e20) {
2506                    if (-djValue>fabs(largest)) {
2507                      largest=djValue;
2508                      iChosen=iSequence;
2509                    }
2510                  } 
2511                }
2512              }
2513            }
2514            if (iChosen>=0) {
2515              int iSave=sequenceIn_;
2516              sequenceIn_=iChosen;
2517              unpack(rowArray_[1]);
2518              sequenceIn_ = iSave;
2519              // if dual infeasibilities then must be free vector so add in dual
2520              if (numberDualInfeasibilities_) {
2521                if (fabs(changeCost)>1.0e-5)
2522                  printf("Odd free/unbounded combo\n");
2523                changeCost += cost_[iChosen];
2524              }
2525              problemStatus_ = checkUnbounded(rowArray_[1],rowArray_[0],
2526                                              changeCost);
2527              rowArray_[1]->clear();
2528            } else {
2529              problemStatus_=-3;
2530            }
2531            if (problemStatus_==2&&perturbation_==101) {
2532              perturbation_=102; // stop any perturbations
2533              cleanDuals=1;
2534              createRim(4);
2535              problemStatus_=-1;
2536            }
2537            if (problemStatus_==2) {
2538              // it is unbounded - restore solution
2539              // but first add in changes to non-basic
2540              int iColumn;
2541              double * original = columnArray_[0]->denseVector();
2542              for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2543                if(getColumnStatus(iColumn)!= basic)
2544                  ray_[iColumn] += 
2545                    columnActivityWork_[iColumn]-original[iColumn];
2546                columnActivityWork_[iColumn] = original[iColumn];
2547              }
2548              ClpDisjointCopyN(rowArray_[2]->denseVector(),numberRows_,
2549                                rowActivityWork_);
2550            }
2551          } else {
2552            doOriginalTolerance=2;
2553            rowArray_[0]->clear();
2554          }
2555        }
2556        ClpFillN(columnArray_[0]->denseVector(),numberColumns_,0.0);
2557        ClpFillN(rowArray_[2]->denseVector(),numberRows_,0.0);
2558      } 
2559      if (problemStatus_==-4) {
2560        // may be infeasible - or may be bounds are wrong
2561        numberChangedBounds=changeBounds(false,NULL,changeCost);
2562        if (perturbation_==101) {
2563          perturbation_=102; // stop any perturbations
2564          cleanDuals=1;
2565          numberChangedBounds=1;
2566          createRim(4);
2567        }
2568        if (numberChangedBounds<=0||dualBound_>1.0e20||
2569            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
2570          problemStatus_=1; // infeasible
2571        } else {
2572          normalType=false;
2573          problemStatus_=-1; //iterate
2574          cleanDuals=1;
2575          if (numberChangedBounds<=0)
2576            doOriginalTolerance=2;
2577          // and delete ray which has been created
2578          delete [] ray_;
2579          ray_ = NULL;
2580        }
2581
2582      }
2583    } else {
2584      cleanDuals=1;
2585    }
2586    if (problemStatus_<0) {
2587      if (doOriginalTolerance==2) {
2588        // put back original tolerance
2589        lastCleaned=numberIterations_;
2590        handler_->message(CLP_DUAL_ORIGINAL,messages_)
2591          <<CoinMessageEol;
2592
2593        perturbation_=102; // stop any perturbations
2594#if 0
2595        double * xcost = new double[numberRows_+numberColumns_];
2596        double * xlower = new double[numberRows_+numberColumns_];
2597        double * xupper = new double[numberRows_+numberColumns_];
2598        double * xdj = new double[numberRows_+numberColumns_];
2599        double * xsolution = new double[numberRows_+numberColumns_];
2600        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
2601        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
2602        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
2603        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
2604        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
2605#endif
2606        createRim(4);
2607        // make sure duals are current
2608        computeDuals(givenDuals);
2609        checkDualSolution();
2610#if 0
2611        int i;
2612        for (i=0;i<numberRows_+numberColumns_;i++) {
2613          if (cost_[i]!=xcost[i])
2614            printf("** %d old cost %g new %g sol %g\n",
2615                   i,xcost[i],cost_[i],solution_[i]);
2616          if (lower_[i]!=xlower[i])
2617            printf("** %d old lower %g new %g sol %g\n",
2618                   i,xlower[i],lower_[i],solution_[i]);
2619          if (upper_[i]!=xupper[i])
2620            printf("** %d old upper %g new %g sol %g\n",
2621                   i,xupper[i],upper_[i],solution_[i]);
2622          if (dj_[i]!=xdj[i])
2623            printf("** %d old dj %g new %g sol %g\n",
2624                   i,xdj[i],dj_[i],solution_[i]);
2625          if (solution_[i]!=xsolution[i])
2626            printf("** %d old solution %g new %g sol %g\n",
2627                   i,xsolution[i],solution_[i],solution_[i]);
2628        }
2629        //delete [] xcost;
2630        //delete [] xupper;
2631        //delete [] xlower;
2632        //delete [] xdj;
2633        //delete [] xsolution;
2634#endif
2635        // put back bounds as they were if was optimal
2636        if (doOriginalTolerance==2) {
2637          changeMade_++; // say something changed
2638          changeBounds(true,NULL,changeCost);
2639          cleanDuals=2;
2640          //cleanDuals=1;
2641        }
2642#if 0
2643        //int i;
2644        for (i=0;i<numberRows_+numberColumns_;i++) {
2645          if (cost_[i]!=xcost[i])
2646            printf("** %d old cost %g new %g sol %g\n",
2647                   i,xcost[i],cost_[i],solution_[i]);
2648          if (lower_[i]!=xlower[i])
2649            printf("** %d old lower %g new %g sol %g\n",
2650                   i,xlower[i],lower_[i],solution_[i]);
2651          if (upper_[i]!=xupper[i])
2652            printf("** %d old upper %g new %g sol %g\n",
2653                   i,xupper[i],upper_[i],solution_[i]);
2654          if (dj_[i]!=xdj[i])
2655            printf("** %d old dj %g new %g sol %g\n",
2656                   i,xdj[i],dj_[i],solution_[i]);
2657          if (solution_[i]!=xsolution[i])
2658            printf("** %d old solution %g new %g sol %g\n",
2659                   i,xsolution[i],solution_[i],solution_[i]);
2660        }
2661        delete [] xcost;
2662        delete [] xupper;
2663        delete [] xlower;
2664        delete [] xdj;
2665        delete [] xsolution;
2666#endif
2667      }
2668      if (cleanDuals==1||(cleanDuals==2&&!numberDualInfeasibilities_)) {
2669        // make sure dual feasible
2670        // look at all rows and columns
2671        rowArray_[0]->clear();
2672        columnArray_[0]->clear();
2673        double objectiveChange=0.0;
2674#if 0
2675        double * xcost = new double[numberRows_+numberColumns_];
2676        double * xlower = new double[numberRows_+numberColumns_];
2677        double * xupper = new double[numberRows_+numberColumns_];
2678        double * xdj = new double[numberRows_+numberColumns_];
2679        double * xsolution = new double[numberRows_+numberColumns_];
2680        memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
2681        memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
2682        memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
2683        memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
2684        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
2685#endif
2686        updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
2687                          0.0,objectiveChange,true);
2688#if 0
2689        int i;
2690        for (i=0;i<numberRows_+numberColumns_;i++) {
2691          if (cost_[i]!=xcost[i])
2692            printf("** %d old cost %g new %g sol %g\n",
2693                   i,xcost[i],cost_[i],solution_[i]);
2694          if (lower_[i]!=xlower[i])
2695            printf("** %d old lower %g new %g sol %g\n",
2696                   i,xlower[i],lower_[i],solution_[i]);
2697          if (upper_[i]!=xupper[i])
2698            printf("** %d old upper %g new %g sol %g\n",
2699                   i,xupper[i],upper_[i],solution_[i]);
2700          if (dj_[i]!=xdj[i])
2701            printf("** %d old dj %g new %g sol %g\n",
2702                   i,xdj[i],dj_[i],solution_[i]);
2703          if (solution_[i]!=xsolution[i])
2704            printf("** %d old solution %g new %g sol %g\n",
2705                   i,xsolution[i],solution_[i],solution_[i]);
2706        }
2707        delete [] xcost;
2708        delete [] xupper;
2709        delete [] xlower;
2710        delete [] xdj;
2711        delete [] xsolution;
2712#endif
2713        // for now - recompute all
2714        gutsOfSolution(NULL,NULL);
2715        //assert(numberDualInfeasibilitiesWithoutFree_==0);
2716
2717        if (numberDualInfeasibilities_||situationChanged==2) 
2718          problemStatus_=-1; // carry on as normal
2719        situationChanged=0;
2720      } else {
2721        // iterate
2722        if (cleanDuals!=2) 
2723          problemStatus_=-1;
2724        else
2725          problemStatus_ = 10; // try primal
2726      }
2727    }
2728  }
2729  if (type==0||type==1) {
2730    if (!type) {
2731      // create save arrays
2732      delete [] saveStatus_;
2733      delete [] savedSolution_;
2734      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
2735      savedSolution_ = new double [numberRows_+numberColumns_];
2736    }
2737    // save arrays
2738    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
2739    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
2740           numberRows_*sizeof(double));
2741    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
2742  }
2743
2744  // restore weights (if saved) - also recompute infeasibility list
2745  if (tentativeStatus>-3) 
2746    dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
2747  else
2748    dualRowPivot_->saveWeights(this,3);
2749  // unflag all variables (we may want to wait a bit?)
2750  if (tentativeStatus!=-2&&unflagVariables) {
2751    int iRow;
2752    for (iRow=0;iRow<numberRows_;iRow++) {
2753      int iPivot=pivotVariable_[iRow];
2754      clearFlagged(iPivot);
2755    }
2756  }
2757  // see if cutoff reached
2758  double limit = 0.0;
2759  getDblParam(ClpDualObjectiveLimit, limit);
2760  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
2761           optimizationDirection_*limit&&
2762           !numberAtFakeBound()&&!numberDualInfeasibilities_) {
2763    problemStatus_=1;
2764    secondaryStatus_ = 1; // and say was on cutoff
2765  }
2766  if (problemStatus_<0&&!changeMade_) {
2767    problemStatus_=4; // unknown
2768  }
2769  lastGoodIteration_ = numberIterations_;
2770#if 0
2771  double thisObj = progress->lastObjective(0);
2772  double lastObj = progress->lastObjective(1);
2773  if (lastObj>thisObj+1.0e-6*max(fabs(thisObj),fabs(lastObj))+1.0e-8
2774      &&givenDuals==NULL) {
2775    int maxFactor = factorization_->maximumPivots();
2776    if (maxFactor>10) {
2777      if (forceFactorization_<0)
2778        forceFactorization_= maxFactor;
2779      forceFactorization_ = max (1,(forceFactorization_>>1));
2780      printf("Reducing factorization frequency\n");
2781    }
2782  }
2783#endif
2784}
2785/* While updateDualsInDual sees what effect is of flip
2786   this does actual flipping.
2787   If change >0.0 then value in array >0.0 => from lower to upper
2788*/
2789void 
2790ClpPdco::flipBounds(CoinIndexedVector * rowArray,
2791                  CoinIndexedVector * columnArray,
2792                  double change)
2793{
2794  int number;
2795  int * which;
2796 
2797  int iSection;
2798
2799  for (iSection=0;iSection<2;iSection++) {
2800    int i;
2801    double * solution = solutionRegion(iSection);
2802    double * lower = lowerRegion(iSection);
2803    double * upper = upperRegion(iSection);
2804    int addSequence;
2805    if (!iSection) {
2806      number = rowArray->getNumElements();
2807      which = rowArray->getIndices();
2808      addSequence = numberColumns_;
2809    } else {
2810      number = columnArray->getNumElements();
2811      which = columnArray->getIndices();
2812      addSequence = 0;
2813    }
2814   
2815    for (i=0;i<number;i++) {
2816      int iSequence = which[i];
2817      Status status = getStatus(iSequence+addSequence);
2818
2819      switch(status) {
2820
2821      case basic:
2822      case isFree:
2823      case superBasic:
2824      case ClpSimplex::isFixed:
2825        break;
2826      case atUpperBound:
2827        // to lower bound
2828        setStatus(iSequence+addSequence,atLowerBound);
2829        solution[iSequence] = lower[iSequence];
2830        break;
2831      case atLowerBound:
2832        // to upper bound
2833        setStatus(iSequence+addSequence,atUpperBound);
2834        solution[iSequence] = upper[iSequence];
2835        break;
2836      }
2837    }
2838  }
2839  rowArray->setNumElements(0);
2840  columnArray->setNumElements(0);
2841}
2842// Restores bound to original bound
2843void 
2844ClpPdco::originalBound( int iSequence)
2845{
2846  if (getFakeBound(iSequence)!=noFake)
2847    numberFake_--;;
2848  if (iSequence>=numberColumns_) {
2849    // rows
2850    int iRow = iSequence-numberColumns_;
2851    rowLowerWork_[iRow]=rowLower_[iRow];
2852    rowUpperWork_[iRow]=rowUpper_[iRow];
2853    if (rowScale_) {
2854      if (rowLowerWork_[iRow]>-1.0e50)
2855        rowLowerWork_[iRow] *= rowScale_[iRow];
2856      if (rowUpperWork_[iRow]<1.0e50)
2857        rowUpperWork_[iRow] *= rowScale_[iRow];
2858    }
2859  } else {
2860    // columns
2861    columnLowerWork_[iSequence]=columnLower_[iSequence];
2862    columnUpperWork_[iSequence]=columnUpper_[iSequence];
2863    if (rowScale_) {
2864      double multiplier = 1.0/columnScale_[iSequence];
2865      if (columnLowerWork_[iSequence]>-1.0e50)
2866        columnLowerWork_[iSequence] *= multiplier;
2867      if (columnUpperWork_[iSequence]<1.0e50)
2868        columnUpperWork_[iSequence] *= multiplier;
2869    }
2870  }
2871  setFakeBound(iSequence,noFake);
2872}
2873/* As changeBounds but just changes new bounds for a single variable.
2874   Returns true if change */
2875bool 
2876ClpPdco::changeBound( int iSequence)
2877{
2878  // old values
2879  double oldLower=lower_[iSequence];
2880  double oldUpper=upper_[iSequence];
2881  double value=solution_[iSequence];
2882  bool modified=false;
2883  originalBound(iSequence);
2884  // original values
2885  double lowerValue=lower_[iSequence];
2886  double upperValue=upper_[iSequence];
2887  // back to altered values
2888  lower_[iSequence] = oldLower;
2889  upper_[iSequence] = oldUpper;
2890  if (getFakeBound(iSequence)!=noFake)
2891    numberFake_--;;
2892  if (value==oldLower) {
2893    if (upperValue > oldLower + dualBound_) {
2894      upper_[iSequence]=oldLower+dualBound_;
2895      setFakeBound(iSequence,upperFake);
2896      modified=true;
2897      numberFake_++;
2898    }
2899  } else if (value==oldUpper) {
2900    if (lowerValue < oldUpper - dualBound_) {
2901      lower_[iSequence]=oldUpper-dualBound_;
2902      setFakeBound(iSequence,lowerFake);
2903      modified=true;
2904      numberFake_++;
2905    }
2906  } else {
2907    assert(value==oldLower||value==oldUpper);
2908  }
2909  return modified;
2910}
2911// Perturbs problem
2912void 
2913ClpPdco::perturb()
2914{
2915  if (perturbation_>100)
2916    return; //perturbed already
2917  if (perturbation_==100)
2918    perturbation_=50; // treat as normal
2919  bool modifyRowCosts=false;
2920  // dual perturbation
2921  double perturbation=1.0e-20;
2922  // maximum fraction of cost to perturb
2923  double maximumFraction = 1.0e-5;
2924  double constantPerturbation = 100.0*dualTolerance_;
2925  const int * lengths = matrix_->getVectorLengths();
2926  int maxLength=0;
2927  int minLength=numberRows_;
2928  if (!numberIterations_&&perturbation_==50) {
2929    // See if we need to perturb
2930    double * sort = new double[numberColumns_];
2931    int i;
2932    for (i=0;i<numberColumns_;i++) {
2933      double value = fabs(objectiveWork_[i]);
2934      sort[i]=value;
2935    }
2936    std::sort(sort,sort+numberColumns_);
2937    int number=1;
2938    double last = sort[0];
2939    for (i=1;i<numberColumns_;i++) {
2940      if (last!=sort[i])
2941        number++;
2942      last=sort[i];
2943    }
2944    delete [] sort;
2945    if (number*4>numberColumns_)
2946      return; // good enough
2947  }
2948  int iColumn;
2949  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2950    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
2951      int length = lengths[iColumn];
2952      if (length>2) {
2953        maxLength = max(maxLength,length);
2954        minLength = min(minLength,length);
2955      }
2956    }
2957  }
2958  // If > 70 then do rows
2959  if (perturbation_>=70) {
2960    modifyRowCosts=true;
2961    perturbation_ -= 20;
2962    printf("Row costs modified, ");
2963  }
2964  bool uniformChange=false;
2965  if (perturbation_>50) {
2966    // Experiment
2967    // maximumFraction could be 1.0e-10 to 1.0
2968    double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
2969    maximumFraction = m[min(perturbation_-51,10)];
2970  }
2971  int iRow;
2972  double smallestNonZero=1.0e100;
2973  if (perturbation_>=50) {
2974    perturbation = 1.0e-8;
2975    bool allSame=true;
2976    double lastValue=0.0;
2977    for (iRow=0;iRow<numberRows_;iRow++) {
2978      double lo = rowLowerWork_[iRow];
2979      double up = rowUpperWork_[iRow];
2980      if (lo<up) {
2981        double value = fabs(rowObjectiveWork_[iRow]);
2982        perturbation = max(perturbation,value);
2983        if (value) {
2984          modifyRowCosts=true;
2985          smallestNonZero = min(smallestNonZero,value);
2986        }
2987      } 
2988      if (lo&&lo>-1.0e10) {
2989        lo=fabs(lo);
2990        if (!lastValue) 
2991          lastValue=lo;
2992        else if (fabs(lo-lastValue)>1.0e-7)
2993          allSame=false;
2994      }
2995      if (up&&up<1.0e10) {
2996        up=fabs(up);
2997        if (!lastValue) 
2998          lastValue=up;
2999        else if (fabs(up-lastValue)>1.0e-7)
3000          allSame=false;
3001      }
3002    }
3003    for (iColumn=0;iColumn<numberColumns_;iColumn++) { 
3004      double lo = columnLowerWork_[iColumn];
3005      double up = columnUpperWork_[iColumn];
3006      if (lo<up) {
3007        double value = 
3008          fabs(objectiveWork_[iColumn]);
3009        perturbation = max(perturbation,value);
3010        if (value) {
3011          smallestNonZero = min(smallestNonZero,value);
3012        }
3013      }
3014      if (lo&&lo>-1.0e10) {
3015        lo=fabs(lo);
3016        if (!lastValue) 
3017          lastValue=lo;
3018        else if (fabs(lo-lastValue)>1.0e-7)
3019          allSame=false;
3020      }
3021      if (up&&up<1.0e10) {
3022        up=fabs(up);
3023        if (!lastValue) 
3024          lastValue=up;
3025        else if (fabs(up-lastValue)>1.0e-7)
3026          allSame=false;
3027      }
3028    }
3029    if (allSame) {
3030      // Really hit perturbation
3031      maximumFraction = max(1.0e-2*lastValue,maximumFraction);
3032    }
3033    perturbation = min(perturbation,smallestNonZero/maximumFraction);
3034  } else {
3035    // user is in charge
3036    maximumFraction = 1.0e-1;
3037    // but some experiments
3038    if (perturbation_<=-900) {
3039      modifyRowCosts=true;
3040      perturbation_ += 1000;
3041      printf("Row costs modified, ");
3042    }
3043    if (perturbation_<=-10) {
3044      perturbation_ += 10; 
3045      maximumFraction = 1.0;
3046      if ((-perturbation_)%100>=10) {
3047        uniformChange=true;
3048        perturbation_+=20;
3049      }
3050      while (perturbation_<-10) {
3051        perturbation_ += 100;
3052        maximumFraction *= 1.0e-1;
3053      }
3054    }
3055    perturbation = pow(10.0,perturbation_);
3056  }
3057  double largestZero=0.0;
3058  double largest=0.0;
3059  double largestPerCent=0.0;
3060  // modify costs
3061  bool printOut=(handler_->logLevel()==63);
3062  printOut=false;
3063  if (modifyRowCosts) {
3064    for (iRow=0;iRow<numberRows_;iRow++) {
3065      if (rowLowerWork_[iRow]<rowUpperWork_[iRow]) {
3066        double value = perturbation;
3067        double currentValue = rowObjectiveWork_[iRow];
3068        value = min(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
3069        if (rowLowerWork_[iRow]>-largeValue_) {
3070          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow])) 
3071            value *= CoinDrand48();
3072          else
3073            value *= -CoinDrand48();
3074        } else if (rowUpperWork_[iRow]<largeValue_) {
3075          value *= -CoinDrand48();
3076        } else {
3077          value=0.0;
3078        }
3079        if (currentValue) {
3080          largest = max(largest,fabs(value));
3081          if (fabs(value)>fabs(currentValue)*largestPerCent)
3082            largestPerCent=fabs(value/currentValue);
3083        } else {
3084          largestZero=max(largestZero,fabs(value));
3085        }
3086        if (printOut)
3087          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
3088        rowObjectiveWork_[iRow] += value;
3089      }
3090    }
3091  }
3092  // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
3093  // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
3094  double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
3095  //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
3096  double extraWeight=10.0;
3097  // Scale back if wanted
3098  double weight2[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
3099  if (constantPerturbation<99.0*dualTolerance_) {
3100    perturbation *= 0.1;
3101    extraWeight=0.5;
3102    memcpy(weight,weight2,sizeof(weight2));
3103  }
3104  // adjust weights if all columns long
3105  double factor=1.0;
3106  if (maxLength) {
3107    factor = 3.0/(double) minLength;
3108  }
3109  // Make variables with more elements more expensive
3110  const double m1 = 0.5;
3111  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3112    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
3113      double value = perturbation;
3114      double currentValue = objectiveWork_[iColumn];
3115      value = min(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
3116      //value = min(value,constantPerturbation;+maximumFraction*fabs(currentValue));
3117      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
3118      if (uniformChange) {
3119        value = maximumFraction;
3120        value2=maximumFraction;
3121      }
3122      if (columnLowerWork_[iColumn]>-largeValue_) {
3123        if (fabs(columnLowerWork_[iColumn])<
3124            fabs(columnUpperWork_[iColumn])) {
3125          value *= (1.0-m1 +m1*CoinDrand48());
3126          value2 *= (1.0-m1 +m1*CoinDrand48());
3127        } else {
3128          value *= -(1.0-m1+m1*CoinDrand48());
3129          value2 *= -(1.0-m1+m1*CoinDrand48());
3130        }
3131      } else if (columnUpperWork_[iColumn]<largeValue_) {
3132        value *= -(1.0-m1+m1*CoinDrand48());
3133        value2 *= -(1.0-m1+m1*CoinDrand48());
3134      } else {
3135        value=0.0;
3136      }
3137      if (value) {
3138        int length = lengths[iColumn];
3139        if (length>3) {
3140          length = (int) ((double) length * factor);
3141          length = max(3,length);
3142        }
3143        double multiplier;
3144#if 1
3145        if (length<10)
3146          multiplier=weight[length];
3147        else
3148          multiplier = weight[10];
3149#else
3150        if (length<10)
3151          multiplier=weight[length];
3152        else
3153          multiplier = weight[10]+extraWeight*(length-10);
3154        multiplier *= 0.5;
3155#endif
3156        value *= multiplier;
3157        value = min (value,value2);
3158        if (fabs(value)<=dualTolerance_)
3159          value=0.0;
3160        if (currentValue) {
3161          largest = max(largest,fabs(value));
3162          if (fabs(value)>fabs(currentValue)*largestPerCent)
3163            largestPerCent=fabs(value/currentValue);
3164        } else {
3165          largestZero=max(largestZero,fabs(value));
3166        }
3167        if (printOut)
3168          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
3169        objectiveWork_[iColumn] += value;
3170      }
3171    }
3172  }
3173  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
3174    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
3175    <<CoinMessageEol;
3176  // and zero changes
3177  //int nTotal = numberRows_+numberColumns_;
3178  //memset(cost_+nTotal,0,nTotal*sizeof(double));
3179  // say perturbed
3180  perturbation_=101;
3181
3182}
3183/* For strong branching.  On input lower and upper are new bounds
3184   while on output they are change in objective function values
3185   (>1.0e50 infeasible).
3186   Return code is 0 if nothing interesting, -1 if infeasible both
3187   ways and +1 if infeasible one way (check values to see which one(s))
3188*/
3189int ClpPdco::strongBranching(int numberVariables,const int * variables,
3190                                    double * newLower, double * newUpper,
3191                                    double ** outputSolution,
3192                                    int * outputStatus, int * outputIterations,
3193                                    bool stopOnFirstInfeasible,
3194                                    bool alwaysFinish)
3195{
3196  int i;
3197  int returnCode=0;
3198  double saveObjectiveValue = objectiveValue_;
3199#if 1
3200  algorithm_ = -1;
3201
3202  //scaling(false);
3203
3204  // put in standard form (and make row copy)
3205  // create modifiable copies of model rim and do optional scaling
3206  createRim(7+8+16,true);
3207
3208  // change newLower and newUpper if scaled
3209
3210  // Do initial factorization
3211  // and set certain stuff
3212  // We can either set increasing rows so ...IsBasic gives pivot row
3213  // or we can just increment iBasic one by one
3214  // for now let ...iBasic give pivot row
3215  factorization_->increasingRows(2);
3216  // row activities have negative sign
3217  factorization_->slackValue(-1.0);
3218  factorization_->zeroTolerance(1.0e-13);
3219
3220  int factorizationStatus = internalFactorize(0);
3221  if (factorizationStatus<0)
3222    return 1; // some error
3223  else if (factorizationStatus)
3224    handler_->message(CLP_SINGULARITIES,messages_)
3225      <<factorizationStatus
3226      <<CoinMessageEol;
3227 
3228  // save stuff
3229  ClpFactorization saveFactorization(*factorization_);
3230  // save basis and solution
3231  double * saveSolution = new double[numberRows_+numberColumns_];
3232  memcpy(saveSolution,solution_,
3233         (numberRows_+numberColumns_)*sizeof(double));
3234  unsigned char * saveStatus = 
3235    new unsigned char [numberRows_+numberColumns_];
3236  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
3237  // save bounds as createRim makes clean copies
3238  double * saveLower = new double[numberRows_+numberColumns_];
3239  memcpy(saveLower,lower_,
3240         (numberRows_+numberColumns_)*sizeof(double));
3241  double * saveUpper = new double[numberRows_+numberColumns_];
3242  memcpy(saveUpper,upper_,
3243         (numberRows_+numberColumns_)*sizeof(double));
3244  double * saveObjective = new double[numberRows_+numberColumns_];
3245  memcpy(saveObjective,cost_,
3246         (numberRows_+numberColumns_)*sizeof(double));
3247  int * savePivot = new int [numberRows_];
3248  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
3249  // need to save/restore weights.
3250
3251  int iSolution = 0;
3252  for (i=0;i<numberVariables;i++) {
3253    int iColumn = variables[i];
3254    double objectiveChange;
3255    double saveBound;
3256   
3257    // try down
3258   
3259    saveBound = columnUpper_[iColumn];
3260    // external view - in case really getting optimal
3261    columnUpper_[iColumn] = newUpper[i];
3262    if (scalingFlag_<=0) 
3263      upper_[iColumn] = newUpper[i];
3264    else 
3265      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
3266    // Start of fast iterations
3267    int status = fastDual(alwaysFinish);
3268    if (status) {
3269      // not finished - might be optimal
3270      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
3271      double limit = 0.0;
3272      getDblParam(ClpDualObjectiveLimit, limit);
3273      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
3274          optimizationDirection_*limit) { 
3275        problemStatus_=0;
3276      } 
3277      status=problemStatus_;
3278    }
3279
3280    if (scalingFlag_<=0) {
3281      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
3282    } else {
3283      int j;
3284      double * sol = outputSolution[iSolution];
3285      for (j=0;j<numberColumns_;j++) 
3286        sol[j] = solution_[j]*columnScale_[j];
3287    }
3288    outputStatus[iSolution]=status;
3289    outputIterations[iSolution]=numberIterations_;
3290    iSolution++;
3291    // restore
3292    memcpy(solution_,saveSolution,
3293           (numberRows_+numberColumns_)*sizeof(double));
3294    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3295    memcpy(lower_,saveLower,
3296           (numberRows_+numberColumns_)*sizeof(double));
3297    memcpy(upper_,saveUpper,
3298           (numberRows_+numberColumns_)*sizeof(double));
3299    memcpy(cost_,saveObjective,
3300         (numberRows_+numberColumns_)*sizeof(double));
3301    columnUpper_[iColumn] = saveBound;
3302    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
3303    delete factorization_;
3304    factorization_ = new ClpFactorization(saveFactorization);
3305
3306    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
3307      objectiveChange = objectiveValue_-saveObjectiveValue;
3308    } else {
3309      objectiveChange = 1.0e100;
3310    }
3311    newUpper[i]=objectiveChange;
3312#ifdef CLP_DEBUG
3313    printf("down on %d costs %g\n",iColumn,objectiveChange);
3314#endif
3315         
3316    // try up
3317   
3318    saveBound = columnLower_[iColumn];
3319    // external view - in case really getting optimal
3320    columnLower_[iColumn] = newLower[i];
3321    if (scalingFlag_<=0) 
3322      lower_[iColumn] = newLower[i];
3323    else 
3324      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
3325    // Start of fast iterations
3326    status = fastDual(alwaysFinish);
3327    if (status) {
3328      // not finished - might be optimal
3329      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
3330      double limit = 0.0;
3331      getDblParam(ClpDualObjectiveLimit, limit);
3332      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
3333          optimizationDirection_*limit) { 
3334        problemStatus_=0;
3335      } 
3336      status=problemStatus_;
3337    }
3338    if (scalingFlag_<=0) {
3339      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
3340    } else {
3341      int j;
3342      double * sol = outputSolution[iSolution];
3343      for (j=0;j<numberColumns_;j++) 
3344        sol[j] = solution_[j]*columnScale_[j];
3345    }
3346    outputStatus[iSolution]=status;
3347    outputIterations[iSolution]=numberIterations_;
3348    iSolution++;
3349
3350    // restore
3351    memcpy(solution_,saveSolution,
3352           (numberRows_+numberColumns_)*sizeof(double));
3353    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3354    memcpy(lower_,saveLower,
3355           (numberRows_+numberColumns_)*sizeof(double));
3356    memcpy(upper_,saveUpper,
3357           (numberRows_+numberColumns_)*sizeof(double));
3358    memcpy(cost_,saveObjective,
3359         (numberRows_+numberColumns_)*sizeof(double));
3360    columnLower_[iColumn] = saveBound;
3361    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
3362    delete factorization_;
3363    factorization_ = new ClpFactorization(saveFactorization);
3364
3365    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
3366      objectiveChange = objectiveValue_-saveObjectiveValue;
3367    } else {
3368      objectiveChange = 1.0e100;
3369    }
3370    newLower[i]=objectiveChange;
3371#ifdef CLP_DEBUG
3372    printf("up on %d costs %g\n",iColumn,objectiveChange);
3373#endif
3374         
3375    /* Possibilities are:
3376       Both sides feasible - store
3377       Neither side feasible - set objective high and exit if desired
3378       One side feasible - change bounds and resolve
3379    */
3380    if (newUpper[i]<1.0e100) {
3381      if(newLower[i]<1.0e100) {
3382        // feasible - no action
3383      } else {
3384        // up feasible, down infeasible
3385        returnCode=1;
3386        if (stopOnFirstInfeasible)
3387          break;
3388      }
3389    } else {
3390      if(newLower[i]<1.0e100) {
3391        // down feasible, up infeasible
3392        returnCode=1;
3393        if (stopOnFirstInfeasible)
3394          break;
3395      } else {
3396        // neither side feasible
3397        returnCode=-1;
3398        break;
3399      }
3400    }
3401  }
3402  delete [] saveSolution;
3403  delete [] saveLower;
3404  delete [] saveUpper;
3405  delete [] saveObjective;
3406  delete [] saveStatus;
3407  delete [] savePivot;
3408
3409  // Get rid of some arrays and empty factorization
3410  deleteRim();
3411#else
3412  // save basis and solution
3413  double * rowSolution = new double[numberRows_];
3414  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
3415  double * columnSolution = new double[numberColumns_];
3416  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
3417  unsigned char * saveStatus = 
3418    new unsigned char [numberRows_+numberColumns_];
3419  if (!status_) {
3420    // odd but anyway setup all slack basis
3421    status_ = new unsigned char [numberColumns_+numberRows_];
3422    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3423  }
3424  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
3425  int iSolution =0;
3426  for (i=0;i<numberVariables;i++) {
3427    int iColumn = variables[i];
3428    double objectiveChange;
3429   
3430    // try down
3431   
3432    double saveUpper = columnUpper_[iColumn];
3433    columnUpper_[iColumn] = newUpper[i];
3434    int status=dual(0);
3435    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
3436    outputStatus[iSolution]=status;
3437    outputIterations[iSolution]=numberIterations_;
3438    iSolution++;
3439
3440    // restore
3441    columnUpper_[iColumn] = saveUpper;
3442    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
3443    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
3444    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3445
3446    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
3447      objectiveChange = objectiveValue_-saveObjectiveValue;
3448    } else {
3449      objectiveChange = 1.0e100;
3450    }
3451    newUpper[i]=objectiveChange;
3452#ifdef CLP_DEBUG
3453    printf("down on %d costs %g\n",iColumn,objectiveChange);
3454#endif
3455         
3456    // try up
3457   
3458    double saveLower = columnLower_[iColumn];
3459    columnLower_[iColumn] = newLower[i];
3460    status=dual(0);
3461    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
3462    outputStatus[iSolution]=status;
3463    outputIterations[iSolution]=numberIterations_;
3464    iSolution++;
3465
3466    // restore
3467    columnLower_[iColumn] = saveLower;
3468    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
3469    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
3470    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
3471
3472    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
3473      objectiveChange = objectiveValue_-saveObjectiveValue;
3474    } else {
3475      objectiveChange = 1.0e100;
3476    }
3477    newLower[i]=objectiveChange;
3478#ifdef CLP_DEBUG
3479    printf("up on %d costs %g\n",iColumn,objectiveChange);
3480#endif
3481         
3482    /* Possibilities are:
3483       Both sides feasible - store
3484       Neither side feasible - set objective high and exit
3485       One side feasible - change bounds and resolve
3486    */
3487    if (newUpper[i]<1.0e100) {
3488      if(newLower[i]<1.0e100) {
3489        // feasible - no action
3490      } else {
3491        // up feasible, down infeasible
3492        returnCode=1;
3493        if (stopOnFirstInfeasible)
3494          break;
3495      }
3496    } else {
3497      if(newLower[i]<1.0e100) {
3498        // down feasible, up infeasible
3499        returnCode=1;
3500        if (stopOnFirstInfeasible)
3501          break;
3502      } else {
3503        // neither side feasible
3504        returnCode=-1;
3505        break;
3506      }
3507    }
3508  }
3509  delete [] rowSolution;
3510  delete [] columnSolution;
3511  delete [] saveStatus;
3512#endif
3513  objectiveValue_ = saveObjectiveValue;
3514  return returnCode;
3515}
3516// treat no pivot as finished (unless interesting)
3517int ClpPdco::fastDual(bool alwaysFinish)
3518{
3519  algorithm_ = -1;
3520  // save data
3521  ClpDataSave data = saveData();
3522  dualTolerance_=dblParam_[ClpDualTolerance];
3523  primalTolerance_=dblParam_[ClpPrimalTolerance];
3524
3525  // save dual bound
3526  double saveDualBound = dualBound_;
3527
3528  double objectiveChange;
3529  // for dual we will change bounds using dualBound_
3530  // for this we need clean basis so it is after factorize
3531  gutsOfSolution(NULL,NULL);
3532  numberFake_ =0; // Number of variables at fake bounds
3533  changeBounds(true,NULL,objectiveChange);
3534
3535  problemStatus_ = -1;
3536  numberIterations_=0;
3537  factorization_->sparseThreshold(0);
3538  factorization_->goSparse();
3539
3540  int lastCleaned=0; // last time objective or bounds cleaned up
3541
3542  // number of times we have declared optimality
3543  numberTimesOptimal_=0;
3544
3545  // This says whether to restore things etc
3546  int factorType=0;
3547  /*
3548    Status of problem:
3549    0 - optimal
3550    1 - infeasible
3551    2 - unbounded
3552    -1 - iterating
3553    -2 - factorization wanted
3554    -3 - redo checking without factorization
3555    -4 - looks infeasible
3556
3557    BUT also from whileIterating return code is:
3558
3559   -1 iterations etc
3560   -2 inaccuracy
3561   -3 slight inaccuracy (and done iterations)
3562   +0 looks optimal (might be unbounded - but we will investigate)
3563   +1 looks infeasible
3564   +3 max iterations
3565
3566  */
3567
3568  int returnCode = 0;
3569
3570  int iRow,iColumn;
3571  while (problemStatus_<0) {
3572    // clear
3573    for (iRow=0;iRow<4;iRow++) {
3574      rowArray_[iRow]->clear();
3575    }   
3576   
3577    for (iColumn=0;iColumn<2;iColumn++) {
3578      columnArray_[iColumn]->clear();
3579    }   
3580
3581    // give matrix (and model costs and bounds a chance to be
3582    // refreshed (normally null)
3583    matrix_->refresh(this);
3584    // may factorize, checks if problem finished
3585    // should be able to speed this up on first time
3586    statusOfProblemInDual(lastCleaned,factorType,progress_,NULL);
3587
3588    // Say good factorization
3589    factorType=1;
3590
3591    // Do iterations
3592    if (problemStatus_<0) {
3593      double * givenPi=NULL;
3594      returnCode = whileIterating(givenPi);
3595      if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
3596        double limit = 0.0;
3597        getDblParam(ClpDualObjectiveLimit, limit);
3598        if(fabs(limit)>1.0e30||objectiveValue()*optimizationDirection_<
3599           optimizationDirection_*limit|| 
3600           numberAtFakeBound()) {
3601          returnCode=1;
3602          secondaryStatus_ = 1; // and say was on cutoff
3603          // can't say anything interesting - might as well return
3604#ifdef CLP_DEBUG
3605          printf("returning from fastDual after %d iterations with code %d\n",
3606                 numberIterations_,returnCode);
3607#endif
3608          break;
3609        }
3610      }
3611      returnCode=0;
3612    }
3613  }
3614
3615  // clear
3616  for (iRow=0;iRow<4;iRow++) {
3617    rowArray_[iRow]->clear();
3618  }   
3619 
3620  for (iColumn=0;iColumn<2;iColumn++) {
3621    columnArray_[iColumn]->clear();
3622  }   
3623  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
3624  // Restore any saved stuff
3625  restoreData(data);
3626  dualBound_ = saveDualBound;
3627  return returnCode;
3628}
3629/* Checks number of variables at fake bounds.  This is used by fastDual
3630   so can exit gracefully before end */
3631int 
3632ClpPdco::numberAtFakeBound()
3633{
3634  int iSequence;
3635  int numberFake=0;
3636 
3637  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3638    FakeBound bound = getFakeBound(iSequence);
3639    switch(getStatus(iSequence)) {
3640
3641    case basic:
3642      break;
3643    case isFree:
3644    case superBasic:
3645    case ClpSimplex::isFixed:
3646      assert (bound==noFake);
3647      break;
3648    case atUpperBound:
3649      if (bound==upperFake||bound==bothFake)
3650        numberFake++;
3651      break;
3652    case atLowerBound:
3653      if (bound==lowerFake||bound==bothFake)
3654        numberFake++;
3655      break;
3656    }
3657  }
3658  numberFake_ = numberFake;
3659  return numberFake;
3660}
3661/* Pivot out a variable and choose an incoing one.  Assumes dual
3662   feasible - will not go through a reduced cost. 
3663   Returns step length in theta
3664   Returns ray in ray_ (or NULL if no pivot)
3665   Return codes as before but -1 means no acceptable pivot
3666*/
3667int 
3668ClpPdco::pivotResult()
3669{
3670  abort();
3671  return 0;
3672}
3673/*
3674   Row array has row part of pivot row
3675   Column array has column part.
3676   This is used in dual values pass
3677*/
3678int
3679ClpPdco::checkPossibleValuesMove(CoinIndexedVector * rowArray,
3680                                        CoinIndexedVector * columnArray,
3681                                        double acceptablePivot,
3682                                        CoinBigIndex * dubiousWeights)
3683{
3684  double * work;
3685  int number;
3686  int * which;
3687  int iSection;
3688
3689  double tolerance = dualTolerance_*1.001;
3690
3691  double thetaDown = 1.0e31;
3692  double changeDown ;
3693  double thetaUp = 1.0e31;
3694  double bestAlphaDown = acceptablePivot*0.99999;
3695  double bestAlphaUp = acceptablePivot*0.99999;
3696  int sequenceDown =-1;
3697  int sequenceUp = sequenceOut_;
3698
3699  double djBasic = dj_[sequenceOut_];
3700  if (djBasic>0.0) {
3701    // basic at lower bound so directionOut_ 1 and -1 in pivot row
3702    // dj will go to zero on other way
3703    thetaUp = djBasic;
3704    changeDown = -lower_[sequenceOut_];
3705  } else {
3706    // basic at upper bound so directionOut_ -1 and 1 in pivot row
3707    // dj will go to zero on other way
3708    thetaUp = -djBasic;
3709    changeDown = upper_[sequenceOut_];
3710  }
3711  bestAlphaUp = 1.0;
3712  int addSequence;
3713
3714  double alphaUp=0.0;
3715  double alphaDown=0.0;
3716
3717  for (iSection=0;iSection<2;iSection++) {
3718
3719    int i;
3720    if (!iSection) {
3721      work = rowArray->denseVector();
3722      number = rowArray->getNumElements();
3723      which = rowArray->getIndices();
3724      addSequence = numberColumns_;
3725    } else {
3726      work = columnArray->denseVector();
3727      number = columnArray->getNumElements();
3728      which = columnArray->getIndices();
3729      addSequence = 0;
3730    }
3731   
3732    for (i=0;i<number;i++) {
3733      int iSequence = which[i];
3734      int iSequence2 = iSequence + addSequence;
3735      double alpha;
3736      double oldValue;
3737      double value;
3738
3739      switch(getStatus(iSequence2)) {
3740         
3741      case basic:
3742        break;
3743      case ClpSimplex::isFixed:
3744        alpha = work[i];
3745        changeDown += alpha*upper_[iSequence2];
3746        break;
3747      case isFree:
3748      case superBasic:
3749        alpha = work[i];
3750        // dj must be effectively zero as dual feasible
3751        if (fabs(alpha)>bestAlphaUp) {
3752          thetaDown = 0.0;
3753          thetaUp = 0.0;
3754          bestAlphaDown = fabs(alpha);
3755          bestAlphaUp = bestAlphaUp;
3756          sequenceDown =iSequence2;
3757          sequenceUp = sequenceDown;
3758          alphaUp = alpha;
3759          alphaDown = alpha;
3760        }
3761        break;
3762      case atUpperBound:
3763        alpha = work[i];
3764        oldValue = dj_[iSequence2];
3765        changeDown += alpha*upper_[iSequence2];
3766        if (alpha>=acceptablePivot) {
3767          // might do other way
3768          value = oldValue+thetaUp*alpha;
3769          if (value>-tolerance) {
3770            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
3771              thetaUp = -oldValue/alpha;
3772              bestAlphaUp = fabs(alpha);
3773              sequenceUp = iSequence2;
3774              alphaUp = alpha;
3775            }
3776          }
3777        } else if (alpha<=-acceptablePivot) {
3778          // might do this way
3779          value = oldValue-thetaDown*alpha;
3780          if (value>-tolerance) {
3781            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
3782              thetaDown = oldValue/alpha;
3783              bestAlphaDown = fabs(alpha);
3784              sequenceDown = iSequence2;
3785              alphaDown = alpha;
3786            }
3787          }
3788        }
3789        break;
3790      case atLowerBound:
3791        alpha = work[i];
3792        oldValue = dj_[iSequence2];
3793        changeDown += alpha*lower_[iSequence2];
3794        if (alpha<=-acceptablePivot) {
3795          // might do other way
3796          value = oldValue+thetaUp*alpha;
3797          if (value<tolerance) {
3798            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
3799              thetaUp = -oldValue/alpha;
3800              bestAlphaUp = fabs(alpha);
3801              sequenceUp = iSequence2;
3802              alphaUp = alpha;
3803            }
3804          }
3805        } else if (alpha>=acceptablePivot) {
3806          // might do this way
3807          value = oldValue-thetaDown*alpha;
3808          if (value<tolerance) {
3809            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
3810              thetaDown = oldValue/alpha;
3811              bestAlphaDown = fabs(alpha);
3812              sequenceDown = iSequence2;
3813              alphaDown = alpha;
3814            }
3815          }
3816        }
3817        break;
3818      }
3819    }
3820  }
3821  int returnCode;
3822  thetaUp *= -1.0;
3823  double changeUp = -thetaUp * changeDown;
3824  changeDown = -thetaDown*changeDown;
3825  if (max(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
3826    // largest
3827    if (fabs(alphaDown)<fabs(alphaUp)) {
3828      sequenceDown =-1;
3829    }
3830  }
3831  // choose
3832  if (changeDown>changeUp&&sequenceDown>=0) {
3833    theta_ = thetaDown;
3834    sequenceIn_ = sequenceDown;
3835    alpha_ = alphaDown;
3836#ifdef CLP_DEBUG
3837    if ((handler_->logLevel()&32))
3838      printf("predicted way - dirout %d, change %g,%g theta %g\n",
3839             directionOut_,changeDown,changeUp,theta_);
3840#endif
3841    returnCode = 0;
3842  } else {
3843    theta_ = thetaUp;
3844    sequenceIn_ = sequenceUp;
3845    alpha_ = alphaUp;
3846    if (sequenceIn_!=sequenceOut_) {
3847#ifdef CLP_DEBUG
3848      if ((handler_->logLevel()&32))
3849        printf("opposite way - dirout %d, change %g,%g theta %g\n",
3850               directionOut_,changeDown,changeUp,theta_);
3851#endif
3852      returnCode = 0;
3853    } else {
3854#ifdef CLP_DEBUG
3855      if ((handler_->logLevel()&32))
3856        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
3857               directionOut_,changeDown,changeUp,theta_);
3858#endif
3859      returnCode = 1;
3860    }
3861  }
3862  return returnCode;
3863}
3864/*
3865   This sees if we can move duals in dual values pass.
3866   This is done before any pivoting
3867*/
3868void ClpPdco::doEasyOnesInValuesPass(double * dj)
3869{
3870  // Get column copy
3871  CoinPackedMatrix * columnCopy = matrix();
3872  // Get a row copy in standard format
3873  CoinPackedMatrix copy;
3874  copy.reverseOrderedCopyOf(*columnCopy);
3875  // get matrix data pointers
3876  const int * column = copy.getIndices();
3877  const CoinBigIndex * rowStart = copy.getVectorStarts();
3878  const int * rowLength = copy.getVectorLengths(); 
3879  const double * elementByRow = copy.getElements();
3880  double tolerance = dualTolerance_*1.001;
3881
3882  int iRow;
3883#ifdef CLP_DEBUG
3884  {
3885    double value5=0.0;
3886    int i;
3887    for (i=0;i<numberRows_+numberColumns_;i++) {
3888      if (dj[i]<-1.0e-6)
3889        value5 += dj[i]*upper_[i];
3890      else if (dj[i] >1.0e-6)
3891        value5 += dj[i]*lower_[i];
3892    }
3893    printf("Values objective Value before %g\n",value5);
3894  }
3895#endif
3896  // for scaled row
3897  double * scaled=NULL;
3898  if (rowScale_)
3899    scaled = new double[numberColumns_];
3900  for (iRow=0;iRow<numberRows_;iRow++) {
3901
3902    int iSequence = iRow + numberColumns_;
3903    double djBasic = dj[iSequence];
3904    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
3905
3906      double changeUp ;
3907      // always -1 in pivot row
3908      if (djBasic>0.0) {
3909        // basic at lower bound
3910        changeUp = -lower_[iSequence];
3911      } else {
3912        // basic at upper bound
3913        changeUp = upper_[iSequence];
3914      }
3915      bool canMove=true;
3916      int i;
3917      const double * thisElements = elementByRow + rowStart[iRow]; 
3918      const int * thisIndices = column+rowStart[iRow];
3919      if (rowScale_) {
3920        // scale row
3921        double scale = rowScale_[iRow];
3922        for (i = 0;i<rowLength[iRow];i++) {
3923          int iColumn = thisIndices[i];
3924          double alpha = thisElements[i];
3925          scaled[i] = scale*alpha*columnScale_[iColumn];
3926        }
3927        thisElements = scaled;
3928      }
3929      for (i = 0;i<rowLength[iRow];i++) {
3930        int iColumn = thisIndices[i];
3931        double alpha = thisElements[i];
3932        double oldValue = dj[iColumn];;
3933        double value;
3934
3935        switch(getStatus(iColumn)) {
3936         
3937        case basic:
3938          if (dj[iColumn]<-tolerance&&
3939              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
3940            // at ub
3941            changeUp += alpha*upper_[iColumn];
3942            // might do other way
3943            value = oldValue+djBasic*alpha;
3944            if (value>tolerance) 
3945              canMove=false;
3946          } else if (dj[iColumn]>tolerance&&
3947              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
3948            changeUp += alpha*lower_[iColumn];
3949            // might do other way
3950            value = oldValue+djBasic*alpha;
3951            if (value<-tolerance)
3952              canMove=false;
3953          } else {
3954            canMove=false;
3955          }
3956          break;
3957        case ClpSimplex::isFixed:
3958          changeUp += alpha*upper_[iColumn];
3959          break;
3960        case isFree:
3961        case superBasic:
3962          canMove=false;
3963        break;
3964      case atUpperBound:
3965        changeUp += alpha*upper_[iColumn];
3966        // might do other way
3967        value = oldValue+djBasic*alpha;
3968        if (value>tolerance) 
3969          canMove=false;
3970        break;
3971      case atLowerBound:
3972        changeUp += alpha*lower_[iColumn];
3973        // might do other way
3974        value = oldValue+djBasic*alpha;
3975        if (value<-tolerance)
3976          canMove=false;
3977        break;
3978        }
3979      }
3980      if (canMove) {
3981        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
3982          // move
3983          for (i = 0;i<rowLength[iRow];i++) {
3984            int iColumn = thisIndices[i];
3985            double alpha = thisElements[i];
3986            dj[iColumn] += djBasic * alpha;
3987          }
3988          dj[iSequence]=0.0;
3989#ifdef CLP_DEBUG
3990          {
3991            double value5=0.0;
3992            int i;
3993            for (i=0;i<numberRows_+numberColumns_;i++) {
3994              if (dj[i]<-1.0e-6)
3995                value5 += dj[i]*upper_[i];
3996              else if (dj[i] >1.0e-6)
3997                value5 += dj[i]*lower_[i];
3998            }
3999            printf("Values objective Value after row %d old dj %g %g\n",
4000                   iRow,djBasic,value5);
4001          }
4002#endif
4003        }
4004      }
4005    }     
4006  }
4007  delete [] scaled;
4008}
4009int
4010ClpPdco::nextSuperBasic()
4011{
4012  if (firstFree_>=0) {
4013    int returnValue=firstFree_;
4014    int iColumn=firstFree_+1;
4015    for (;iColumn<numberRows_+numberColumns_;iColumn++) {
4016      if (getStatus(iColumn)==isFree) 
4017        if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
4018          break;
4019    }
4020    firstFree_ = iColumn;
4021    if (firstFree_==numberRows_+numberColumns_)
4022      firstFree_=-1;
4023    return returnValue;
4024  } else {
4025    return -1;
4026  }
4027}
Note: See TracBrowser for help on using the repository browser.