source: trunk/ClpSimplex.cpp @ 437

Last change on this file since 437 was 432, checked in by forrest, 16 years ago

change to allSlackBasis

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 190.0 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpPackedMatrix.hpp"
15#include "CoinIndexedVector.hpp"
16#include "ClpDualRowDantzig.hpp"
17#include "ClpDualRowSteepest.hpp"
18#include "ClpPrimalColumnDantzig.hpp"
19#include "ClpPrimalColumnSteepest.hpp"
20#include "ClpNonLinearCost.hpp"
21#include "ClpMessage.hpp"
22#include "ClpEventHandler.hpp"
23#include "ClpLinearObjective.hpp"
24#include "ClpHelperFunctions.hpp"
25#include <cfloat>
26
27#include <string>
28#include <stdio.h>
29#include <iostream>
30//#############################################################################
31
32ClpSimplex::ClpSimplex () :
33
34  ClpModel(),
35  columnPrimalInfeasibility_(0.0),
36  rowPrimalInfeasibility_(0.0),
37  columnPrimalSequence_(-2),
38  rowPrimalSequence_(-2), 
39  columnDualInfeasibility_(0.0),
40  rowDualInfeasibility_(0.0),
41  columnDualSequence_(-2),
42  rowDualSequence_(-2),
43  primalToleranceToGetOptimal_(-1.0),
44  remainingDualInfeasibility_(0.0),
45  largeValue_(1.0e15),
46  largestPrimalError_(0.0),
47  largestDualError_(0.0),
48  largestSolutionError_(0.0),
49  dualBound_(1.0e10),
50  alpha_(0.0),
51  theta_(0.0),
52  lowerIn_(0.0),
53  valueIn_(0.0),
54  upperIn_(0.0),
55  dualIn_(0.0),
56  lowerOut_(-1),
57  valueOut_(-1),
58  upperOut_(-1),
59  dualOut_(-1),
60  dualTolerance_(0.0),
61  primalTolerance_(0.0),
62  sumDualInfeasibilities_(0.0),
63  sumPrimalInfeasibilities_(0.0),
64  infeasibilityCost_(1.0e10),
65  sumOfRelaxedDualInfeasibilities_(0.0),
66  sumOfRelaxedPrimalInfeasibilities_(0.0),
67  lower_(NULL),
68  rowLowerWork_(NULL),
69  columnLowerWork_(NULL),
70  upper_(NULL),
71  rowUpperWork_(NULL),
72  columnUpperWork_(NULL),
73  cost_(NULL),
74  rowObjectiveWork_(NULL),
75  objectiveWork_(NULL),
76  sequenceIn_(-1),
77  directionIn_(-1),
78  sequenceOut_(-1),
79  directionOut_(-1),
80  pivotRow_(-1),
81  lastGoodIteration_(-100),
82  dj_(NULL),
83  rowReducedCost_(NULL),
84  reducedCostWork_(NULL),
85  solution_(NULL),
86  rowActivityWork_(NULL),
87  columnActivityWork_(NULL),
88  numberDualInfeasibilities_(0),
89  numberDualInfeasibilitiesWithoutFree_(0),
90  numberPrimalInfeasibilities_(100),
91  numberRefinements_(0),
92  pivotVariable_(NULL),
93  factorization_(NULL),
94  savedSolution_(NULL),
95  numberTimesOptimal_(0),
96  changeMade_(1),
97  algorithm_(0),
98  forceFactorization_(-1),
99  perturbation_(100),
100  nonLinearCost_(NULL),
101  specialOptions_(0),
102  lastBadIteration_(-999999),
103  lastFlaggedIteration_(-999999),
104  numberFake_(0),
105  progressFlag_(0),
106  firstFree_(-1),
107  numberExtraRows_(0),
108  maximumBasic_(0),
109  incomingInfeasibility_(1.0),
110  allowedInfeasibility_(10.0),
111  automaticScale_(0),
112  progress_(NULL)
113{
114  int i;
115  for (i=0;i<6;i++) {
116    rowArray_[i]=NULL;
117    columnArray_[i]=NULL;
118  }
119  saveStatus_=NULL;
120  // get an empty factorization so we can set tolerances etc
121  factorization_ = new ClpFactorization();
122  // Say sparse
123  factorization_->sparseThreshold(1);
124  // say Steepest pricing
125  dualRowPivot_ = new ClpDualRowSteepest();
126  // say Steepest pricing
127  primalColumnPivot_ = new ClpPrimalColumnSteepest();
128  solveType_=1; // say simplex based life form
129 
130}
131
132// Subproblem constructor
133ClpSimplex::ClpSimplex ( const ClpModel * rhs,
134                     int numberRows, const int * whichRow,
135                     int numberColumns, const int * whichColumn,
136                     bool dropNames, bool dropIntegers)
137  : ClpModel(rhs, numberRows, whichRow,
138             numberColumns,whichColumn,dropNames,dropIntegers),
139  columnPrimalInfeasibility_(0.0),
140  rowPrimalInfeasibility_(0.0),
141  columnPrimalSequence_(-2),
142  rowPrimalSequence_(-2), 
143  columnDualInfeasibility_(0.0),
144  rowDualInfeasibility_(0.0),
145  columnDualSequence_(-2),
146  rowDualSequence_(-2),
147  primalToleranceToGetOptimal_(-1.0),
148  remainingDualInfeasibility_(0.0),
149  largeValue_(1.0e15),
150  largestPrimalError_(0.0),
151  largestDualError_(0.0),
152  largestSolutionError_(0.0),
153  dualBound_(1.0e10),
154  alpha_(0.0),
155  theta_(0.0),
156  lowerIn_(0.0),
157  valueIn_(0.0),
158  upperIn_(0.0),
159  dualIn_(0.0),
160  lowerOut_(-1),
161  valueOut_(-1),
162  upperOut_(-1),
163  dualOut_(-1),
164  dualTolerance_(0.0),
165  primalTolerance_(0.0),
166  sumDualInfeasibilities_(0.0),
167  sumPrimalInfeasibilities_(0.0),
168  infeasibilityCost_(1.0e10),
169  sumOfRelaxedDualInfeasibilities_(0.0),
170  sumOfRelaxedPrimalInfeasibilities_(0.0),
171  lower_(NULL),
172  rowLowerWork_(NULL),
173  columnLowerWork_(NULL),
174  upper_(NULL),
175  rowUpperWork_(NULL),
176  columnUpperWork_(NULL),
177  cost_(NULL),
178  rowObjectiveWork_(NULL),
179  objectiveWork_(NULL),
180  sequenceIn_(-1),
181  directionIn_(-1),
182  sequenceOut_(-1),
183  directionOut_(-1),
184  pivotRow_(-1),
185  lastGoodIteration_(-100),
186  dj_(NULL),
187  rowReducedCost_(NULL),
188  reducedCostWork_(NULL),
189  solution_(NULL),
190  rowActivityWork_(NULL),
191  columnActivityWork_(NULL),
192  numberDualInfeasibilities_(0),
193  numberDualInfeasibilitiesWithoutFree_(0),
194  numberPrimalInfeasibilities_(100),
195  numberRefinements_(0),
196  pivotVariable_(NULL),
197  factorization_(NULL),
198  savedSolution_(NULL),
199  numberTimesOptimal_(0),
200  changeMade_(1),
201  algorithm_(0),
202  forceFactorization_(-1),
203  perturbation_(100),
204  nonLinearCost_(NULL),
205  specialOptions_(0),
206  lastBadIteration_(-999999),
207  lastFlaggedIteration_(-999999),
208  numberFake_(0),
209  progressFlag_(0),
210  firstFree_(-1),
211  numberExtraRows_(0),
212  maximumBasic_(0),
213  incomingInfeasibility_(1.0),
214  allowedInfeasibility_(10.0),
215  automaticScale_(0),
216  progress_(NULL)
217{
218  int i;
219  for (i=0;i<6;i++) {
220    rowArray_[i]=NULL;
221    columnArray_[i]=NULL;
222  }
223  saveStatus_=NULL;
224  // get an empty factorization so we can set tolerances etc
225  factorization_ = new ClpFactorization();
226  // say Steepest pricing
227  dualRowPivot_ = new ClpDualRowSteepest();
228  // say Steepest pricing
229  primalColumnPivot_ = new ClpPrimalColumnSteepest();
230  solveType_=1; // say simplex based life form
231 
232}
233
234//-----------------------------------------------------------------------------
235
236ClpSimplex::~ClpSimplex ()
237{
238  gutsOfDelete(0);
239  delete nonLinearCost_;
240}
241//#############################################################################
242void ClpSimplex::setLargeValue( double value) 
243{
244  if (value>0.0&&value<COIN_DBL_MAX)
245    largeValue_=value;
246}
247int
248ClpSimplex::gutsOfSolution ( double * givenDuals,
249                             const double * givenPrimals,
250                             bool valuesPass)
251{
252
253
254  // if values pass, save values of basic variables
255  double * save = NULL;
256  double oldValue=0.0;
257  if (valuesPass) {
258    assert(algorithm_>0); // only primal at present
259    assert(nonLinearCost_);
260    int iRow;
261    checkPrimalSolution( rowActivityWork_, columnActivityWork_);
262    // get correct bounds on all variables
263    nonLinearCost_->checkInfeasibilities(primalTolerance_);
264    oldValue = nonLinearCost_->largestInfeasibility();
265    save = new double[numberRows_];
266    for (iRow=0;iRow<numberRows_;iRow++) {
267      int iPivot=pivotVariable_[iRow];
268      save[iRow] = solution_[iPivot];
269    }
270  }
271  // do work
272  computePrimals(rowActivityWork_, columnActivityWork_);
273  // If necessary - override results
274  if (givenPrimals) {
275    memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double));
276    memset(rowActivityWork_,0,numberRows_*sizeof(double));
277    times(-1.0,columnActivityWork_,rowActivityWork_);
278  }
279  double objectiveModification = 0.0;
280  if (algorithm_>0&&nonLinearCost_!=NULL) {
281    // primal algorithm
282    // get correct bounds on all variables
283    // If  4 bit set - Force outgoing variables to exact bound (primal)
284    if ((specialOptions_&4)==0)
285      nonLinearCost_->checkInfeasibilities(primalTolerance_);
286    else
287      nonLinearCost_->checkInfeasibilities(0.0);
288    objectiveModification += nonLinearCost_->changeInCost();
289    if (nonLinearCost_->numberInfeasibilities())
290      handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
291        <<nonLinearCost_->changeInCost()
292        <<nonLinearCost_->numberInfeasibilities()
293        <<CoinMessageEol;
294  }
295  if (valuesPass) {
296#ifdef CLP_DEBUG
297    std::cout<<"Largest given infeasibility "<<oldValue
298             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
299#endif
300    int numberOut=0;
301    if (oldValue<incomingInfeasibility_
302        &&nonLinearCost_->largestInfeasibility()>
303        CoinMax(incomingInfeasibility_,allowedInfeasibility_)||
304        largestPrimalError_>1.0e-3) {
305      //printf("Original largest infeas %g, now %g, primalError %g\n",
306      //     oldValue,nonLinearCost_->largestInfeasibility(),
307      //     largestPrimalError_);
308      // throw out up to 1000 structurals
309      int iRow;
310      int * sort = new int[numberRows_];
311      // first put back solution and store difference
312      for (iRow=0;iRow<numberRows_;iRow++) {
313        int iPivot=pivotVariable_[iRow];
314        double difference = fabs(solution_[iPivot]-save[iRow]);
315        solution_[iPivot]=save[iRow];
316        save[iRow]=difference;
317      }
318      for (iRow=0;iRow<numberRows_;iRow++) {
319        int iPivot=pivotVariable_[iRow];
320
321        if (iPivot<numberColumns_) {
322          // column
323          double difference= save[iRow];
324          if (difference>1.0e-4) {
325            sort[numberOut]=iPivot;
326            save[numberOut++]=difference;
327          }
328        }
329      }
330      CoinSort_2(save, save + numberOut, sort,
331                 CoinFirstGreater_2<double, int>());
332      numberOut = CoinMin(1000,numberOut);
333      for (iRow=0;iRow<numberOut;iRow++) {
334        int iColumn=sort[iRow];
335        setColumnStatus(iColumn,superBasic);
336        if (fabs(solution_[iColumn])>1.0e10) {
337          if (upper_[iColumn]<0.0) {
338            solution_[iColumn]=upper_[iColumn];
339          } else if (lower_[iColumn]>0.0) {
340            solution_[iColumn]=lower_[iColumn];
341          } else {
342            solution_[iColumn]=0.0;
343          }
344        }
345      }
346      delete [] sort;
347    }
348    delete [] save;
349    if (numberOut)
350      return numberOut;
351  }
352
353  computeDuals(givenDuals);
354
355  // now check solutions
356  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
357  objectiveValue_ += objectiveModification/(objectiveScale_*rhsScale_);
358  checkDualSolution();
359  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
360                               largestDualError_>1.0e-2)) 
361    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
362      <<largestPrimalError_
363      <<largestDualError_
364      <<CoinMessageEol;
365  // Switch off false values pass indicator
366  if (!valuesPass&&algorithm_>0)
367    firstFree_ = -1;
368  return 0;
369}
370void
371ClpSimplex::computePrimals ( const double * rowActivities,
372                                     const double * columnActivities)
373{
374
375  //work space
376  CoinIndexedVector  * workSpace = rowArray_[0];
377
378  CoinIndexedVector arrayVector;
379  arrayVector.reserve(numberRows_+1);
380  CoinIndexedVector previousVector;
381  previousVector.reserve(numberRows_+1);
382
383  // accumulate non basic stuff
384
385  int iRow;
386  // order is this way for scaling
387  if (columnActivities!=columnActivityWork_)
388    ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
389  if (rowActivities!=rowActivityWork_)
390    ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
391  double * array = arrayVector.denseVector();
392  int * index = arrayVector.getIndices();
393  int number=0;
394  const double * rhsOffset = matrix_->rhsOffset(this,false,true);
395  if (!rhsOffset) {
396    // Use whole matrix every time to make it easier for ClpMatrixBase
397    // So zero out basic
398    for (iRow=0;iRow<numberRows_;iRow++) {
399      int iPivot=pivotVariable_[iRow];
400      solution_[iPivot] = 0.0;
401    }
402    // Extended solution before "update"
403    matrix_->primalExpanded(this,0);
404    times(-1.0,columnActivityWork_,array);
405    for (iRow=0;iRow<numberRows_;iRow++) {
406      double value = array[iRow] + rowActivityWork_[iRow];
407      if (value) {
408        array[iRow]=value;
409        index[number++]=iRow;
410      } else {
411        array[iRow]=0.0;
412      }
413    }
414  } else {
415    // we have an effective rhs lying around
416    // zero out basic (really just for slacks)
417    for (iRow=0;iRow<numberRows_;iRow++) {
418      int iPivot=pivotVariable_[iRow];
419      solution_[iPivot] = 0.0;
420    }
421    for (iRow=0;iRow<numberRows_;iRow++) {
422      double value = rhsOffset[iRow] + rowActivityWork_[iRow];
423      if (value) {
424        array[iRow]=value;
425        index[number++]=iRow;
426      } else {
427        array[iRow]=0.0;
428      }
429    }
430  }
431  arrayVector.setNumElements(number);
432#ifdef CLP_DEBUG
433  if (numberIterations_==-3840) {
434    int i;
435    for (i=0;i<numberRows_+numberColumns_;i++)
436      printf("%d status %d\n",i,status_[i]);
437    printf("xxxxx1\n");
438    for (i=0;i<numberRows_;i++)
439      if (array[i])
440        printf("%d rhs %g\n",i,array[i]);
441    printf("xxxxx2\n");
442    for (i=0;i<numberRows_+numberColumns_;i++)
443      if (getStatus(i)!=basic)
444        printf("%d non basic %g %g %g\n",i,lower_[i],solution_[i],upper_[i]);
445    printf("xxxxx3\n");
446  }
447#endif
448  // Ftran adjusted RHS and iterate to improve accuracy
449  double lastError=COIN_DBL_MAX;
450  int iRefine;
451  double * work = workSpace->denseVector();
452  CoinIndexedVector * thisVector = &arrayVector;
453  CoinIndexedVector * lastVector = &previousVector;
454  factorization_->updateColumn(workSpace,thisVector);
455#ifdef CLP_DEBUG
456  if (numberIterations_==-3840) {
457    int i;
458    for (i=0;i<numberRows_;i++)
459      if (array[i])
460        printf("%d after rhs %g\n",i,array[i]);
461    printf("xxxxx4\n");
462  }
463#endif
464  bool goodSolution=true;
465  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
466
467    int numberIn = thisVector->getNumElements();
468    int * indexIn = thisVector->getIndices();
469    double * arrayIn = thisVector->denseVector();
470    // put solution in correct place
471    if (!rhsOffset) {
472      int j;
473      for (j=0;j<numberIn;j++) {
474        iRow = indexIn[j];
475        int iPivot=pivotVariable_[iRow];
476        solution_[iPivot] = arrayIn[iRow];
477      }
478    } else {
479      for (iRow=0;iRow<numberRows_;iRow++) {
480        int iPivot=pivotVariable_[iRow];
481        solution_[iPivot] = arrayIn[iRow];
482      }
483    }
484    // Extended solution after "update"
485    matrix_->primalExpanded(this,1);
486    // check Ax == b  (for all)
487    // signal column generated matrix to just do basic (and gub)
488    unsigned int saveOptions = specialOptions();
489    setSpecialOptions(16);
490    times(-1.0,columnActivityWork_,work);
491    setSpecialOptions(saveOptions);
492    largestPrimalError_=0.0;
493    double multiplier = 131072.0;
494    for (iRow=0;iRow<numberRows_;iRow++) {
495      double value = work[iRow] + rowActivityWork_[iRow];
496      work[iRow] = value*multiplier;
497      if (fabs(value)>largestPrimalError_) {
498        largestPrimalError_=fabs(value);
499      }
500    }
501    if (largestPrimalError_>=lastError) {
502      // restore
503      CoinIndexedVector * temp = thisVector;
504      thisVector = lastVector;
505      lastVector=temp;
506      goodSolution=false;
507      break;
508    }
509    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-10) {
510      // try and make better
511      // save this
512      CoinIndexedVector * temp = thisVector;
513      thisVector = lastVector;
514      lastVector=temp;
515      int * indexOut = thisVector->getIndices();
516      int number=0;
517      array = thisVector->denseVector();
518      thisVector->clear();
519      for (iRow=0;iRow<numberRows_;iRow++) {
520        double value = work[iRow];
521        if (value) {
522          array[iRow]=value;
523          indexOut[number++]=iRow;
524          work[iRow]=0.0;
525        }
526      }
527      thisVector->setNumElements(number);
528      lastError=largestPrimalError_;
529      factorization_->updateColumn(workSpace,thisVector);
530      multiplier = 1.0/multiplier;
531      double * previous = lastVector->denseVector();
532      number=0;
533      for (iRow=0;iRow<numberRows_;iRow++) {
534        double value = previous[iRow] + multiplier*array[iRow];
535        if (value) {
536          array[iRow]=value;
537          indexOut[number++]=iRow;
538        } else {
539          array[iRow]=0.0;
540        }
541      }
542      thisVector->setNumElements(number);
543    } else {
544      break;
545    }
546  }
547
548  // solution as accurate as we are going to get
549  ClpFillN(work,numberRows_,0.0);
550  if (!goodSolution) {
551    array = thisVector->denseVector();
552    // put solution in correct place
553    for (iRow=0;iRow<numberRows_;iRow++) {
554      int iPivot=pivotVariable_[iRow];
555      solution_[iPivot] = array[iRow];
556    }
557  }
558#ifdef CLP_DEBUG
559  if (numberIterations_==-3840) {
560    exit(77);
561  }
562#endif
563}
564// now dual side
565void
566ClpSimplex::computeDuals(double * givenDjs)
567{
568  if (objective_->type()==1||!objective_->activated()) {
569    // Linear
570    //work space
571    CoinIndexedVector  * workSpace = rowArray_[0];
572   
573    CoinIndexedVector arrayVector;
574    arrayVector.reserve(numberRows_+1);
575    CoinIndexedVector previousVector;
576    previousVector.reserve(numberRows_+1);
577   
578   
579    int iRow;
580#ifdef CLP_DEBUG
581    workSpace->checkClear();
582#endif
583    double * array = arrayVector.denseVector();
584    int * index = arrayVector.getIndices();
585    int number=0;
586    if (!givenDjs) {
587      for (iRow=0;iRow<numberRows_;iRow++) {
588        int iPivot=pivotVariable_[iRow];
589        double value = cost_[iPivot];
590        if (value) {
591          array[iRow]=value;
592          index[number++]=iRow;
593        }
594      }
595    } else {
596      // dual values pass - djs may not be zero
597      for (iRow=0;iRow<numberRows_;iRow++) {
598        int iPivot=pivotVariable_[iRow];
599        // make sure zero if done
600        if (!pivoted(iPivot))
601          givenDjs[iPivot]=0.0;
602        double value =cost_[iPivot]-givenDjs[iPivot];
603        if (value) {
604          array[iRow]=value;
605          index[number++]=iRow;
606        }
607      }
608    }
609    arrayVector.setNumElements(number);
610    // Extended duals before "updateTranspose"
611    matrix_->dualExpanded(this,&arrayVector,givenDjs,0);
612   
613    // Btran basic costs and get as accurate as possible
614    double lastError=COIN_DBL_MAX;
615    int iRefine;
616    double * work = workSpace->denseVector();
617    CoinIndexedVector * thisVector = &arrayVector;
618    CoinIndexedVector * lastVector = &previousVector;
619    factorization_->updateColumnTranspose(workSpace,thisVector);
620   
621    for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
622      // check basic reduced costs zero
623      largestDualError_=0.0;
624      // would be faster to do just for basic but this reduces code
625      ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
626      transposeTimes(-1.0,array,reducedCostWork_);
627      // update by duals on sets
628      matrix_->dualExpanded(this,NULL,NULL,1);
629      if (!givenDjs) {
630        for (iRow=0;iRow<numberRows_;iRow++) {
631          int iPivot=pivotVariable_[iRow];
632          double value;
633          if (iPivot>=numberColumns_) {
634            // slack
635            value = rowObjectiveWork_[iPivot-numberColumns_]
636              + array[iPivot-numberColumns_];
637          } else {
638            // column
639            value = reducedCostWork_[iPivot];
640          }
641          work[iRow]=value;
642          if (fabs(value)>largestDualError_) {
643            largestDualError_=fabs(value);
644          }
645        }
646      } else {
647        for (iRow=0;iRow<numberRows_;iRow++) {
648          int iPivot=pivotVariable_[iRow];
649          if (iPivot>=numberColumns_) {
650            // slack
651            work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
652              + array[iPivot-numberColumns_]-givenDjs[iPivot];
653          } else {
654            // column
655            work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
656          }
657          if (fabs(work[iRow])>largestDualError_) {
658            largestDualError_=fabs(work[iRow]);
659            //assert (largestDualError_<1.0e-7);
660            //if (largestDualError_>1.0e-7)
661            //printf("large dual error %g\n",largestDualError_);
662          }
663        }
664      }
665      if (largestDualError_>=lastError) {
666        // restore
667        CoinIndexedVector * temp = thisVector;
668        thisVector = lastVector;
669        lastVector=temp;
670        break;
671      }
672      if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
673          &&!givenDjs) {
674        // try and make better
675        // save this
676        CoinIndexedVector * temp = thisVector;
677        thisVector = lastVector;
678        lastVector=temp;
679        int * indexOut = thisVector->getIndices();
680        int number=0;
681        array = thisVector->denseVector();
682        thisVector->clear();
683        double multiplier = 131072.0;
684        for (iRow=0;iRow<numberRows_;iRow++) {
685          double value = multiplier*work[iRow];
686          if (value) {
687            array[iRow]=value;
688            indexOut[number++]=iRow;
689            work[iRow]=0.0;
690          }
691          work[iRow]=0.0;
692        }
693        thisVector->setNumElements(number);
694        lastError=largestDualError_;
695        factorization_->updateColumnTranspose(workSpace,thisVector);
696        multiplier = 1.0/multiplier;
697        double * previous = lastVector->denseVector();
698        number=0;
699        for (iRow=0;iRow<numberRows_;iRow++) {
700          double value = previous[iRow] + multiplier*array[iRow];
701          if (value) {
702            array[iRow]=value;
703            indexOut[number++]=iRow;
704          } else {
705            array[iRow]=0.0;
706          }
707        }
708        thisVector->setNumElements(number);
709      } else {
710        break;
711      }
712    }
713    ClpFillN(work,numberRows_,0.0);
714    // now look at dual solution
715    array = thisVector->denseVector();
716    for (iRow=0;iRow<numberRows_;iRow++) {
717      // slack
718      double value = array[iRow];
719      dual_[iRow]=value;
720      value += rowObjectiveWork_[iRow];
721      rowReducedCost_[iRow]=value;
722    }
723    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
724    transposeTimes(-1.0,dual_,reducedCostWork_);
725    // Extended duals and check dual infeasibility
726    if (!matrix_->skipDualCheck()||algorithm_<0||problemStatus_!=-2) 
727      matrix_->dualExpanded(this,NULL,NULL,2);
728    // If necessary - override results
729    if (givenDjs) {
730      // restore accurate duals
731      memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
732    }
733  } else {
734    // Nonlinear
735    objective_->reducedGradient(this,dj_,false);
736  }
737}
738/* Given an existing factorization computes and checks
739   primal and dual solutions.  Uses input arrays for variables at
740   bounds.  Returns feasibility states */
741int ClpSimplex::getSolution ( const double * rowActivities,
742                               const double * columnActivities)
743{
744  if (!factorization_->status()) {
745    // put in standard form
746    createRim(7+8+16+32);
747    // do work
748    gutsOfSolution ( NULL,NULL);
749    // release extra memory
750    deleteRim(0);
751  }
752  return factorization_->status();
753}
754/* Given an existing factorization computes and checks
755   primal and dual solutions.  Uses current problem arrays for
756   bounds.  Returns feasibility states */
757int ClpSimplex::getSolution ( )
758{
759  double * rowActivities = new double[numberRows_];
760  double * columnActivities = new double[numberColumns_];
761  ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
762  ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
763  int status = getSolution( rowActivities, columnActivities);
764  delete [] rowActivities;
765  delete [] columnActivities;
766  return status;
767}
768// Factorizes using current basis.  This is for external use
769// Return codes are as from ClpFactorization
770int ClpSimplex::factorize ()
771{
772  // put in standard form
773  createRim(7+8+16+32,false);
774  // do work
775  int status = internalFactorize(-1);
776  // release extra memory
777  deleteRim(0);
778
779  return status;
780}
781// Clean up status
782void 
783ClpSimplex::cleanStatus()
784{
785  int iRow,iColumn;
786  int numberBasic=0;
787  // make row activities correct
788  memset(rowActivityWork_,0,numberRows_*sizeof(double));
789  times(1.0,columnActivityWork_,rowActivityWork_);
790  if (!status_)
791    createStatus();
792  for (iRow=0;iRow<numberRows_;iRow++) {
793    if (getRowStatus(iRow)==basic) 
794      numberBasic++;
795    else {
796      setRowStatus(iRow,superBasic);
797      // but put to bound if close
798      if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
799          <=primalTolerance_) {
800        rowActivityWork_[iRow]=rowLowerWork_[iRow];
801        setRowStatus(iRow,atLowerBound);
802      } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
803                 <=primalTolerance_) {
804        rowActivityWork_[iRow]=rowUpperWork_[iRow];
805        setRowStatus(iRow,atUpperBound);
806      }
807    }
808  }
809  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
810    if (getColumnStatus(iColumn)==basic) {
811      if (numberBasic==numberRows_) {
812        // take out of basis
813        setColumnStatus(iColumn,superBasic);
814        // but put to bound if close
815        if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
816            <=primalTolerance_) {
817          columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
818          setColumnStatus(iColumn,atLowerBound);
819        } else if (fabs(columnActivityWork_[iColumn]
820                        -columnUpperWork_[iColumn])
821                   <=primalTolerance_) {
822          columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
823          setColumnStatus(iColumn,atUpperBound);
824        }
825      } else 
826        numberBasic++;
827    } else {
828      setColumnStatus(iColumn,superBasic);
829      // but put to bound if close
830      if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
831          <=primalTolerance_) {
832        columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
833        setColumnStatus(iColumn,atLowerBound);
834      } else if (fabs(columnActivityWork_[iColumn]
835                      -columnUpperWork_[iColumn])
836                 <=primalTolerance_) {
837        columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
838        setColumnStatus(iColumn,atUpperBound);
839      }
840    }
841  }
842}
843
844/* Factorizes using current basis. 
845      solveType - 1 iterating, 0 initial, -1 external
846      - 2 then iterating but can throw out of basis
847      If 10 added then in primal values pass
848*/
849/* Return codes are as from ClpFactorization unless initial factorization
850   when total number of singularities is returned */
851int ClpSimplex::internalFactorize ( int solveType)
852{
853  int iRow,iColumn;
854  int totalSlacks=numberRows_;
855  if (!status_)
856    createStatus();
857
858  bool valuesPass=false;
859  if (solveType>=10) {
860    valuesPass=true;
861    solveType -= 10;
862  }
863#ifdef CLP_DEBUG
864  if (solveType>0) {
865    int numberFreeIn=0,numberFreeOut=0;
866    double biggestDj=0.0;
867    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
868      switch(getColumnStatus(iColumn)) {
869
870      case basic:
871        if (columnLower_[iColumn]<-largeValue_
872            &&columnUpper_[iColumn]>largeValue_) 
873          numberFreeIn++;
874        break;
875      default:
876        if (columnLower_[iColumn]<-largeValue_
877            &&columnUpper_[iColumn]>largeValue_) {
878          numberFreeOut++;
879          biggestDj = CoinMax(fabs(dj_[iColumn]),biggestDj);
880        }
881        break;
882      }
883    }
884    if (numberFreeIn+numberFreeOut)
885      printf("%d in basis, %d out - largest dj %g\n",
886             numberFreeIn,numberFreeOut,biggestDj);
887  }
888#endif
889  if (solveType<=0) {
890    // Make sure everything is clean
891    for (iRow=0;iRow<numberRows_;iRow++) {
892      if(getRowStatus(iRow)==isFixed) {
893        // double check fixed
894        if (rowUpperWork_[iRow]>rowLowerWork_[iRow])
895          setRowStatus(iRow,atLowerBound);
896      } else if (getRowStatus(iRow)==isFree) {
897        // may not be free after all
898        if (rowLowerWork_[iRow]>-largeValue_||rowUpperWork_[iRow]<largeValue_)
899          setRowStatus(iRow,superBasic);
900      }
901    }
902    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
903      if(getColumnStatus(iColumn)==isFixed) {
904        // double check fixed
905        if (columnUpperWork_[iColumn]>columnLowerWork_[iColumn])
906          setColumnStatus(iColumn,atLowerBound);
907      } else if (getColumnStatus(iColumn)==isFree) {
908        // may not be free after all
909        if (columnLowerWork_[iColumn]>-largeValue_||columnUpperWork_[iColumn]<largeValue_)
910          setColumnStatus(iColumn,superBasic);
911      }
912    }
913    if (!valuesPass) {
914      // not values pass so set to bounds
915      bool allSlack=true;
916      if (status_) {
917        for (iRow=0;iRow<numberRows_;iRow++) {
918          if (getRowStatus(iRow)!=basic) {
919            allSlack=false;
920            break;
921          }
922        }
923      }
924      if (!allSlack) {
925        // set values from warm start (if sensible)
926        int numberBasic=0;
927        for (iRow=0;iRow<numberRows_;iRow++) {
928          switch(getRowStatus(iRow)) {
929           
930          case basic:
931            numberBasic++;
932            break;
933          case atUpperBound:
934            rowActivityWork_[iRow]=rowUpperWork_[iRow];
935            if (rowActivityWork_[iRow]>largeValue_) {
936              if (rowLowerWork_[iRow]>-largeValue_) {
937                rowActivityWork_[iRow]=rowLowerWork_[iRow];
938                setRowStatus(iRow,atLowerBound);
939              } else {
940                // say free
941                setRowStatus(iRow,isFree);
942                rowActivityWork_[iRow]=0.0;
943              }
944            }
945            break;
946          case ClpSimplex::isFixed:
947          case atLowerBound:
948            rowActivityWork_[iRow]=rowLowerWork_[iRow];
949            if (rowActivityWork_[iRow]<-largeValue_) {
950              if (rowUpperWork_[iRow]<largeValue_) {
951                rowActivityWork_[iRow]=rowUpperWork_[iRow];
952                setRowStatus(iRow,atUpperBound);
953              } else {
954                // say free
955                setRowStatus(iRow,isFree);
956                rowActivityWork_[iRow]=0.0;
957              }
958            }
959            break;
960          case isFree:
961              break;
962            // not really free - fall through to superbasic
963          case superBasic:
964            if (rowUpperWork_[iRow]>largeValue_) {
965              if (rowLowerWork_[iRow]>-largeValue_) {
966                rowActivityWork_[iRow]=rowLowerWork_[iRow];
967                setRowStatus(iRow,atLowerBound);
968              } else {
969                // say free
970                setRowStatus(iRow,isFree);
971                rowActivityWork_[iRow]=0.0;
972              }
973            } else {
974              if (rowLowerWork_[iRow]>-largeValue_) {
975                // set to nearest
976                if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
977                    <fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])) {
978                  rowActivityWork_[iRow]=rowLowerWork_[iRow];
979                  setRowStatus(iRow,atLowerBound);
980                } else {
981                  rowActivityWork_[iRow]=rowUpperWork_[iRow];
982                  setRowStatus(iRow,atUpperBound);
983                }
984              } else {
985                rowActivityWork_[iRow]=rowUpperWork_[iRow];
986                setRowStatus(iRow,atUpperBound);
987              }
988            }
989            break;
990          }
991        }
992        totalSlacks=numberBasic;
993
994        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
995          switch(getColumnStatus(iColumn)) {
996           
997          case basic:
998            if (numberBasic==maximumBasic_) {
999              // take out of basis
1000              if (columnLowerWork_[iColumn]>-largeValue_) {
1001                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
1002                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
1003                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1004                  setColumnStatus(iColumn,atLowerBound);
1005                } else {
1006                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1007                  setColumnStatus(iColumn,atUpperBound);
1008                }
1009              } else if (columnUpperWork_[iColumn]<largeValue_) {
1010                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1011                setColumnStatus(iColumn,atUpperBound);
1012              } else {
1013                columnActivityWork_[iColumn]=0.0;
1014                setColumnStatus(iColumn,isFree);
1015              }
1016            } else {
1017              numberBasic++;
1018            }
1019            break;
1020          case atUpperBound:
1021            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1022            if (columnActivityWork_[iColumn]>largeValue_) {
1023              if (columnLowerWork_[iColumn]<-largeValue_) {
1024                columnActivityWork_[iColumn]=0.0;
1025                setColumnStatus(iColumn,isFree);
1026              } else {
1027                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1028                setColumnStatus(iColumn,atLowerBound);
1029              }
1030            }
1031            break;
1032          case isFixed:
1033          case atLowerBound:
1034            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1035            if (columnActivityWork_[iColumn]<-largeValue_) {
1036              if (columnUpperWork_[iColumn]>largeValue_) {
1037                columnActivityWork_[iColumn]=0.0;
1038                setColumnStatus(iColumn,isFree);
1039              } else {
1040                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1041                setColumnStatus(iColumn,atUpperBound);
1042              }
1043            }
1044            break;
1045          case isFree:
1046              break;
1047            // not really free - fall through to superbasic
1048          case superBasic:
1049            if (columnUpperWork_[iColumn]>largeValue_) {
1050              if (columnLowerWork_[iColumn]>-largeValue_) {
1051                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1052                setColumnStatus(iColumn,atLowerBound);
1053              } else {
1054                // say free
1055                setColumnStatus(iColumn,isFree);
1056                columnActivityWork_[iColumn]=0.0;
1057              }
1058            } else {
1059              if (columnLowerWork_[iColumn]>-largeValue_) {
1060                // set to nearest
1061                if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
1062                    <fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])) {
1063                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1064                  setColumnStatus(iColumn,atLowerBound);
1065                } else {
1066                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1067                  setColumnStatus(iColumn,atUpperBound);
1068                }
1069              } else {
1070                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1071                setColumnStatus(iColumn,atUpperBound);
1072              }
1073            }
1074            break;
1075          }
1076        }
1077      } else {
1078        // all slack basis
1079        int numberBasic=0;
1080        if (!status_) {
1081          createStatus();
1082        }
1083        for (iRow=0;iRow<numberRows_;iRow++) {
1084          double lower=rowLowerWork_[iRow];
1085          double upper=rowUpperWork_[iRow];
1086          if (lower>-largeValue_||upper<largeValue_) {
1087            if (fabs(lower)<=fabs(upper)) {
1088              rowActivityWork_[iRow]=lower;
1089            } else {
1090              rowActivityWork_[iRow]=upper;
1091            }
1092          } else {
1093            rowActivityWork_[iRow]=0.0;
1094          }
1095          setRowStatus(iRow,basic);
1096          numberBasic++;
1097        }
1098        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1099          double lower=columnLowerWork_[iColumn];
1100          double upper=columnUpperWork_[iColumn];
1101          double big_bound = largeValue_;
1102          if (lower>-big_bound||upper<big_bound) {
1103            if ((getColumnStatus(iColumn)==atLowerBound&&
1104                 columnActivityWork_[iColumn]==lower)||
1105                (getColumnStatus(iColumn)==atUpperBound&&
1106                 columnActivityWork_[iColumn]==upper)) {
1107              // status looks plausible
1108            } else {
1109              // set to sensible
1110              if (fabs(lower)<=fabs(upper)) {
1111                setColumnStatus(iColumn,atLowerBound);
1112                columnActivityWork_[iColumn]=lower;
1113              } else {
1114                setColumnStatus(iColumn,atUpperBound);
1115                columnActivityWork_[iColumn]=upper;
1116              }
1117            }
1118          } else {
1119            setColumnStatus(iColumn,isFree);
1120            columnActivityWork_[iColumn]=0.0;
1121          }
1122        }
1123      }
1124    } else {
1125      // values pass has less coding
1126      // make row activities correct and clean basis a bit
1127      cleanStatus();
1128      if (status_) {
1129        int numberBasic=0;
1130        for (iRow=0;iRow<numberRows_;iRow++) {
1131          if (getRowStatus(iRow)==basic) 
1132            numberBasic++;
1133        }
1134        totalSlacks=numberBasic;
1135        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1136          if (getColumnStatus(iColumn)==basic) 
1137            numberBasic++;
1138        }
1139      } else {
1140        // all slack basis
1141        int numberBasic=0;
1142        if (!status_) {
1143          createStatus();
1144        }
1145        for (iRow=0;iRow<numberRows_;iRow++) {
1146          setRowStatus(iRow,basic);
1147          numberBasic++;
1148        }
1149        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1150          setColumnStatus(iColumn,superBasic);
1151          // but put to bound if close
1152          if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
1153              <=primalTolerance_) {
1154            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1155            setColumnStatus(iColumn,atLowerBound);
1156          } else if (fabs(columnActivityWork_[iColumn]
1157                        -columnUpperWork_[iColumn])
1158                     <=primalTolerance_) {
1159            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
1160            setColumnStatus(iColumn,atUpperBound);
1161          }
1162        }
1163      }
1164    }
1165    numberRefinements_=1;
1166    // set fixed if they are
1167    for (iRow=0;iRow<numberRows_;iRow++) {
1168      if (getRowStatus(iRow)!=basic ) {
1169        if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
1170          rowActivityWork_[iRow]=rowLowerWork_[iRow];
1171          setRowStatus(iRow,isFixed);
1172        }
1173      }
1174    }
1175    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1176      if (getColumnStatus(iColumn)!=basic ) {
1177        if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
1178          columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1179          setColumnStatus(iColumn,isFixed);
1180        }
1181      }
1182    }
1183  }
1184  if (0)  {
1185    static int k=0;
1186    printf("start basis\n");
1187    int i;
1188    for (i=0;i<numberRows_;i++)
1189      printf ("xx %d %d\n",i,pivotVariable_[i]);
1190    for (i=0;i<numberRows_+numberColumns_;i++)
1191      if (getColumnStatus(i)==basic)
1192        printf ("yy %d basic\n",i);
1193    if (k>20)
1194      exit(0);
1195    k++;
1196  }
1197  int status = factorization_->factorize(this, solveType,valuesPass);
1198  if (status) {
1199    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
1200      <<status
1201      <<CoinMessageEol;
1202    return -1;
1203  } else if (!solveType) {
1204    // Initial basis - return number of singularities
1205    int numberSlacks=0;
1206    for (iRow=0;iRow<numberRows_;iRow++) {
1207      if (getRowStatus(iRow) == basic)
1208        numberSlacks++;
1209    }
1210    status= CoinMax(numberSlacks-totalSlacks,0);
1211  }
1212
1213  // sparse methods
1214  //if (factorization_->sparseThreshold()) {
1215    // get default value
1216    factorization_->sparseThreshold(0);
1217    factorization_->goSparse();
1218    //}
1219 
1220  return status;
1221}
1222/*
1223   This does basis housekeeping and does values for in/out variables.
1224   Can also decide to re-factorize
1225*/
1226int 
1227ClpSimplex::housekeeping(double objectiveChange)
1228{
1229  // save value of incoming and outgoing
1230  double oldIn = solution_[sequenceIn_];
1231  double oldOut = solution_[sequenceOut_];
1232  numberIterations_++;
1233  changeMade_++; // something has happened
1234  // incoming variable
1235  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
1236    <<directionOut_
1237    <<directionIn_<<theta_
1238    <<dualOut_<<dualIn_<<alpha_
1239    <<CoinMessageEol;
1240  if (getStatus(sequenceIn_)==isFree) {
1241    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
1242      <<sequenceIn_
1243      <<CoinMessageEol;
1244  }
1245  // change of incoming
1246  char rowcol[]={'R','C'};
1247  if (pivotRow_>=0)
1248    pivotVariable_[pivotRow_]=sequenceIn();
1249  if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20)
1250    progressFlag_ |= 2; // making real progress
1251  solution_[sequenceIn_]=valueIn_;
1252  if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
1253    progressFlag_ |= 1; // making real progress
1254  if (sequenceIn_!=sequenceOut_) {
1255    //assert( getStatus(sequenceOut_)== basic);
1256    setStatus(sequenceIn_,basic);
1257    if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
1258      // As Nonlinear costs may have moved bounds (to more feasible)
1259      // Redo using value
1260      if (fabs(valueOut_-lower_[sequenceOut_])<fabs(valueOut_-upper_[sequenceOut_])) {
1261        // going to lower
1262        setStatus(sequenceOut_,atLowerBound);
1263        oldOut = lower_[sequenceOut_];
1264      } else {
1265        // going to upper
1266        setStatus(sequenceOut_,atUpperBound);
1267        oldOut = upper_[sequenceOut_];
1268      }
1269    } else {
1270      // fixed
1271      setStatus(sequenceOut_,isFixed);
1272    }
1273    solution_[sequenceOut_]=valueOut_;
1274  } else {
1275    //if (objective_->type()<2)
1276    assert (fabs(theta_)>1.0e-13);
1277    // flip from bound to bound
1278    // As Nonlinear costs may have moved bounds (to more feasible)
1279    // Redo using value
1280    if (fabs(valueIn_-lower_[sequenceIn_])<fabs(valueIn_-upper_[sequenceIn_])) {
1281      // as if from upper bound
1282      setStatus(sequenceIn_, atLowerBound);
1283    } else {
1284      // as if from lower bound
1285      setStatus(sequenceIn_, atUpperBound);
1286    }
1287  }
1288 
1289  // Update hidden stuff e.g. effective RHS and gub
1290  matrix_->updatePivot(this,oldIn,oldOut);
1291  objectiveValue_ += objectiveChange/(objectiveScale_*rhsScale_);
1292  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
1293    <<numberIterations_<<objectiveValue()
1294    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
1295    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
1296  handler_->printing(algorithm_<0)<<dualOut_<<theta_;
1297  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
1298  handler_->message()<<CoinMessageEol;
1299  if (hitMaximumIterations())
1300    return 2;
1301#if 1
1302  //if (numberIterations_>14000)
1303  //handler_->setLogLevel(63);
1304  //if (numberIterations_>24000)
1305  //exit(77);
1306  // check for small cycles
1307  int in = sequenceIn_;
1308  int out = sequenceOut_;
1309  matrix_->correctSequence(in,out);
1310  int cycle=progress_->cycle(in,out,
1311                            directionIn_,directionOut_);
1312  if (cycle>0&&objective_->type()<2) {
1313    //if (cycle>0) {
1314    if (handler_->logLevel()>=63)
1315      printf("Cycle of %d\n",cycle);
1316    // reset
1317    progress_->startCheck();
1318    if (factorization_->pivots()>cycle) {
1319      forceFactorization_=cycle-1;
1320    } else {
1321      // need to reject something
1322      int iSequence;
1323      if (algorithm_<0)
1324        iSequence = sequenceIn_;
1325      else
1326        iSequence = sequenceOut_;
1327      char x = isColumn(iSequence) ? 'C' :'R';
1328      if (handler_->logLevel()>=63)
1329        handler_->message(CLP_SIMPLEX_FLAG,messages_)
1330          <<x<<sequenceWithin(iSequence)
1331          <<CoinMessageEol;
1332      setFlagged(iSequence);
1333      //printf("flagging %d\n",iSequence);
1334    }
1335    return 1;
1336  }
1337#endif
1338  // only time to re-factorize if one before real time
1339  // this is so user won't be surprised that maximumPivots has exact meaning
1340  if (factorization_->pivots()==factorization_->maximumPivots()||
1341      factorization_->maximumPivots()<2) {
1342    return 1;
1343  } else {
1344    if (forceFactorization_>0&&
1345        factorization_->pivots()==forceFactorization_) {
1346      // relax
1347      forceFactorization_ = (3+5*forceFactorization_)/4;
1348      if (forceFactorization_>factorization_->maximumPivots())
1349        forceFactorization_ = -1; //off
1350      return 1;
1351    } else {
1352      // carry on iterating
1353      return 0;
1354    }
1355  }
1356}
1357// Copy constructor.
1358ClpSimplex::ClpSimplex(const ClpSimplex &rhs,int scalingMode) :
1359  ClpModel(rhs,scalingMode),
1360  columnPrimalInfeasibility_(0.0),
1361  rowPrimalInfeasibility_(0.0),
1362  columnPrimalSequence_(-2),
1363  rowPrimalSequence_(-2), 
1364  columnDualInfeasibility_(0.0),
1365  rowDualInfeasibility_(0.0),
1366  columnDualSequence_(-2),
1367  rowDualSequence_(-2),
1368  primalToleranceToGetOptimal_(-1.0),
1369  remainingDualInfeasibility_(0.0),
1370  largeValue_(1.0e15),
1371  largestPrimalError_(0.0),
1372  largestDualError_(0.0),
1373  largestSolutionError_(0.0),
1374  dualBound_(1.0e10),
1375  alpha_(0.0),
1376  theta_(0.0),
1377  lowerIn_(0.0),
1378  valueIn_(0.0),
1379  upperIn_(0.0),
1380  dualIn_(0.0),
1381  lowerOut_(-1),
1382  valueOut_(-1),
1383  upperOut_(-1),
1384  dualOut_(-1),
1385  dualTolerance_(0.0),
1386  primalTolerance_(0.0),
1387  sumDualInfeasibilities_(0.0),
1388  sumPrimalInfeasibilities_(0.0),
1389  infeasibilityCost_(1.0e10),
1390  sumOfRelaxedDualInfeasibilities_(0.0),
1391  sumOfRelaxedPrimalInfeasibilities_(0.0),
1392  lower_(NULL),
1393  rowLowerWork_(NULL),
1394  columnLowerWork_(NULL),
1395  upper_(NULL),
1396  rowUpperWork_(NULL),
1397  columnUpperWork_(NULL),
1398  cost_(NULL),
1399  rowObjectiveWork_(NULL),
1400  objectiveWork_(NULL),
1401  sequenceIn_(-1),
1402  directionIn_(-1),
1403  sequenceOut_(-1),
1404  directionOut_(-1),
1405  pivotRow_(-1),
1406  lastGoodIteration_(-100),
1407  dj_(NULL),
1408  rowReducedCost_(NULL),
1409  reducedCostWork_(NULL),
1410  solution_(NULL),
1411  rowActivityWork_(NULL),
1412  columnActivityWork_(NULL),
1413  numberDualInfeasibilities_(0),
1414  numberDualInfeasibilitiesWithoutFree_(0),
1415  numberPrimalInfeasibilities_(100),
1416  numberRefinements_(0),
1417  pivotVariable_(NULL),
1418  factorization_(NULL),
1419  savedSolution_(NULL),
1420  numberTimesOptimal_(0),
1421  changeMade_(1),
1422  algorithm_(0),
1423  forceFactorization_(-1),
1424  perturbation_(100),
1425  nonLinearCost_(NULL),
1426  specialOptions_(0),
1427  lastBadIteration_(-999999),
1428  lastFlaggedIteration_(-999999),
1429  numberFake_(0),
1430  progressFlag_(0),
1431  firstFree_(-1),
1432  numberExtraRows_(0),
1433  maximumBasic_(0),
1434  incomingInfeasibility_(1.0),
1435  allowedInfeasibility_(10.0),
1436  automaticScale_(0),
1437  progress_(NULL)
1438{
1439  int i;
1440  for (i=0;i<6;i++) {
1441    rowArray_[i]=NULL;
1442    columnArray_[i]=NULL;
1443  }
1444  saveStatus_=NULL;
1445  factorization_ = NULL;
1446  dualRowPivot_ = NULL;
1447  primalColumnPivot_ = NULL;
1448  gutsOfDelete(0);
1449  specialOptions_ =0;
1450  delete nonLinearCost_;
1451  nonLinearCost_ = NULL;
1452  gutsOfCopy(rhs);
1453  solveType_=1; // say simplex based life form
1454}
1455// Copy constructor from model
1456ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
1457  ClpModel(rhs,scalingMode),
1458  columnPrimalInfeasibility_(0.0),
1459  rowPrimalInfeasibility_(0.0),
1460  columnPrimalSequence_(-2),
1461  rowPrimalSequence_(-2), 
1462  columnDualInfeasibility_(0.0),
1463  rowDualInfeasibility_(0.0),
1464  columnDualSequence_(-2),
1465  rowDualSequence_(-2),
1466  primalToleranceToGetOptimal_(-1.0),
1467  remainingDualInfeasibility_(0.0),
1468  largeValue_(1.0e15),
1469  largestPrimalError_(0.0),
1470  largestDualError_(0.0),
1471  largestSolutionError_(0.0),
1472  dualBound_(1.0e10),
1473  alpha_(0.0),
1474  theta_(0.0),
1475  lowerIn_(0.0),
1476  valueIn_(0.0),
1477  upperIn_(0.0),
1478  dualIn_(0.0),
1479  lowerOut_(-1),
1480  valueOut_(-1),
1481  upperOut_(-1),
1482  dualOut_(-1),
1483  dualTolerance_(0.0),
1484  primalTolerance_(0.0),
1485  sumDualInfeasibilities_(0.0),
1486  sumPrimalInfeasibilities_(0.0),
1487  infeasibilityCost_(1.0e10),
1488  sumOfRelaxedDualInfeasibilities_(0.0),
1489  sumOfRelaxedPrimalInfeasibilities_(0.0),
1490  lower_(NULL),
1491  rowLowerWork_(NULL),
1492  columnLowerWork_(NULL),
1493  upper_(NULL),
1494  rowUpperWork_(NULL),
1495  columnUpperWork_(NULL),
1496  cost_(NULL),
1497  rowObjectiveWork_(NULL),
1498  objectiveWork_(NULL),
1499  sequenceIn_(-1),
1500  directionIn_(-1),
1501  sequenceOut_(-1),
1502  directionOut_(-1),
1503  pivotRow_(-1),
1504  lastGoodIteration_(-100),
1505  dj_(NULL),
1506  rowReducedCost_(NULL),
1507  reducedCostWork_(NULL),
1508  solution_(NULL),
1509  rowActivityWork_(NULL),
1510  columnActivityWork_(NULL),
1511  numberDualInfeasibilities_(0),
1512  numberDualInfeasibilitiesWithoutFree_(0),
1513  numberPrimalInfeasibilities_(100),
1514  numberRefinements_(0),
1515  pivotVariable_(NULL),
1516  factorization_(NULL),
1517  savedSolution_(NULL),
1518  numberTimesOptimal_(0),
1519  changeMade_(1),
1520  algorithm_(0),
1521  forceFactorization_(-1),
1522  perturbation_(100),
1523  nonLinearCost_(NULL),
1524  specialOptions_(0),
1525  lastBadIteration_(-999999),
1526  lastFlaggedIteration_(-999999),
1527  numberFake_(0),
1528  progressFlag_(0),
1529  firstFree_(-1),
1530  numberExtraRows_(0),
1531  maximumBasic_(0),
1532  incomingInfeasibility_(1.0),
1533  allowedInfeasibility_(10.0),
1534  automaticScale_(0),
1535  progress_(NULL)
1536{
1537  int i;
1538  for (i=0;i<6;i++) {
1539    rowArray_[i]=NULL;
1540    columnArray_[i]=NULL;
1541  }
1542  saveStatus_=NULL;
1543  // get an empty factorization so we can set tolerances etc
1544  factorization_ = new ClpFactorization();
1545  // say Steepest pricing
1546  dualRowPivot_ = new ClpDualRowSteepest();
1547  // say Steepest pricing
1548  primalColumnPivot_ = new ClpPrimalColumnSteepest();
1549  solveType_=1; // say simplex based life form
1550 
1551}
1552// Assignment operator. This copies the data
1553ClpSimplex & 
1554ClpSimplex::operator=(const ClpSimplex & rhs)
1555{
1556  if (this != &rhs) {
1557    gutsOfDelete(0);
1558    specialOptions_=0;
1559    delete nonLinearCost_;
1560    nonLinearCost_ = NULL;
1561    ClpModel::operator=(rhs);
1562    gutsOfCopy(rhs);
1563  }
1564  return *this;
1565}
1566void 
1567ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1568{
1569  assert (numberRows_==rhs.numberRows_);
1570  assert (numberColumns_==rhs.numberColumns_);
1571  numberExtraRows_ = rhs.numberExtraRows_;
1572  maximumBasic_ = rhs.maximumBasic_;
1573  int numberRows2 = numberRows_+numberExtraRows_;
1574  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows2);
1575  rowLowerWork_ = lower_+numberColumns_;
1576  columnLowerWork_ = lower_;
1577  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows2);
1578  rowUpperWork_ = upper_+numberColumns_;
1579  columnUpperWork_ = upper_;
1580  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
1581  cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows2));
1582  objectiveWork_ = cost_;
1583  rowObjectiveWork_ = cost_+numberColumns_;
1584  dj_ = ClpCopyOfArray(rhs.dj_,numberRows2+numberColumns_);
1585  if (dj_) {
1586    reducedCostWork_ = dj_;
1587    rowReducedCost_ = dj_+numberColumns_;
1588  }
1589  solution_ = ClpCopyOfArray(rhs.solution_,numberRows2+numberColumns_);
1590  if (solution_) {
1591    columnActivityWork_ = solution_;
1592    rowActivityWork_ = solution_+numberColumns_;
1593  }
1594  if (rhs.pivotVariable_) {
1595    pivotVariable_ = new int[numberRows2];
1596    CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
1597  } else {
1598    pivotVariable_=NULL;
1599  }
1600  if (rhs.factorization_) {
1601    factorization_ = new ClpFactorization(*rhs.factorization_);
1602  } else {
1603    factorization_=NULL;
1604  }
1605  savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1606  int i;
1607  for (i=0;i<6;i++) {
1608    rowArray_[i]=NULL;
1609    if (rhs.rowArray_[i]) 
1610      rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
1611    columnArray_[i]=NULL;
1612    if (rhs.columnArray_[i]) 
1613      columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
1614  }
1615  if (rhs.saveStatus_) {
1616    saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1617  }
1618  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1619  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1620  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1621  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1622  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1623  columnDualSequence_ = rhs.columnDualSequence_;
1624  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1625  rowDualSequence_ = rhs.rowDualSequence_;
1626  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1627  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1628  largeValue_ = rhs.largeValue_;
1629  largestPrimalError_ = rhs.largestPrimalError_;
1630  largestDualError_ = rhs.largestDualError_;
1631  largestSolutionError_ = rhs.largestSolutionError_;
1632  dualBound_ = rhs.dualBound_;
1633  alpha_ = rhs.alpha_;
1634  theta_ = rhs.theta_;
1635  lowerIn_ = rhs.lowerIn_;
1636  valueIn_ = rhs.valueIn_;
1637  upperIn_ = rhs.upperIn_;
1638  dualIn_ = rhs.dualIn_;
1639  sequenceIn_ = rhs.sequenceIn_;
1640  directionIn_ = rhs.directionIn_;
1641  lowerOut_ = rhs.lowerOut_;
1642  valueOut_ = rhs.valueOut_;
1643  upperOut_ = rhs.upperOut_;
1644  dualOut_ = rhs.dualOut_;
1645  sequenceOut_ = rhs.sequenceOut_;
1646  directionOut_ = rhs.directionOut_;
1647  pivotRow_ = rhs.pivotRow_;
1648  lastGoodIteration_ = rhs.lastGoodIteration_;
1649  numberRefinements_ = rhs.numberRefinements_;
1650  dualTolerance_ = rhs.dualTolerance_;
1651  primalTolerance_ = rhs.primalTolerance_;
1652  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1653  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1654  numberDualInfeasibilitiesWithoutFree_ = 
1655    rhs.numberDualInfeasibilitiesWithoutFree_;
1656  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1657  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1658  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1659  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1660  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1661  changeMade_ = rhs.changeMade_;
1662  algorithm_ = rhs.algorithm_;
1663  forceFactorization_ = rhs.forceFactorization_;
1664  perturbation_ = rhs.perturbation_;
1665  infeasibilityCost_ = rhs.infeasibilityCost_;
1666  specialOptions_ = rhs.specialOptions_;
1667  lastBadIteration_ = rhs.lastBadIteration_;
1668  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
1669  numberFake_ = rhs.numberFake_;
1670  progressFlag_ = rhs.progressFlag_;
1671  firstFree_ = rhs.firstFree_;
1672  incomingInfeasibility_ = rhs.incomingInfeasibility_;
1673  allowedInfeasibility_ = rhs.allowedInfeasibility_;
1674  automaticScale_ = rhs.automaticScale_;
1675  if (rhs.progress_)
1676    progress_ = new ClpSimplexProgress(*rhs.progress_);
1677  else
1678    progress_=NULL;
1679  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
1680  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
1681  if (rhs.nonLinearCost_!=NULL)
1682    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1683  else
1684    nonLinearCost_=NULL;
1685  solveType_=rhs.solveType_;
1686}
1687// type == 0 do everything, most + pivot data, 2 factorization data as well
1688void 
1689ClpSimplex::gutsOfDelete(int type)
1690{
1691  delete [] lower_;
1692  lower_=NULL;
1693  rowLowerWork_=NULL;
1694  columnLowerWork_=NULL;
1695  delete [] upper_;
1696  upper_=NULL;
1697  rowUpperWork_=NULL;
1698  columnUpperWork_=NULL;
1699  delete [] cost_;
1700  cost_=NULL;
1701  objectiveWork_=NULL;
1702  rowObjectiveWork_=NULL;
1703  delete [] dj_;
1704  dj_=NULL;
1705  reducedCostWork_=NULL;
1706  rowReducedCost_=NULL;
1707  delete [] solution_;
1708  solution_=NULL;
1709  rowActivityWork_=NULL;
1710  columnActivityWork_=NULL;
1711  delete [] savedSolution_;
1712  savedSolution_ = NULL;
1713  if ((specialOptions_&2)==0) {
1714    delete nonLinearCost_;
1715    nonLinearCost_ = NULL;
1716  }
1717  int i;
1718  for (i=0;i<6;i++) {
1719    delete rowArray_[i];
1720    rowArray_[i]=NULL;
1721    delete columnArray_[i];
1722    columnArray_[i]=NULL;
1723  }
1724  delete rowCopy_;
1725  rowCopy_=NULL;
1726  delete [] saveStatus_;
1727  saveStatus_=NULL;
1728  if (!type) {
1729    // delete everything
1730    delete factorization_;
1731    factorization_ = NULL;
1732    delete [] pivotVariable_;
1733    pivotVariable_=NULL;
1734    delete dualRowPivot_;
1735    dualRowPivot_ = NULL;
1736    delete primalColumnPivot_;
1737    primalColumnPivot_ = NULL;
1738    delete progress_;
1739    progress_=NULL;
1740  } else {
1741    // delete any size information in methods
1742    if (type>1) {
1743      factorization_->clearArrays();
1744      delete [] pivotVariable_;
1745      pivotVariable_=NULL;
1746    }
1747    dualRowPivot_->clearArrays();
1748    primalColumnPivot_->clearArrays();
1749  }
1750}
1751// This sets largest infeasibility and most infeasible
1752void 
1753ClpSimplex::checkPrimalSolution(const double * rowActivities,
1754                                        const double * columnActivities)
1755{
1756  double * solution;
1757  int iRow,iColumn;
1758
1759  objectiveValue_ = 0.0;
1760  // now look at primal solution
1761  columnPrimalInfeasibility_=0.0;
1762  columnPrimalSequence_=-1;
1763  rowPrimalInfeasibility_=0.0;
1764  rowPrimalSequence_=-1;
1765  largestSolutionError_=0.0;
1766  solution = rowActivityWork_;
1767  sumPrimalInfeasibilities_=0.0;
1768  numberPrimalInfeasibilities_=0;
1769  double primalTolerance = primalTolerance_;
1770  double relaxedTolerance=primalTolerance_;
1771  // we can't really trust infeasibilities if there is primal error
1772  double error = CoinMin(1.0e-3,largestPrimalError_);
1773  // allow tolerance at least slightly bigger than standard
1774  relaxedTolerance = relaxedTolerance +  error;
1775  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
1776
1777  for (iRow=0;iRow<numberRows_;iRow++) {
1778    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
1779    double infeasibility=0.0;
1780    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1781    if (solution[iRow]>rowUpperWork_[iRow]) {
1782      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1783    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1784      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1785    }
1786    if (infeasibility>primalTolerance) {
1787      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1788      if (infeasibility>relaxedTolerance) 
1789        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1790      numberPrimalInfeasibilities_ ++;
1791    }
1792    if (infeasibility>rowPrimalInfeasibility_) {
1793      rowPrimalInfeasibility_=infeasibility;
1794      rowPrimalSequence_=iRow;
1795    }
1796    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1797    if (infeasibility>largestSolutionError_)
1798      largestSolutionError_=infeasibility;
1799  }
1800  // Check any infeasibilities from dynamic rows
1801  matrix_->primalExpanded(this,2);
1802  solution = columnActivityWork_;
1803  if (!matrix_->rhsOffset(this)) {
1804    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1805      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
1806      double infeasibility=0.0;
1807      objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1808      if (solution[iColumn]>columnUpperWork_[iColumn]) {
1809        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1810      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1811        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1812      }
1813      if (infeasibility>columnPrimalInfeasibility_) {
1814        columnPrimalInfeasibility_=infeasibility;
1815        columnPrimalSequence_=iColumn;
1816      }
1817      if (infeasibility>primalTolerance) {
1818        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1819        if (infeasibility>relaxedTolerance) 
1820          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1821        numberPrimalInfeasibilities_ ++;
1822      }
1823      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1824      if (infeasibility>largestSolutionError_)
1825        largestSolutionError_=infeasibility;
1826    }
1827  } else {
1828    // as we are using effective rhs we only check basics
1829    // But we do need to get objective
1830    objectiveValue_ += innerProduct(objectiveWork_,numberColumns_,solution);
1831    for (int j=0;j<numberRows_;j++) {
1832      int iColumn = pivotVariable_[j];
1833      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
1834      double infeasibility=0.0;
1835      if (solution[iColumn]>columnUpperWork_[iColumn]) {
1836        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1837      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1838        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1839      }
1840      if (infeasibility>columnPrimalInfeasibility_) {
1841        columnPrimalInfeasibility_=infeasibility;
1842        columnPrimalSequence_=iColumn;
1843      }
1844      if (infeasibility>primalTolerance) {
1845        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1846        if (infeasibility>relaxedTolerance) 
1847          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1848        numberPrimalInfeasibilities_ ++;
1849      }
1850      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1851      if (infeasibility>largestSolutionError_)
1852        largestSolutionError_=infeasibility;
1853    }
1854  }
1855  objectiveValue_ += objective_->nonlinearOffset();
1856  objectiveValue_ /= (objectiveScale_*rhsScale_);
1857}
1858void 
1859ClpSimplex::checkDualSolution()
1860{
1861
1862  int iRow,iColumn;
1863  sumDualInfeasibilities_=0.0;
1864  numberDualInfeasibilities_=0;
1865  numberDualInfeasibilitiesWithoutFree_=0;
1866  columnDualInfeasibility_=0.0;
1867  columnDualSequence_=-1;
1868  if (matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2) {
1869    // pretend we found dual infeasibilities
1870    sumOfRelaxedDualInfeasibilities_ = 1.0;
1871    sumDualInfeasibilities_=1.0;
1872    numberDualInfeasibilities_=1;
1873    return;
1874  }
1875  rowDualInfeasibility_=0.0;
1876  rowDualSequence_=-1;
1877  int firstFreePrimal = -1;
1878  int firstFreeDual = -1;
1879  int numberSuperBasicWithDj=0;
1880  primalToleranceToGetOptimal_=CoinMax(rowPrimalInfeasibility_,
1881                                   columnPrimalInfeasibility_);
1882  remainingDualInfeasibility_=0.0;
1883  double relaxedTolerance=dualTolerance_;
1884  // we can't really trust infeasibilities if there is dual error
1885  double error = CoinMin(1.0e-3,largestDualError_);
1886  // allow tolerance at least slightly bigger than standard
1887  relaxedTolerance = relaxedTolerance +  error;
1888  sumOfRelaxedDualInfeasibilities_ = 0.0;
1889
1890  // Check any djs from dynamic rows
1891  matrix_->dualExpanded(this,NULL,NULL,3);
1892  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
1893  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1894    if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
1895      // not basic
1896      double distanceUp = columnUpperWork_[iColumn]-
1897        columnActivityWork_[iColumn];
1898      double distanceDown = columnActivityWork_[iColumn] -
1899        columnLowerWork_[iColumn];
1900      if (distanceUp>primalTolerance_) {
1901        double value = reducedCostWork_[iColumn];
1902        // Check if "free"
1903        if (distanceDown>primalTolerance_) {
1904          if (fabs(value)>1.0e2*relaxedTolerance) {
1905            numberSuperBasicWithDj++;
1906            if (firstFreeDual<0)
1907              firstFreeDual = iColumn;
1908          }
1909          if (firstFreePrimal<0)
1910            firstFreePrimal = iColumn;
1911        }
1912        // should not be negative
1913        if (value<0.0) {
1914          value = - value;
1915          if (value>columnDualInfeasibility_) {
1916            columnDualInfeasibility_=value;
1917            columnDualSequence_=iColumn;
1918          }
1919          if (value>dualTolerance_) {
1920            if (getColumnStatus(iColumn) != isFree) {
1921              numberDualInfeasibilitiesWithoutFree_ ++;
1922              sumDualInfeasibilities_ += value-dualTolerance_;
1923              if (value>relaxedTolerance) 
1924                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1925              numberDualInfeasibilities_ ++;
1926            } else {
1927              // free so relax a lot
1928              value *= 0.01;
1929              if (value>dualTolerance_) {
1930                sumDualInfeasibilities_ += value-dualTolerance_;
1931                if (value>relaxedTolerance) 
1932                  sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1933                numberDualInfeasibilities_ ++;
1934              }
1935            }
1936            // maybe we can make feasible by increasing tolerance
1937            if (distanceUp<largeValue_) {
1938              if (distanceUp>primalToleranceToGetOptimal_)
1939                primalToleranceToGetOptimal_=distanceUp;
1940            } else {
1941              //gap too big for any tolerance
1942              remainingDualInfeasibility_=
1943                CoinMax(remainingDualInfeasibility_,value);
1944            }
1945          }
1946        }
1947      }
1948      if (distanceDown>primalTolerance_) {
1949        double value = reducedCostWork_[iColumn];
1950        // should not be positive
1951        if (value>0.0) {
1952          if (value>columnDualInfeasibility_) {
1953            columnDualInfeasibility_=value;
1954            columnDualSequence_=iColumn;
1955          }
1956          if (value>dualTolerance_) {
1957            sumDualInfeasibilities_ += value-dualTolerance_;
1958            if (value>relaxedTolerance) 
1959              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1960            numberDualInfeasibilities_ ++;
1961            if (getColumnStatus(iColumn) != isFree) 
1962              numberDualInfeasibilitiesWithoutFree_ ++;
1963            // maybe we can make feasible by increasing tolerance
1964            if (distanceDown<largeValue_&&
1965                distanceDown>primalToleranceToGetOptimal_)
1966              primalToleranceToGetOptimal_=distanceDown;
1967          }
1968        }
1969      }
1970    }
1971  }
1972  for (iRow=0;iRow<numberRows_;iRow++) {
1973    if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
1974      // not basic
1975      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
1976      double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
1977      if (distanceUp>primalTolerance_) {
1978        double value = rowReducedCost_[iRow];
1979        // Check if "free"
1980        if (distanceDown>primalTolerance_) {
1981          if (fabs(value)>1.0e2*relaxedTolerance) {
1982            numberSuperBasicWithDj++;
1983            if (firstFreeDual<0)
1984              firstFreeDual = iRow+numberColumns_;
1985          }
1986          if (firstFreePrimal<0)
1987            firstFreePrimal = iRow+numberColumns_;
1988        }
1989        // should not be negative
1990        if (value<0.0) {
1991          value = - value;
1992          if (value>rowDualInfeasibility_) {
1993            rowDualInfeasibility_=value;
1994            rowDualSequence_=iRow;
1995          }
1996          if (value>dualTolerance_) {
1997            sumDualInfeasibilities_ += value-dualTolerance_;
1998            if (value>relaxedTolerance) 
1999              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2000            numberDualInfeasibilities_ ++;
2001            if (getRowStatus(iRow) != isFree) 
2002              numberDualInfeasibilitiesWithoutFree_ ++;
2003            // maybe we can make feasible by increasing tolerance
2004            if (distanceUp<largeValue_) {
2005              if (distanceUp>primalToleranceToGetOptimal_)
2006                primalToleranceToGetOptimal_=distanceUp;
2007            } else {
2008              //gap too big for any tolerance
2009              remainingDualInfeasibility_=
2010                CoinMax(remainingDualInfeasibility_,value);
2011            }
2012          }
2013        }
2014      }
2015      if (distanceDown>primalTolerance_) {
2016        double value = rowReducedCost_[iRow];
2017        // should not be positive
2018        if (value>0.0) {
2019          if (value>rowDualInfeasibility_) {
2020            rowDualInfeasibility_=value;
2021            rowDualSequence_=iRow;
2022          }
2023          if (value>dualTolerance_) {
2024            sumDualInfeasibilities_ += value-dualTolerance_;
2025            if (value>relaxedTolerance) 
2026              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2027            numberDualInfeasibilities_ ++;
2028            if (getRowStatus(iRow) != isFree) 
2029              numberDualInfeasibilitiesWithoutFree_ ++;
2030            // maybe we can make feasible by increasing tolerance
2031            if (distanceDown<largeValue_&&
2032                distanceDown>primalToleranceToGetOptimal_)
2033              primalToleranceToGetOptimal_=distanceDown;
2034          }
2035        }
2036      }
2037    }
2038  }
2039  if (algorithm_<0&&firstFreeDual>=0) {
2040    // dual
2041    firstFree_ = firstFreeDual;
2042  } else if (numberSuperBasicWithDj||
2043             (progress_&&progress_->lastIterationNumber(0)<=0)) {
2044    firstFree_=firstFreePrimal;
2045  }
2046}
2047/* Adds multiple of a column into an array */
2048void 
2049ClpSimplex::add(double * array,
2050                int sequence, double multiplier) const
2051{
2052  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2053    //slack
2054    array [sequence-numberColumns_] -= multiplier;
2055  } else {
2056    // column
2057    matrix_->add(this,array,sequence,multiplier);
2058  }
2059}
2060/*
2061  Unpacks one column of the matrix into indexed array
2062*/
2063void 
2064ClpSimplex::unpack(CoinIndexedVector * rowArray) const
2065{
2066  rowArray->clear();
2067  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2068    //slack
2069    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
2070  } else {
2071    // column
2072    matrix_->unpack(this,rowArray,sequenceIn_);
2073  }
2074}
2075void 
2076ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence) const
2077{
2078  rowArray->clear();
2079  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2080    //slack
2081    rowArray->insert(sequence-numberColumns_,-1.0);
2082  } else {
2083    // column
2084    matrix_->unpack(this,rowArray,sequence);
2085  }
2086}
2087/*
2088  Unpacks one column of the matrix into indexed array
2089*/
2090void 
2091ClpSimplex::unpackPacked(CoinIndexedVector * rowArray) 
2092{
2093  rowArray->clear();
2094  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2095    //slack
2096    int * index = rowArray->getIndices();
2097    double * array = rowArray->denseVector();
2098    array[0]=-1.0;
2099    index[0]=sequenceIn_-numberColumns_;
2100    rowArray->setNumElements(1);
2101    rowArray->setPackedMode(true);
2102  } else {
2103    // column
2104    matrix_->unpackPacked(this,rowArray,sequenceIn_);
2105  }
2106}
2107void 
2108ClpSimplex::unpackPacked(CoinIndexedVector * rowArray,int sequence)
2109{
2110  rowArray->clear();
2111  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2112    //slack
2113    int * index = rowArray->getIndices();
2114    double * array = rowArray->denseVector();
2115    array[0]=-1.0;
2116    index[0]=sequence-numberColumns_;
2117    rowArray->setNumElements(1);
2118    rowArray->setPackedMode(true);
2119  } else {
2120    // column
2121    matrix_->unpackPacked(this,rowArray,sequence);
2122  }
2123}
2124bool
2125ClpSimplex::createRim(int what,bool makeRowCopy)
2126{
2127  bool goodMatrix=true;
2128  int saveLevel=handler_->logLevel();
2129  if (problemStatus_==10) {
2130    handler_->setLogLevel(0); // switch off messages
2131  } else if (factorization_) {
2132    // match up factorization messages
2133    if (handler_->logLevel()<3)
2134      factorization_->messageLevel(0);
2135    else
2136      factorization_->messageLevel(CoinMax(3,factorization_->messageLevel()));
2137  }
2138  numberExtraRows_ = matrix_->generalExpanded(this,2,maximumBasic_);
2139  if (numberExtraRows_) {
2140    // make sure status array large enough
2141    assert (status_);
2142    int numberOld = numberRows_+numberColumns_;
2143    int numberNew = numberRows_+numberColumns_+numberExtraRows_;
2144    unsigned char * newStatus = new unsigned char [numberNew];
2145    memset(newStatus+numberOld,0,numberExtraRows_);
2146    memcpy(newStatus,status_,numberOld);
2147    delete [] status_;
2148    status_=newStatus;
2149  }
2150  int numberRows2 = numberRows_+numberExtraRows_;
2151  int i;
2152  if ((what&1)!=0) {
2153    delete [] lower_;
2154    delete [] upper_;
2155    lower_ = new double[numberColumns_+numberRows2];
2156    upper_ = new double[numberColumns_+numberRows2];
2157    rowLowerWork_ = lower_+numberColumns_;
2158    columnLowerWork_ = lower_;
2159    rowUpperWork_ = upper_+numberColumns_;
2160    columnUpperWork_ = upper_;
2161    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
2162    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
2163    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
2164    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
2165    // clean up any mismatches on infinity
2166    for (i=0;i<numberColumns_;i++) {
2167      if (columnLowerWork_[i]<-1.0e30)
2168        columnLowerWork_[i] = -COIN_DBL_MAX;
2169      if (columnUpperWork_[i]>1.0e30)
2170        columnUpperWork_[i] = COIN_DBL_MAX;
2171    }
2172    // clean up any mismatches on infinity
2173    for (i=0;i<numberRows_;i++) {
2174      if (rowLowerWork_[i]<-1.0e30)
2175        rowLowerWork_[i] = -COIN_DBL_MAX;
2176      if (rowUpperWork_[i]>1.0e30)
2177        rowUpperWork_[i] = COIN_DBL_MAX;
2178    }
2179  }
2180  if ((what&(16+32))!=0) {
2181    // move information to work arrays
2182    double direction = optimizationDirection_;
2183    // direction is actually scale out not scale in
2184    if (direction)
2185      direction = 1.0/direction;
2186    if (direction!=1.0) {
2187      // reverse all dual signs
2188      for (i=0;i<numberColumns_;i++) 
2189        reducedCost_[i] *= direction;
2190      for (i=0;i<numberRows_;i++) 
2191        dual_[i] *= direction;
2192    }
2193    // row reduced costs
2194    if (!dj_) {
2195      dj_ = new double[numberRows2+numberColumns_];
2196      reducedCostWork_ = dj_;
2197      rowReducedCost_ = dj_+numberColumns_;
2198      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
2199      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
2200    }
2201    if (!solution_||(what&32)!=0) {
2202      if (!solution_)
2203        solution_ = new double[numberRows2+numberColumns_];
2204      columnActivityWork_ = solution_;
2205      rowActivityWork_ = solution_+numberColumns_;
2206      if (status_) {
2207        for (i=0;i<numberColumns_;i++) {
2208          columnActivityWork_[i] = columnActivity_[i];
2209          setFakeBound(i,noFake);
2210        }
2211        for (i=0;i<numberRows_;i++) {
2212          rowActivityWork_[i] = rowActivity_[i];
2213          setFakeBound(i+numberColumns_,noFake);
2214        }
2215      } else {
2216        for (i=0;i<numberColumns_;i++) {
2217          columnActivityWork_[i] = columnActivity_[i];
2218        }
2219        for (i=0;i<numberRows_;i++) {
2220          rowActivityWork_[i] = rowActivity_[i];
2221        }
2222      }
2223      for (i=0;i<numberColumns_+numberRows_;i++) {
2224        if (fabs(solution_[i])>1.0e20) {
2225          //printf("bad value of %g for column %d\n",solution_[i],i);
2226          setColumnStatus(i,superBasic);
2227          if (upper_[i]<0.0) {
2228            solution_[i]=upper_[i];
2229          } else if (lower_[i]>0.0) {
2230            solution_[i]=lower_[i];
2231          } else {
2232            solution_[i]=0.0;
2233          }
2234        }
2235      }
2236    }
2237  }
2238  if ((what&16)!=0) {
2239    //check matrix
2240    if (!matrix_)
2241      matrix_=new ClpPackedMatrix();
2242    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20)) {
2243      problemStatus_=4;
2244      goodMatrix= false;
2245    }
2246    if (makeRowCopy) {
2247      delete rowCopy_;
2248      // may return NULL if can't give row copy
2249      rowCopy_ = matrix_->reverseOrderedCopy();
2250#ifdef TAKEOUT
2251      {
2252
2253        ClpPackedMatrix* rowCopy =
2254          dynamic_cast< ClpPackedMatrix*>(rowCopy_);
2255        const int * column = rowCopy->getIndices();
2256        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2257        const double * element = rowCopy->getElements();
2258        int i;
2259        for (i=133;i<numberRows_;i++) {
2260          if (rowStart[i+1]-rowStart[i]==10||rowStart[i+1]-rowStart[i]==15)
2261            printf("Row %d has %d elements\n",i,rowStart[i+1]-rowStart[i]);
2262        }
2263      } 
2264#endif
2265    }
2266  }
2267  if ((what&4)!=0) {
2268    delete [] cost_;
2269    // extra copy with original costs
2270    int nTotal = numberRows2+numberColumns_;
2271    //cost_ = new double[2*nTotal];
2272    cost_ = new double[nTotal];
2273    objectiveWork_ = cost_;
2274    rowObjectiveWork_ = cost_+numberColumns_;
2275    memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double));
2276    if (rowObjective_)
2277      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
2278    else
2279      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2280    // and initialize changes to zero
2281    //memset(cost_+nTotal,0,nTotal*sizeof(double));
2282  }
2283  // do scaling if needed
2284  if (scalingFlag_>0&&!rowScale_) {
2285    if (matrix_->scale(this))
2286      scalingFlag_=-scalingFlag_; // not scaled after all
2287    if (rowScale_&&automaticScale_) {
2288      // try automatic scaling
2289      double smallestObj=1.0e100;
2290      double largestObj=0.0;
2291      double largestRhs=0.0;
2292      for (i=0;i<numberColumns_+numberRows_;i++) {
2293        double value = fabs(cost_[i]);
2294        if (i<numberColumns_)
2295          value *= columnScale_[i];
2296        if (value&&lower_[i]!=upper_[i]) {
2297          smallestObj = CoinMin(smallestObj,value);
2298          largestObj = CoinMax(largestObj,value);
2299        }
2300        if (lower_[i]>0.0||upper_[i]<0.0) {
2301          double scale;
2302          if (i<numberColumns_)
2303            scale = 1.0/columnScale_[i];
2304          else
2305            scale = rowScale_[i-numberColumns_];
2306          //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2307          if (lower_[i]>0)
2308            largestRhs=CoinMax(largestRhs,lower_[i]*scale);
2309          if (upper_[i]<0.0)
2310            largestRhs=CoinMax(largestRhs,-upper_[i]*scale);
2311        }
2312      }
2313      printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2314      bool scalingDone=false;
2315      // look at element range
2316      double smallestNegative;
2317      double largestNegative;
2318      double smallestPositive;
2319      double largestPositive;
2320      matrix_->rangeOfElements(smallestNegative, largestNegative,
2321                               smallestPositive, largestPositive);
2322      smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2323      largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2324      if (largestObj) {
2325        double ratio = largestObj/smallestObj;
2326        double scale=1.0;
2327        if (ratio<1.0e8) {
2328          // reasonable
2329          if (smallestObj<1.0e-4) {
2330            // may as well scale up
2331            scalingDone=true;
2332            scale=1.0e-3/smallestObj;
2333          } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2334            //done=true;
2335          } else {
2336            scalingDone=true;
2337            if (algorithm_<0) {
2338              scale = 1.0e6/largestObj;
2339            } else {
2340              scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2341            }
2342          }
2343        } else if (ratio<1.0e12) {
2344          // not so good
2345          if (smallestObj<1.0e-7) {
2346            // may as well scale up
2347            scalingDone=true;
2348            scale=1.0e-6/smallestObj;
2349          } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2350            //done=true;
2351          } else {
2352            scalingDone=true;
2353            if (algorithm_<0) {
2354              scale = 1.0e7/largestObj;
2355            } else {
2356              scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2357            }
2358          }
2359        } else {
2360          // Really nasty problem
2361          if (smallestObj<1.0e-8) {
2362            // may as well scale up
2363            scalingDone=true;
2364            scale=1.0e-7/smallestObj;
2365            largestObj *= scale;
2366          } 
2367          if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2368            //done=true;
2369          } else {
2370            scalingDone=true;
2371            if (algorithm_<0) {
2372              scale = 1.0e7/largestObj;
2373            } else {
2374              scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2375            }
2376          }
2377        }
2378        objectiveScale_=scale;
2379      }
2380      if (largestRhs>1.0e12) {
2381        scalingDone=true;
2382        rhsScale_=1.0e9/largestRhs;
2383      } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
2384        scalingDone=true;
2385        rhsScale_=1.0e6/largestRhs;
2386      } else {
2387        rhsScale_=1.0;
2388      }
2389      if (scalingDone) {
2390        handler_->message(CLP_RIM_SCALE,messages_)
2391          <<objectiveScale_<<rhsScale_
2392          <<CoinMessageEol;
2393      }
2394    }
2395  } else if (makeRowCopy&&(what&16)!=0&&scalingFlag_>0) {
2396    matrix_->scaleRowCopy(this);
2397  }
2398  if ((what&4)!=0) {
2399    double direction = optimizationDirection_*objectiveScale_;
2400    // but also scale by scale factors
2401    // not really sure about row scaling
2402    if (!rowScale_) {
2403      if (direction!=1.0) {
2404        for (i=0;i<numberRows_;i++)
2405          rowObjectiveWork_[i] *= direction;
2406        for (i=0;i<numberColumns_;i++)
2407          objectiveWork_[i] *= direction;
2408      }
2409    } else {
2410      for (i=0;i<numberRows_;i++)
2411        rowObjectiveWork_[i] *= direction/rowScale_[i];
2412      for (i=0;i<numberColumns_;i++)
2413        objectiveWork_[i] *= direction*columnScale_[i];
2414    }
2415  }
2416  if ((what&(1+32))!=0) {
2417    if(rowScale_) {
2418      for (i=0;i<numberColumns_;i++) {
2419        double multiplier = rhsScale_/columnScale_[i];
2420        if (columnLowerWork_[i]>-1.0e50)
2421          columnLowerWork_[i] *= multiplier;
2422        if (columnUpperWork_[i]<1.0e50)
2423          columnUpperWork_[i] *= multiplier;
2424       
2425      }
2426      for (i=0;i<numberRows_;i++) {
2427        double multiplier = rhsScale_*rowScale_[i];
2428        if (rowLowerWork_[i]>-1.0e50)
2429          rowLowerWork_[i] *= multiplier;
2430        if (rowUpperWork_[i]<1.0e50)
2431          rowUpperWork_[i] *= multiplier;
2432      }
2433    } else if (rhsScale_!=1.0) {
2434      for (i=0;i<numberColumns_+numberRows_;i++) {
2435        if (lower_[i]>-1.0e50)
2436          lower_[i] *= rhsScale_;
2437        if (upper_[i]<1.0e50)
2438          upper_[i] *= rhsScale_;
2439       
2440      }
2441    }
2442  }
2443  if ((what&(8+32))!=0) {
2444    if (rowScale_) {
2445      // on entry
2446      for (i=0;i<numberColumns_;i++) {
2447        columnActivityWork_[i] /= columnScale_[i];
2448        columnActivityWork_[i] *= rhsScale_;
2449        reducedCostWork_[i] *= objectiveScale_*columnScale_[i];
2450      }
2451      for (i=0;i<numberRows_;i++) {
2452        rowActivityWork_[i] *= rowScale_[i]*rhsScale_;
2453        dual_[i] /= rowScale_[i];
2454        dual_[i] *= objectiveScale_;
2455        rowReducedCost_[i] = dual_[i];
2456      }
2457    } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
2458      // on entry
2459      for (i=0;i<numberColumns_;i++) {
2460        columnActivityWork_[i] *= rhsScale_;
2461        reducedCostWork_[i] *= objectiveScale_;
2462      }
2463      for (i=0;i<numberRows_;i++) {
2464        rowActivityWork_[i] *= rhsScale_;
2465        dual_[i] *= objectiveScale_;
2466        rowReducedCost_[i] = dual_[i];
2467      }
2468    }
2469  }
2470 
2471  if ((what&16)!=0) {
2472    // check rim of problem okay
2473    if (!sanityCheck())
2474      goodMatrix=false;
2475  } 
2476  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
2477  // maybe we need to move scales to SimplexModel for factorization?
2478  if ((what&8)!=0&&!pivotVariable_) {
2479    pivotVariable_=new int[numberRows2];
2480    for (int i=0;i<numberRows2;i++)
2481      pivotVariable_[i]=-1;
2482  }
2483
2484  if ((what&16)!=0&&!rowArray_[2]) {
2485    // get some arrays
2486    int iRow,iColumn;
2487    // these are "indexed" arrays so we always know where nonzeros are
2488    /**********************************************************
2489      rowArray_[3] is long enough for rows+columns
2490    *********************************************************/
2491    for (iRow=0;iRow<4;iRow++) {
2492      if (!rowArray_[iRow]) {
2493        rowArray_[iRow]=new CoinIndexedVector();
2494        int length =numberRows2+factorization_->maximumPivots();
2495        if (iRow==3||objective_->type()>1)
2496          length += numberColumns_;
2497        rowArray_[iRow]->reserve(length);
2498      }
2499    }
2500   
2501    for (iColumn=0;iColumn<2;iColumn++) {
2502      if (!columnArray_[iColumn]) {
2503        columnArray_[iColumn]=new CoinIndexedVector();
2504        if (!iColumn)
2505          columnArray_[iColumn]->reserve(numberColumns_);
2506        else
2507          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
2508      }
2509    }   
2510  }
2511  double primalTolerance=dblParam_[ClpPrimalTolerance];
2512  if ((what&1)!=0) {
2513    // fix any variables with tiny gaps
2514    for (i=0;i<numberColumns_;i++) {
2515      if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2516        if (columnLowerWork_[i]>=0.0) {
2517          columnUpperWork_[i] = columnLowerWork_[i];
2518        } else if (columnUpperWork_[i]<=0.0) {
2519          columnLowerWork_[i] = columnUpperWork_[i];
2520        } else {
2521          columnUpperWork_[i] = 0.0;
2522          columnLowerWork_[i] = 0.0;
2523        }
2524      }
2525    }
2526    for (i=0;i<numberRows_;i++) {
2527      if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2528        if (rowLowerWork_[i]>=0.0) {
2529          rowUpperWork_[i] = rowLowerWork_[i];
2530        } else if (rowUpperWork_[i]<=0.0) {
2531          rowLowerWork_[i] = rowUpperWork_[i];
2532        } else {
2533          rowUpperWork_[i] = 0.0;
2534          rowLowerWork_[i] = 0.0;
2535        }
2536      }
2537    }
2538  }
2539  if (problemStatus_==10) {
2540    problemStatus_=-1;
2541    handler_->setLogLevel(saveLevel); // switch back messages
2542  }
2543  if ((what&5)!=0) 
2544    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
2545  return goodMatrix;
2546}
2547void
2548ClpSimplex::deleteRim(int getRidOfFactorizationData)
2549{
2550  int i;
2551  if (problemStatus_!=1&&problemStatus_!=2) {
2552    delete [] ray_;
2553    ray_=NULL;
2554  }
2555  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
2556  upperOut_=1.0;
2557  // ray may be null if in branch and bound
2558  if (rowScale_) {
2559    // Collect infeasibilities
2560    int numberPrimalScaled=0;
2561    int numberPrimalUnscaled=0;
2562    int numberDualScaled=0;
2563    int numberDualUnscaled=0;
2564    double scaleC = 1.0/objectiveScale_;
2565    double scaleR = 1.0/rhsScale_;
2566    for (i=0;i<numberColumns_;i++) {
2567      double scaleFactor = columnScale_[i];
2568      double valueScaled = columnActivityWork_[i];
2569      double lowerScaled = columnLowerWork_[i];
2570      double upperScaled = columnUpperWork_[i];
2571      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2572        if (valueScaled<lowerScaled-primalTolerance_||
2573            valueScaled>upperScaled+primalTolerance_)
2574          numberPrimalScaled++;
2575        else
2576          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2577      }
2578      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
2579      double value = columnActivity_[i];
2580      if (value<columnLower_[i]-primalTolerance_)
2581        numberPrimalUnscaled++;
2582      else if (value>columnUpper_[i]+primalTolerance_)
2583        numberPrimalUnscaled++;
2584      double valueScaledDual = reducedCostWork_[i];
2585      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2586        numberDualScaled++;
2587      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2588        numberDualScaled++;
2589      reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
2590      double valueDual = reducedCost_[i];
2591      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2592        numberDualUnscaled++;
2593      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2594        numberDualUnscaled++;
2595    }
2596    for (i=0;i<numberRows_;i++) {
2597      double scaleFactor = rowScale_[i];
2598      double valueScaled = rowActivityWork_[i];
2599      double lowerScaled = rowLowerWork_[i];
2600      double upperScaled = rowUpperWork_[i];
2601      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2602        if (valueScaled<lowerScaled-primalTolerance_||
2603            valueScaled>upperScaled+primalTolerance_)
2604          numberPrimalScaled++;
2605        else
2606          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2607      }
2608      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
2609      double value = rowActivity_[i];
2610      if (value<rowLower_[i]-primalTolerance_)
2611        numberPrimalUnscaled++;
2612      else if (value>rowUpper_[i]+primalTolerance_)
2613        numberPrimalUnscaled++;
2614      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2615      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2616        numberDualScaled++;
2617      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2618        numberDualScaled++;
2619      dual_[i] = valueScaledDual*scaleFactor*scaleC;
2620      double valueDual = dual_[i]; 
2621      if (rowObjective_)
2622        valueDual += rowObjective_[i];
2623      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2624        numberDualUnscaled++;
2625      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2626        numberDualUnscaled++;
2627    }
2628    if (!problemStatus_&&!secondaryStatus_) {
2629      // See if we need to set secondary status
2630      if (numberPrimalUnscaled) {
2631        if (numberDualUnscaled) 
2632          secondaryStatus_=4;
2633        else
2634          secondaryStatus_=2;
2635      } else {
2636        if (numberDualUnscaled) 
2637          secondaryStatus_=3;
2638      }
2639    }
2640    if (problemStatus_==2) {
2641      for (i=0;i<numberColumns_;i++) {
2642        ray_[i] *= columnScale_[i];
2643      }
2644    } else if (problemStatus_==1&&ray_) {
2645      for (i=0;i<numberRows_;i++) {
2646        ray_[i] *= rowScale_[i];
2647      }
2648    }
2649  } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
2650    // Collect infeasibilities
2651    int numberPrimalScaled=0;
2652    int numberPrimalUnscaled=0;
2653    int numberDualScaled=0;
2654    int numberDualUnscaled=0;
2655    double scaleC = 1.0/objectiveScale_;
2656    double scaleR = 1.0/rhsScale_;
2657    for (i=0;i<numberColumns_;i++) {
2658      double valueScaled = columnActivityWork_[i];
2659      double lowerScaled = columnLowerWork_[i];
2660      double upperScaled = columnUpperWork_[i];
2661      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2662        if (valueScaled<lowerScaled-primalTolerance_||
2663            valueScaled>upperScaled+primalTolerance_)
2664          numberPrimalScaled++;
2665        else
2666          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2667      }
2668      columnActivity_[i] = valueScaled*scaleR;
2669      double value = columnActivity_[i];
2670      if (value<columnLower_[i]-primalTolerance_)
2671        numberPrimalUnscaled++;
2672      else if (value>columnUpper_[i]+primalTolerance_)
2673        numberPrimalUnscaled++;
2674      double valueScaledDual = reducedCostWork_[i];
2675      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2676        numberDualScaled++;
2677      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2678        numberDualScaled++;
2679      reducedCost_[i] = valueScaledDual*scaleC;
2680      double valueDual = reducedCost_[i];
2681      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2682        numberDualUnscaled++;
2683      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2684        numberDualUnscaled++;
2685    }
2686    for (i=0;i<numberRows_;i++) {
2687      double valueScaled = rowActivityWork_[i];
2688      double lowerScaled = rowLowerWork_[i];
2689      double upperScaled = rowUpperWork_[i];
2690      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2691        if (valueScaled<lowerScaled-primalTolerance_||
2692            valueScaled>upperScaled+primalTolerance_)
2693          numberPrimalScaled++;
2694        else
2695          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2696      }
2697      rowActivity_[i] = valueScaled*scaleR;
2698      double value = rowActivity_[i];
2699      if (value<rowLower_[i]-primalTolerance_)
2700        numberPrimalUnscaled++;
2701      else if (value>rowUpper_[i]+primalTolerance_)
2702        numberPrimalUnscaled++;
2703      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2704      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2705        numberDualScaled++;
2706      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2707        numberDualScaled++;
2708      dual_[i] = valueScaledDual*scaleC;
2709      double valueDual = dual_[i]; 
2710      if (rowObjective_)
2711        valueDual += rowObjective_[i];
2712      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2713        numberDualUnscaled++;
2714      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2715        numberDualUnscaled++;
2716    }
2717    if (!problemStatus_&&!secondaryStatus_) {
2718      // See if we need to set secondary status
2719      if (numberPrimalUnscaled) {
2720        if (numberDualUnscaled) 
2721          secondaryStatus_=4;
2722        else
2723          secondaryStatus_=2;
2724      } else {
2725        if (numberDualUnscaled) 
2726          secondaryStatus_=3;
2727      }
2728    }
2729  } else {
2730    if (columnActivityWork_) {
2731      for (i=0;i<numberColumns_;i++) {
2732        double value = columnActivityWork_[i];
2733        double lower = columnLowerWork_[i];
2734        double upper = columnUpperWork_[i];
2735        if (lower>-1.0e20||upper<1.0e20) {
2736          if (value>lower&&value<upper)
2737            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2738        }
2739        columnActivity_[i] = columnActivityWork_[i];
2740        reducedCost_[i] = reducedCostWork_[i];
2741      }
2742      for (i=0;i<numberRows_;i++) {
2743        double value = rowActivityWork_[i];
2744        double lower = rowLowerWork_[i];
2745        double upper = rowUpperWork_[i];
2746        if (lower>-1.0e20||upper<1.0e20) {
2747          if (value>lower&&value<upper)
2748            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2749        }
2750        rowActivity_[i] = rowActivityWork_[i];
2751      }
2752    }
2753  }
2754  // switch off scalefactor if auto
2755  if (automaticScale_) {
2756    rhsScale_=1.0;
2757    objectiveScale_=1.0;
2758  }
2759  if (optimizationDirection_!=1.0) {
2760    // and modify all dual signs
2761    for (i=0;i<numberColumns_;i++) 
2762      reducedCost_[i] *= optimizationDirection_;
2763    for (i=0;i<numberRows_;i++) 
2764      dual_[i] *= optimizationDirection_;
2765  }
2766  // scaling may have been turned off
2767  scalingFlag_ = abs(scalingFlag_);
2768  if(getRidOfFactorizationData>=0)
2769    gutsOfDelete(getRidOfFactorizationData+1);
2770  // get rid of data
2771  matrix_->generalExpanded(this,13,scalingFlag_);
2772}
2773void 
2774ClpSimplex::setDualBound(double value)
2775{
2776  if (value>0.0)
2777    dualBound_=value;
2778}
2779void 
2780ClpSimplex::setInfeasibilityCost(double value)
2781{
2782  if (value>0.0)
2783    infeasibilityCost_=value;
2784}
2785void ClpSimplex::setNumberRefinements( int value) 
2786{
2787  if (value>=0&&value<10)
2788    numberRefinements_=value;
2789}
2790// Sets row pivot choice algorithm in dual
2791void 
2792ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
2793{
2794  delete dualRowPivot_;
2795  dualRowPivot_ = choice.clone(true);
2796}
2797// Sets row pivot choice algorithm in dual
2798void 
2799ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
2800{
2801  delete primalColumnPivot_;
2802  primalColumnPivot_ = choice.clone(true);
2803}
2804// Passes in factorization
2805void 
2806ClpSimplex::setFactorization( ClpFactorization & factorization)
2807{
2808  delete factorization_;
2809  factorization_= new ClpFactorization(factorization);
2810}
2811/* Perturbation:
2812   -50 to +50 - perturb by this power of ten (-6 sounds good)
2813   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
2814   101 - we are perturbed
2815   102 - don't try perturbing again
2816   default is 100
2817*/
2818void 
2819ClpSimplex::setPerturbation(int value)
2820{
2821  if(value<=100&&value >=-1000) {
2822    perturbation_=value;
2823  } 
2824}
2825// Sparsity on or off
2826bool 
2827ClpSimplex::sparseFactorization() const
2828{
2829  return factorization_->sparseThreshold()!=0;
2830}
2831void 
2832ClpSimplex::setSparseFactorization(bool value)
2833{
2834  if (value) {
2835    if (!factorization_->sparseThreshold())
2836      factorization_->goSparse();
2837  } else {
2838    factorization_->sparseThreshold(0);
2839  }
2840}
2841void checkCorrect(ClpSimplex * model,int iRow,
2842                  const double * element,const int * rowStart,const int * rowLength,
2843                  const int * column,
2844                  const double * columnLower_, const double * columnUpper_,
2845                  int infiniteUpperC,
2846                  int infiniteLowerC,
2847                  double &maximumUpC,
2848                  double &maximumDownC)
2849{
2850  int infiniteUpper = 0;
2851  int infiniteLower = 0;
2852  double maximumUp = 0.0;
2853  double maximumDown = 0.0;
2854  CoinBigIndex rStart = rowStart[iRow];
2855  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2856  CoinBigIndex j;
2857  double large=1.0e15;
2858  int iColumn;
2859  // Compute possible lower and upper ranges
2860 
2861  for (j = rStart; j < rEnd; ++j) {
2862    double value=element[j];
2863    iColumn = column[j];
2864    if (value > 0.0) {
2865      if (columnUpper_[iColumn] >= large) {
2866        ++infiniteUpper;
2867      } else {
2868        maximumUp += columnUpper_[iColumn] * value;
2869      }
2870      if (columnLower_[iColumn] <= -large) {
2871        ++infiniteLower;
2872      } else {
2873        maximumDown += columnLower_[iColumn] * value;
2874      }
2875    } else if (value<0.0) {
2876      if (columnUpper_[iColumn] >= large) {
2877        ++infiniteLower;
2878      } else {
2879        maximumDown += columnUpper_[iColumn] * value;
2880      }
2881      if (columnLower_[iColumn] <= -large) {
2882        ++infiniteUpper;
2883      } else {
2884        maximumUp += columnLower_[iColumn] * value;
2885      }
2886    }
2887  }
2888  assert (infiniteLowerC==infiniteLower);
2889  assert (infiniteUpperC==infiniteUpper);
2890  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
2891    printf("row %d comp up %g, true up %g\n",iRow,
2892           maximumUpC,maximumUp);
2893  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
2894    printf("row %d comp down %g, true down %g\n",iRow,
2895           maximumDownC,maximumDown);
2896  maximumUpC=maximumUp;
2897  maximumDownC=maximumDown;
2898}
2899
2900/* Tightens primal bounds to make dual faster.  Unless
2901   fixed, bounds are slightly looser than they could be.
2902   This is to make dual go faster and is probably not needed
2903   with a presolve.  Returns non-zero if problem infeasible
2904
2905   Fudge for branch and bound - put bounds on columns of factor *
2906   largest value (at continuous) - should improve stability
2907   in branch and bound on infeasible branches (0.0 is off)
2908*/
2909int 
2910ClpSimplex::tightenPrimalBounds(double factor)
2911{
2912 
2913  // Get a row copy in standard format
2914  CoinPackedMatrix copy;
2915  copy.reverseOrderedCopyOf(*matrix());
2916  // get matrix data pointers
2917  const int * column = copy.getIndices();
2918  const CoinBigIndex * rowStart = copy.getVectorStarts();
2919  const int * rowLength = copy.getVectorLengths(); 
2920  const double * element = copy.getElements();
2921  int numberChanged=1,iPass=0;
2922  double large = largeValue(); // treat bounds > this as infinite
2923#ifndef NDEBUG
2924  double large2= 1.0e10*large;
2925#endif
2926  int numberInfeasible=0;
2927  int totalTightened = 0;
2928
2929  double tolerance = primalTolerance();
2930
2931
2932  // Save column bounds
2933  double * saveLower = new double [numberColumns_];
2934  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
2935  double * saveUpper = new double [numberColumns_];
2936  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
2937
2938  int iRow, iColumn;
2939
2940  // If wanted - tighten column bounds using solution
2941  if (factor) {
2942    assert (factor>1.0);
2943    double largest=0.0;
2944    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2945      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
2946        largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
2947      }
2948    }
2949    largest *= factor;
2950    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2951      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
2952        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
2953        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
2954      }
2955    }
2956  }
2957#define MAXPASS 10
2958
2959  // Loop round seeing if we can tighten bounds
2960  // Would be faster to have a stack of possible rows
2961  // and we put altered rows back on stack
2962  int numberCheck=-1;
2963  while(numberChanged>numberCheck) {
2964
2965    numberChanged = 0; // Bounds tightened this pass
2966   
2967    if (iPass==MAXPASS) break;
2968    iPass++;
2969   
2970    for (iRow = 0; iRow < numberRows_; iRow++) {
2971
2972      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
2973
2974        // possible row
2975        int infiniteUpper = 0;
2976        int infiniteLower = 0;
2977        double maximumUp = 0.0;
2978        double maximumDown = 0.0;
2979        double newBound;
2980        CoinBigIndex rStart = rowStart[iRow];
2981        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2982        CoinBigIndex j;
2983        // Compute possible lower and upper ranges
2984     
2985        for (j = rStart; j < rEnd; ++j) {
2986          double value=element[j];
2987          iColumn = column[j];
2988          if (value > 0.0) {
2989            if (columnUpper_[iColumn] >= large) {
2990              ++infiniteUpper;
2991            } else {
2992              maximumUp += columnUpper_[iColumn] * value;
2993            }
2994            if (columnLower_[iColumn] <= -large) {
2995              ++infiniteLower;
2996            } else {
2997              maximumDown += columnLower_[iColumn] * value;
2998            }
2999          } else if (value<0.0) {
3000            if (columnUpper_[iColumn] >= large) {
3001              ++infiniteLower;
3002            } else {
3003              maximumDown += columnUpper_[iColumn] * value;
3004            }
3005            if (columnLower_[iColumn] <= -large) {
3006              ++infiniteUpper;
3007            } else {
3008              maximumUp += columnLower_[iColumn] * value;
3009            }
3010          }
3011        }
3012        // Build in a margin of error
3013        maximumUp += 1.0e-8*fabs(maximumUp);
3014        maximumDown -= 1.0e-8*fabs(maximumDown);
3015        double maxUp = maximumUp+infiniteUpper*1.0e31;
3016        double maxDown = maximumDown-infiniteLower*1.0e31;
3017        if (maxUp <= rowUpper_[iRow] + tolerance && 
3018            maxDown >= rowLower_[iRow] - tolerance) {
3019         
3020          // Row is redundant - make totally free
3021          // NO - can't do this for postsolve
3022          // rowLower_[iRow]=-COIN_DBL_MAX;
3023          // rowUpper_[iRow]=COIN_DBL_MAX;
3024          //printf("Redundant row in presolveX %d\n",iRow);
3025
3026        } else {
3027          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
3028              maxDown > rowUpper_[iRow]+100.0*tolerance) {
3029            // problem is infeasible - exit at once
3030            numberInfeasible++;
3031            break;
3032          }
3033          double lower = rowLower_[iRow];
3034          double upper = rowUpper_[iRow];
3035          for (j = rStart; j < rEnd; ++j) {
3036            double value=element[j];
3037            iColumn = column[j];
3038            double nowLower = columnLower_[iColumn];
3039            double nowUpper = columnUpper_[iColumn];
3040            if (value > 0.0) {
3041              // positive value
3042              if (lower>-large) {
3043                if (!infiniteUpper) {
3044                  assert(nowUpper < large2);
3045                  newBound = nowUpper + 
3046                    (lower - maximumUp) / value;
3047                  // relax if original was large
3048                  if (fabs(maximumUp)>1.0e8)
3049                    newBound -= 1.0e-12*fabs(maximumUp);
3050                } else if (infiniteUpper==1&&nowUpper>large) {
3051                  newBound = (lower -maximumUp) / value;
3052                  // relax if original was large
3053                  if (fabs(maximumUp)>1.0e8)
3054                    newBound -= 1.0e-12*fabs(maximumUp);
3055                } else {
3056                  newBound = -COIN_DBL_MAX;
3057                }
3058                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3059                  // Tighten the lower bound
3060                  columnLower_[iColumn] = newBound;
3061                  numberChanged++;
3062                  // check infeasible (relaxed)
3063                  if (nowUpper - newBound < 
3064                      -100.0*tolerance) {
3065                    numberInfeasible++;
3066                  }
3067                  // adjust
3068                  double now;
3069                  if (nowLower<-large) {
3070                    now=0.0;
3071                    infiniteLower--;
3072                  } else {
3073                    now = nowLower;
3074                  }
3075                  maximumDown += (newBound-now) * value;
3076                  nowLower = newBound;
3077#ifdef DEBUG
3078                  checkCorrect(this,iRow,
3079                               element, rowStart, rowLength,
3080                               column,
3081                               columnLower_,  columnUpper_,
3082                               infiniteUpper,
3083                               infiniteLower,
3084                               maximumUp,
3085                               maximumDown);
3086#endif
3087                }
3088              } 
3089              if (upper <large) {
3090                if (!infiniteLower) {
3091                  assert(nowLower >- large2);
3092                  newBound = nowLower + 
3093                    (upper - maximumDown) / value;
3094                  // relax if original was large
3095                  if (fabs(maximumDown)>1.0e8)
3096                    newBound += 1.0e-12*fabs(maximumDown);
3097                } else if (infiniteLower==1&&nowLower<-large) {
3098                  newBound =   (upper - maximumDown) / value;
3099                  // relax if original was large
3100                  if (fabs(maximumDown)>1.0e8)
3101                    newBound += 1.0e-12*fabs(maximumDown);
3102                } else {
3103                  newBound = COIN_DBL_MAX;
3104                }
3105                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3106                  // Tighten the upper bound
3107                  columnUpper_[iColumn] = newBound;
3108                  numberChanged++;
3109                  // check infeasible (relaxed)
3110                  if (newBound - nowLower < 
3111                      -100.0*tolerance) {
3112                    numberInfeasible++;
3113                  }
3114                  // adjust
3115                  double now;
3116                  if (nowUpper>large) {
3117                    now=0.0;
3118                    infiniteUpper--;
3119                  } else {
3120                    now = nowUpper;
3121                  }
3122                  maximumUp += (newBound-now) * value;
3123                  nowUpper = newBound;
3124#ifdef DEBUG
3125                  checkCorrect(this,iRow,
3126                               element, rowStart, rowLength,
3127                               column,
3128                               columnLower_,  columnUpper_,
3129                               infiniteUpper,
3130                               infiniteLower,
3131                               maximumUp,
3132                               maximumDown);
3133#endif
3134                }
3135              }
3136            } else {
3137              // negative value
3138              if (lower>-large) {
3139                if (!infiniteUpper) {
3140                  assert(nowLower < large2);
3141                  newBound = nowLower + 
3142                    (lower - maximumUp) / value;
3143                  // relax if original was large
3144                  if (fabs(maximumUp)>1.0e8)
3145                    newBound += 1.0e-12*fabs(maximumUp);
3146                } else if (infiniteUpper==1&&nowLower<-large) {
3147                  newBound = (lower -maximumUp) / value;
3148                  // relax if original was large
3149                  if (fabs(maximumUp)>1.0e8)
3150                    newBound += 1.0e-12*fabs(maximumUp);
3151                } else {
3152                  newBound = COIN_DBL_MAX;
3153                }
3154                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3155                  // Tighten the upper bound
3156                  columnUpper_[iColumn] = newBound;
3157                  numberChanged++;
3158                  // check infeasible (relaxed)
3159                  if (newBound - nowLower < 
3160                      -100.0*tolerance) {
3161                    numberInfeasible++;
3162                  }
3163                  // adjust
3164                  double now;
3165                  if (nowUpper>large) {
3166                    now=0.0;
3167                    infiniteLower--;
3168                  } else {
3169                    now = nowUpper;
3170                  }
3171                  maximumDown += (newBound-now) * value;
3172                  nowUpper = newBound;
3173#ifdef DEBUG
3174                  checkCorrect(this,iRow,
3175                               element, rowStart, rowLength,
3176                               column,
3177                               columnLower_,  columnUpper_,
3178                               infiniteUpper,
3179                               infiniteLower,
3180                               maximumUp,
3181                               maximumDown);
3182#endif
3183                }
3184              }
3185              if (upper <large) {
3186                if (!infiniteLower) {
3187                  assert(nowUpper < large2);
3188                  newBound = nowUpper + 
3189                    (upper - maximumDown) / value;
3190                  // relax if original was large
3191                  if (fabs(maximumDown)>1.0e8)
3192                    newBound -= 1.0e-12*fabs(maximumDown);
3193                } else if (infiniteLower==1&&nowUpper>large) {
3194                  newBound =   (upper - maximumDown) / value;
3195                  // relax if original was large
3196                  if (fabs(maximumDown)>1.0e8)
3197                    newBound -= 1.0e-12*fabs(maximumDown);
3198                } else {
3199                  newBound = -COIN_DBL_MAX;
3200                }
3201                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3202                  // Tighten the lower bound
3203                  columnLower_[iColumn] = newBound;
3204                  numberChanged++;
3205                  // check infeasible (relaxed)
3206                  if (nowUpper - newBound < 
3207                      -100.0*tolerance) {
3208                    numberInfeasible++;
3209                  }
3210                  // adjust
3211                  double now;
3212                  if (nowLower<-large) {
3213                    now=0.0;
3214                    infiniteUpper--;
3215                  } else {
3216                    now = nowLower;
3217                  }
3218                  maximumUp += (newBound-now) * value;
3219                  nowLower = newBound;
3220#ifdef DEBUG
3221                  checkCorrect(this,iRow,
3222                               element, rowStart, rowLength,
3223                               column,
3224                               columnLower_,  columnUpper_,
3225                               infiniteUpper,
3226                               infiniteLower,
3227                               maximumUp,
3228                               maximumDown);
3229#endif
3230                }
3231              }
3232            }
3233          }
3234        }
3235      }
3236    }
3237    totalTightened += numberChanged;
3238    if (iPass==1)
3239      numberCheck=numberChanged>>4;
3240    if (numberInfeasible) break;
3241  }
3242  if (!numberInfeasible) {
3243    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
3244      <<totalTightened
3245      <<CoinMessageEol;
3246    // Set bounds slightly loose
3247    double useTolerance = 1.0e-3;
3248    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3249      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
3250        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) {
3251          // relax enough so will have correct dj
3252#if 1
3253          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3254                                    columnLower_[iColumn]-100.0*useTolerance);
3255          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3256                                    columnUpper_[iColumn]+100.0*useTolerance);
3257#else
3258          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
3259            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
3260              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
3261            } else {
3262              columnLower_[iColumn]=saveLower[iColumn];
3263              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3264                                        saveLower[iColumn]+100.0*useTolerance);
3265            }
3266          } else {
3267            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
3268              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
3269            } else {
3270              columnUpper_[iColumn]=saveUpper[iColumn];
3271              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3272                                        saveUpper[iColumn]-100.0*useTolerance);
3273            }
3274          }
3275#endif
3276        } else {
3277          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
3278            // relax a bit
3279            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
3280                                        saveUpper[iColumn]);
3281          }
3282          if (columnLower_[iColumn]>saveLower[iColumn]) {
3283            // relax a bit
3284            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
3285                                        saveLower[iColumn]);
3286          }
3287        }
3288      }
3289    }
3290  } else {
3291    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
3292      <<numberInfeasible
3293      <<CoinMessageEol;
3294    // restore column bounds
3295    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
3296    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
3297  }
3298  delete [] saveLower;
3299  delete [] saveUpper;
3300  return (numberInfeasible);
3301}
3302// dual
3303#include "ClpSimplexDual.hpp"
3304#include "ClpSimplexPrimal.hpp"
3305int ClpSimplex::dual (int ifValuesPass )
3306{
3307  int saveQuadraticActivated = objective_->activated();
3308  objective_->setActivated(0);
3309  assert (ifValuesPass>=0&&ifValuesPass<3);
3310  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3311      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3312      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3313      additional data and have no destructor or (non-default) constructor.
3314
3315      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3316      in primal and dual.
3317
3318      As far as I can see this is perfectly safe.
3319  */
3320  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass);
3321  if (problemStatus_==10) {
3322    //printf("Cleaning up with primal\n");
3323    int savePerturbation = perturbation_;
3324    perturbation_=100;
3325    bool denseFactorization = initialDenseFactorization();
3326    // It will be safe to allow dense
3327    setInitialDenseFactorization(true);
3328    // Allow for catastrophe
3329    int saveMax = intParam_[ClpMaxNumIteration];
3330    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
3331      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
3332    // check which algorithms allowed
3333    int dummy;
3334    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
3335      returnCode = ((ClpSimplexPrimal *) this)->primal(1);
3336    else
3337      returnCode = ((ClpSimplexDual *) this)->dual(0);
3338    if (problemStatus_==3&&numberIterations_<saveMax) {
3339      printf("looks like trouble - too many iterations in clean up - trying again\n");
3340      // flatten solution and try again
3341      int iRow,iColumn;
3342      for (iRow=0;iRow<numberRows_;iRow++) {
3343        if (getRowStatus(iRow)!=basic) {
3344          setRowStatus(iRow,superBasic);
3345          // but put to bound if close
3346          if (fabs(rowActivity_[iRow]-rowLower_[iRow])
3347              <=primalTolerance_) {
3348            rowActivity_[iRow]=rowLower_[iRow];
3349            setRowStatus(iRow,atLowerBound);
3350          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])
3351                     <=primalTolerance_) {
3352            rowActivity_[iRow]=rowUpper_[iRow];
3353            setRowStatus(iRow,atUpperBound);
3354          }
3355        }
3356      }
3357      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3358        if (getColumnStatus(iColumn)!=basic) {
3359          setColumnStatus(iColumn,superBasic);
3360          // but put to bound if close
3361          if (fabs(columnActivity_[iColumn]-columnLower_[iColumn])
3362              <=primalTolerance_) {
3363            columnActivity_[iColumn]=columnLower_[iColumn];
3364            setColumnStatus(iColumn,atLowerBound);
3365          } else if (fabs(columnActivity_[iColumn]
3366                          -columnUpper_[iColumn])
3367                     <=primalTolerance_) {
3368            columnActivity_[iColumn]=columnUpper_[iColumn];
3369            setColumnStatus(iColumn,atUpperBound);
3370          }
3371        }
3372      }
3373      problemStatus_=-1;
3374      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 
3375                                          2*numberRows_+numberColumns_,saveMax);
3376      perturbation_=savePerturbation;
3377      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
3378      if (problemStatus_==3&&numberIterations_<saveMax) 
3379        printf("looks like real trouble - too many iterations in second clean up - giving up\n");
3380    }
3381    intParam_[ClpMaxNumIteration] = saveMax;
3382
3383    setInitialDenseFactorization(denseFactorization);
3384    perturbation_=savePerturbation;
3385    if (problemStatus_==10) 
3386      problemStatus_=0;
3387  }
3388  objective_->setActivated(saveQuadraticActivated);
3389  return returnCode;
3390}
3391#include "ClpQuadraticObjective.hpp"
3392// primal
3393int ClpSimplex::primal (int ifValuesPass )
3394{
3395  // See if nonlinear
3396  if (objective_->type()>1&&objective_->activated()) 
3397    return reducedGradient();
3398 
3399  assert (ifValuesPass>=0&&ifValuesPass<3);
3400  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3401      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3402      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3403      additional data and have no destructor or (non-default) constructor.
3404
3405      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3406      in primal and dual.
3407
3408      As far as I can see this is perfectly safe.
3409  */
3410  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
3411  if (problemStatus_==10) {
3412    //printf("Cleaning up with dual\n");
3413    int savePerturbation = perturbation_;
3414    perturbation_=100;
3415    bool denseFactorization = initialDenseFactorization();
3416    // It will be safe to allow dense
3417    setInitialDenseFactorization(true);
3418    // check which algorithms allowed
3419    int dummy;
3420    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3421      double saveBound = dualBound_;
3422      // upperOut_ has largest away from bound
3423      dualBound_=CoinMin(2.0*upperOut_,dualBound_);
3424      returnCode = ((ClpSimplexDual *) this)->dual(0);
3425      dualBound_=saveBound;
3426    } else {
3427      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
3428    }
3429    setInitialDenseFactorization(denseFactorization);
3430    perturbation_=savePerturbation;
3431    if (problemStatus_==10) 
3432      problemStatus_=0;
3433  }
3434  return returnCode;
3435}
3436/* Dual ranging.
3437   This computes increase/decrease in cost for each given variable and corresponding
3438   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3439   and numberColumns.. for artificials/slacks.
3440   For non-basic variables the sequence number will be that of the non-basic variables.
3441   
3442   Up to user to provide correct length arrays.
3443   
3444   Returns non-zero if infeasible unbounded etc
3445*/
3446#include "ClpSimplexOther.hpp"
3447int ClpSimplex::dualRanging(int numberCheck,const int * which,
3448                            double * costIncrease, int * sequenceIncrease,
3449                            double * costDecrease, int * sequenceDecrease)
3450{
3451  int savePerturbation = perturbation_;
3452  perturbation_=100;
3453  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3454  if (problemStatus_==10) {
3455    //printf("Cleaning up with dual\n");
3456    bool denseFactorization = initialDenseFactorization();
3457    // It will be safe to allow dense
3458    setInitialDenseFactorization(true);
3459    // check which algorithms allowed
3460    int dummy;
3461    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3462      // upperOut_ has largest away from bound
3463      double saveBound=dualBound_;
3464      if (upperOut_>0.0)
3465        dualBound_=2.0*upperOut_;
3466      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3467      dualBound_=saveBound;
3468    } else {
3469      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3470    }
3471    setInitialDenseFactorization(denseFactorization);
3472    if (problemStatus_==10) 
3473      problemStatus_=0;
3474  }
3475  perturbation_=savePerturbation;
3476  if (problemStatus_) {
3477    finish(); // get rid of arrays
3478    return 1; // odd status
3479  }
3480  ((ClpSimplexOther *) this)->dualRanging(numberCheck,which,
3481                                          costIncrease,sequenceIncrease,
3482                                          costDecrease,sequenceDecrease);
3483  finish(); // get rid of arrays
3484  return 0;
3485}
3486/* Primal ranging.
3487   This computes increase/decrease in value for each given variable and corresponding
3488   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3489   and numberColumns.. for artificials/slacks.
3490   For basic variables the sequence number will be that of the basic variables.
3491   
3492   Up to user to providen correct length arrays.
3493   
3494   Returns non-zero if infeasible unbounded etc
3495*/
3496int ClpSimplex::primalRanging(int numberCheck,const int * which,
3497                  double * valueIncrease, int * sequenceIncrease,
3498                  double * valueDecrease, int * sequenceDecrease)
3499{
3500  int savePerturbation = perturbation_;
3501  perturbation_=100;
3502  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3503  if (problemStatus_==10) {
3504    //printf("Cleaning up with dual\n");
3505    bool denseFactorization = initialDenseFactorization();
3506    // It will be safe to allow dense
3507    setInitialDenseFactorization(true);
3508    // check which algorithms allowed
3509    int dummy;
3510    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3511      // upperOut_ has largest away from bound
3512      double saveBound=dualBound_;
3513      if (upperOut_>0.0)
3514        dualBound_=2.0*upperOut_;
3515      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3516      dualBound_=saveBound;
3517    } else {
3518      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3519    }
3520    setInitialDenseFactorization(denseFactorization);
3521    if (problemStatus_==10) 
3522      problemStatus_=0;
3523  }
3524  perturbation_=savePerturbation;
3525  if (problemStatus_) {
3526    finish(); // get rid of arrays
3527    return 1; // odd status
3528  }
3529  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
3530                                          valueIncrease,sequenceIncrease,
3531                                          valueDecrease,sequenceDecrease);
3532  finish(); // get rid of arrays
3533  return 0;
3534}
3535#include "ClpSimplexNonlinear.hpp"
3536/* Solves nonlinear problem using SLP - may be used as crash
3537   for other algorithms when number of iterations small
3538*/
3539int 
3540ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
3541{
3542  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
3543}
3544// Solves non-linear using reduced gradient
3545int ClpSimplex::reducedGradient(int phase)
3546{
3547  if (objective_->type()<2||!objective_->activated()) {
3548    // no quadratic part
3549    return primal(0);
3550  }
3551  // get feasible
3552  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
3553    objective_->setActivated(0);
3554    double saveDirection = optimizationDirection();
3555    setOptimizationDirection(0.0);
3556    primal(1);
3557    setOptimizationDirection(saveDirection);
3558    objective_->setActivated(1);
3559    // still infeasible
3560    if (numberPrimalInfeasibilities())
3561      return 0;
3562  }
3563  // Now enter method
3564  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
3565  return returnCode;
3566}
3567#include "ClpPredictorCorrector.hpp"
3568#include "ClpCholeskyBase.hpp"
3569// Preference is WSSMP, TAUCS, UFL (just ordering) then base
3570#if WSSMP_BARRIER
3571#include "ClpCholeskyWssmp.hpp"
3572#include "ClpCholeskyWssmpKKT.hpp"
3573#elif UFL_BARRIER
3574#include "ClpCholeskyUfl.hpp"
3575#elif TAUCS_BARRIER
3576#include "ClpCholeskyTaucs.hpp"
3577#endif
3578#include "ClpPresolve.hpp"
3579/* Solves using barrier (assumes you have good cholesky factor code).
3580   Does crossover to simplex if asked*/
3581int 
3582ClpSimplex::barrier(bool crossover)
3583{
3584  ClpSimplex * model2 = this;
3585  int savePerturbation=perturbation_;
3586  ClpInterior barrier;
3587  barrier.borrowModel(*model2);
3588// Preference is WSSMP, UFL, TAUCS then base
3589#ifdef WSSMP_BARRIER
3590  ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
3591  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
3592  //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
3593#elif UFL_BARRIER
3594  ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
3595#elif TAUCS_BARRIER
3596  ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
3597#else
3598  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
3599#endif
3600  barrier.setCholesky(cholesky);
3601  int numberRows = model2->numberRows();
3602  int numberColumns = model2->numberColumns();
3603  int saveMaxIts = model2->maximumIterations();
3604  if (saveMaxIts<1000) {
3605    barrier.setMaximumBarrierIterations(saveMaxIts);
3606    model2->setMaximumIterations(1000000);
3607  }
3608  barrier.primalDual();
3609  int barrierStatus=barrier.status();
3610  double gap = barrier.complementarityGap();
3611  // get which variables are fixed
3612  double * saveLower=NULL;
3613  double * saveUpper=NULL;
3614  ClpPresolve pinfo2;
3615  ClpSimplex * saveModel2=NULL;
3616  int numberFixed = barrier.numberFixed();
3617  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
3618    // may as well do presolve
3619    int numberRows = barrier.numberRows();
3620    int numberColumns = barrier.numberColumns();
3621    int numberTotal = numberRows+numberColumns;
3622    saveLower = new double [numberTotal];
3623    saveUpper = new double [numberTotal];
3624    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
3625    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
3626    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
3627    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
3628    barrier.fixFixed();
3629    saveModel2=model2;
3630  }
3631  barrier.returnModel(*model2);
3632  double * rowPrimal = new double [numberRows];
3633  double * columnPrimal = new double [numberColumns];
3634  double * rowDual = new double [numberRows];
3635  double * columnDual = new double [numberColumns];
3636  // move solutions other way
3637  CoinMemcpyN(model2->primalRowSolution(),
3638              numberRows,rowPrimal);
3639  CoinMemcpyN(model2->dualRowSolution(),
3640              numberRows,rowDual);
3641  CoinMemcpyN(model2->primalColumnSolution(),
3642              numberColumns,columnPrimal);
3643  CoinMemcpyN(model2->dualColumnSolution(),
3644              numberColumns,columnDual);
3645  if (saveModel2) {
3646    // do presolve
3647    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
3648                                   false,5,true);
3649  }
3650  if (barrierStatus<4&&crossover) {
3651    // make sure no status left
3652    model2->createStatus();
3653    // solve
3654    model2->setPerturbation(100);
3655    // throw some into basis
3656    {
3657      int numberRows = model2->numberRows();
3658      int numberColumns = model2->numberColumns();
3659      double * dsort = new double[numberColumns];
3660      int * sort = new int[numberColumns];
3661      int n=0;
3662      const double * columnLower = model2->columnLower();
3663      const double * columnUpper = model2->columnUpper();
3664      const double * primalSolution = model2->primalColumnSolution();
3665      double tolerance = 10.0*primalTolerance_;
3666      int i;
3667      for ( i=0;i<numberRows;i++) 
3668        model2->setRowStatus(i,superBasic);
3669      for ( i=0;i<numberColumns;i++) {
3670        double distance = CoinMin(columnUpper[i]-primalSolution[i],
3671                              primalSolution[i]-columnLower[i]);
3672        if (distance>tolerance) {
3673          dsort[n]=-distance;
3674          sort[n++]=i;
3675          model2->setStatus(i,superBasic);
3676        } else if (distance>primalTolerance_) {
3677          model2->setStatus(i,superBasic);
3678        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
3679          model2->setStatus(i,atLowerBound);
3680        } else {
3681          model2->setStatus(i,atUpperBound);
3682        }
3683      }
3684      CoinSort_2(dsort,dsort+n,sort);
3685      n = CoinMin(numberRows,n);
3686      for ( i=0;i<n;i++) {
3687        int iColumn = sort[i];
3688        model2->setStatus(iColumn,basic);
3689      }
3690      delete [] sort;
3691      delete [] dsort;
3692    }
3693    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
3694      int numberRows = model2->numberRows();
3695      int numberColumns = model2->numberColumns();
3696      // just primal values pass
3697      double saveScale = model2->objectiveScale();
3698      model2->setObjectiveScale(1.0e-3);
3699      model2->primal(2);
3700      model2->setObjectiveScale(saveScale);
3701      // save primal solution and copy back dual
3702      CoinMemcpyN(model2->primalRowSolution(),
3703                  numberRows,rowPrimal);
3704      CoinMemcpyN(rowDual,
3705                  numberRows,model2->dualRowSolution());
3706      CoinMemcpyN(model2->primalColumnSolution(),
3707                  numberColumns,columnPrimal);
3708      CoinMemcpyN(columnDual,
3709                  numberColumns,model2->dualColumnSolution());
3710      //model2->primal(1);
3711      // clean up reduced costs and flag variables
3712      {
3713        double * dj = model2->dualColumnSolution();
3714        double * cost = model2->objective();
3715        double * saveCost = new double[numberColumns];
3716        memcpy(saveCost,cost,numberColumns*sizeof(double));
3717        double * saveLower = new double[numberColumns];
3718        double * lower = model2->columnLower();
3719        memcpy(saveLower,lower,numberColumns*sizeof(double));
3720        double * saveUpper = new double[numberColumns];
3721        double * upper = model2->columnUpper();
3722        memcpy(saveUpper,upper,numberColumns*sizeof(double));
3723        int i;
3724        double tolerance = 10.0*dualTolerance_;
3725        for ( i=0;i<numberColumns;i++) {
3726          if (model2->getStatus(i)==basic) {
3727            dj[i]=0.0;
3728          } else if (model2->getStatus(i)==atLowerBound) {
3729            if (optimizationDirection_*dj[i]<tolerance) {
3730              if (optimizationDirection_*dj[i]<0.0) {
3731                //if (dj[i]<-1.0e-3)
3732                //printf("bad dj at lb %d %g\n",i,dj[i]);
3733                cost[i] -= dj[i];
3734                dj[i]=0.0;
3735              }
3736            } else {
3737              upper[i]=lower[i];
3738            }
3739          } else if (model2->getStatus(i)==atUpperBound) {
3740            if (optimizationDirection_*dj[i]>tolerance) {
3741              if (optimizationDirection_*dj[i]>0.0) {
3742                //if (dj[i]>1.0e-3)
3743                //printf("bad dj at ub %d %g\n",i,dj[i]);
3744                cost[i] -= dj[i];
3745                dj[i]=0.0;
3746              }
3747            } else {
3748              lower[i]=upper[i];
3749            }
3750          }
3751        }
3752        // just dual values pass
3753        //model2->setLogLevel(63);
3754        //model2->setFactorizationFrequency(1);
3755        model2->dual(2);
3756        memcpy(cost,saveCost,numberColumns*sizeof(double));
3757        delete [] saveCost;
3758        memcpy(lower,saveLower,numberColumns*sizeof(double));
3759        delete [] saveLower;
3760        memcpy(upper,saveUpper,numberColumns*sizeof(double));
3761        delete [] saveUpper;
3762      }
3763      // and finish
3764      // move solutions
3765      CoinMemcpyN(rowPrimal,
3766                  numberRows,model2->primalRowSolution());
3767      CoinMemcpyN(columnPrimal,
3768                  numberColumns,model2->primalColumnSolution());
3769    }
3770//     double saveScale = model2->objectiveScale();
3771//     model2->setObjectiveScale(1.0e-3);
3772//     model2->primal(2);
3773//    model2->setObjectiveScale(saveScale);
3774    model2->primal(1);
3775  } else if (barrierStatus==4&&crossover) {
3776    // memory problems
3777    model2->setPerturbation(savePerturbation);
3778    model2->createStatus();
3779    model2->dual();
3780  }
3781  model2->setMaximumIterations(saveMaxIts);
3782  delete [] rowPrimal;
3783  delete [] columnPrimal;
3784  delete [] rowDual;
3785  delete [] columnDual;
3786  if (saveLower) {
3787    pinfo2.postsolve(true);
3788    delete model2;
3789    model2=saveModel2;
3790    int numberRows = model2->numberRows();
3791    int numberColumns = model2->numberColumns();
3792    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
3793    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
3794    delete [] saveLower;
3795    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
3796    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
3797    delete [] saveUpper;
3798    model2->primal(1);
3799  }
3800  model2->setPerturbation(savePerturbation);
3801  return model2->status();
3802}
3803/* For strong branching.  On input lower and upper are new bounds
3804   while on output they are objective function values (>1.0e50 infeasible).
3805   Return code is 0 if nothing interesting, -1 if infeasible both
3806   ways and +1 if infeasible one way (check values to see which one(s))
3807*/
3808int ClpSimplex::strongBranching(int numberVariables,const int * variables,
3809                                double * newLower, double * newUpper,
3810                                double ** outputSolution,
3811                                int * outputStatus, int * outputIterations,
3812                                bool stopOnFirstInfeasible,
3813                                bool alwaysFinish)
3814{
3815  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
3816                                                    newLower,  newUpper,outputSolution,
3817                                                    outputStatus, outputIterations,
3818                                                    stopOnFirstInfeasible,
3819                                                    alwaysFinish);
3820}
3821/* Borrow model.  This is so we dont have to copy large amounts
3822   of data around.  It assumes a derived class wants to overwrite
3823   an empty model with a real one - while it does an algorithm.
3824   This is same as ClpModel one, but sets scaling on etc. */
3825void 
3826ClpSimplex::borrowModel(ClpModel & otherModel) 
3827{
3828  ClpModel::borrowModel(otherModel);
3829  createStatus();
3830  //ClpDualRowSteepest steep1;
3831  //setDualRowPivotAlgorithm(steep1);
3832  //ClpPrimalColumnSteepest steep2;
3833  //setPrimalColumnPivotAlgorithm(steep2);
3834}
3835void 
3836ClpSimplex::borrowModel(ClpSimplex & otherModel) 
3837{
3838  ClpModel::borrowModel(otherModel);
3839  createStatus();
3840  dualBound_ = otherModel.dualBound_;
3841  dualTolerance_ = otherModel.dualTolerance_;
3842  primalTolerance_ = otherModel.primalTolerance_;
3843  delete dualRowPivot_;
3844  dualRowPivot_ = otherModel.dualRowPivot_->clone(true);
3845  delete primalColumnPivot_;
3846  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(true);
3847  perturbation_ = otherModel.perturbation_;
3848  specialOptions_ = otherModel.specialOptions_;
3849  automaticScale_ = otherModel.automaticScale_;
3850}
3851typedef struct {
3852  double optimizationDirection;
3853  double dblParam[ClpLastDblParam];
3854  double objectiveValue;
3855  double dualBound;
3856  double dualTolerance;
3857  double primalTolerance;
3858  double sumDualInfeasibilities;
3859  double sumPrimalInfeasibilities;
3860  double infeasibilityCost;
3861  int numberRows;
3862  int numberColumns;
3863  int intParam[ClpLastIntParam];
3864  int numberIterations;
3865  int problemStatus;
3866  int maximumIterations;
3867  int lengthNames;
3868  int numberDualInfeasibilities;
3869  int numberDualInfeasibilitiesWithoutFree;
3870  int numberPrimalInfeasibilities;
3871  int numberRefinements;
3872  int scalingFlag;
3873  int algorithm;
3874  unsigned int specialOptions;
3875  int dualPivotChoice;
3876  int primalPivotChoice;
3877  int matrixStorageChoice;
3878} Clp_scalars;
3879
3880int outDoubleArray(double * array, int length, FILE * fp)
3881{
3882  int numberWritten;
3883  if (array&&length) {
3884    numberWritten = fwrite(&length,sizeof(int),1,fp);
3885    if (numberWritten!=1)
3886      return 1;
3887    numberWritten = fwrite(array,sizeof(double),length,fp);
3888    if (numberWritten!=length)
3889      return 1;
3890  } else {
3891    length = 0;
3892    numberWritten = fwrite(&length,sizeof(int),1,fp);
3893    if (numberWritten!=1)
3894      return 1;
3895  }
3896  return 0;
3897}
3898// Save model to file, returns 0 if success
3899int
3900ClpSimplex::saveModel(const char * fileName)
3901{
3902  FILE * fp = fopen(fileName,"wb");
3903  if (fp) {
3904    Clp_scalars scalars;
3905    int i;
3906    CoinBigIndex numberWritten;
3907    // Fill in scalars
3908    scalars.optimizationDirection = optimizationDirection_;
3909    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
3910    scalars.objectiveValue = objectiveValue_;
3911    scalars.dualBound = dualBound_;
3912    scalars.dualTolerance = dualTolerance_;
3913    scalars.primalTolerance = primalTolerance_;
3914    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
3915    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
3916    scalars.infeasibilityCost = infeasibilityCost_;
3917    scalars.numberRows = numberRows_;
3918    scalars.numberColumns = numberColumns_;
3919    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
3920    scalars.numberIterations = numberIterations_;
3921    scalars.problemStatus = problemStatus_;
3922    scalars.maximumIterations = maximumIterations();
3923    scalars.lengthNames = lengthNames_;
3924    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
3925    scalars.numberDualInfeasibilitiesWithoutFree
3926      = numberDualInfeasibilitiesWithoutFree_;
3927    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
3928    scalars.numberRefinements = numberRefinements_;
3929    scalars.scalingFlag = scalingFlag_;
3930    scalars.algorithm = algorithm_;
3931    scalars.specialOptions = specialOptions_;
3932    scalars.dualPivotChoice = dualRowPivot_->type();
3933    scalars.primalPivotChoice = primalColumnPivot_->type();
3934    scalars.matrixStorageChoice = matrix_->type();
3935
3936    // put out scalars
3937    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
3938    if (numberWritten!=1)
3939      return 1;
3940    // strings
3941    CoinBigIndex length;
3942    for (i=0;i<ClpLastStrParam;i++) {
3943      length = strParam_[i].size();
3944      numberWritten = fwrite(&length,sizeof(int),1,fp);
3945      if (numberWritten!=1)
3946        return 1;
3947      if (length) {
3948        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
3949        if (numberWritten!=1)
3950          return 1;
3951      }
3952    }
3953    // arrays - in no particular order
3954    if (outDoubleArray(rowActivity_,numberRows_,fp))
3955        return 1;
3956    if (outDoubleArray(columnActivity_,numberColumns_,fp))
3957        return 1;
3958    if (outDoubleArray(dual_,numberRows_,fp))
3959        return 1;
3960    if (outDoubleArray(reducedCost_,numberColumns_,fp))
3961        return 1;
3962    if (outDoubleArray(rowLower_,numberRows_,fp))
3963        return 1;
3964    if (outDoubleArray(rowUpper_,numberRows_,fp))
3965        return 1;
3966    if (outDoubleArray(objective(),numberColumns_,fp))
3967        return 1;
3968    if (outDoubleArray(rowObjective_,numberRows_,fp))
3969        return 1;
3970    if (outDoubleArray(columnLower_,numberColumns_,fp))
3971        return 1;
3972    if (outDoubleArray(columnUpper_,numberColumns_,fp))
3973        return 1;
3974    if (ray_) {
3975      if (problemStatus_==1)
3976        if (outDoubleArray(ray_,numberRows_,fp))
3977          return 1;
3978      else if (problemStatus_==2)
3979        if (outDoubleArray(ray_,numberColumns_,fp))
3980          return 1;
3981      else
3982        if (outDoubleArray(NULL,0,fp))
3983          return 1;
3984    } else {
3985      if (outDoubleArray(NULL,0,fp))
3986        return 1;
3987    }
3988    if (status_&&(numberRows_+numberColumns_)>0) {
3989      length = numberRows_+numberColumns_;
3990      numberWritten = fwrite(&length,sizeof(int),1,fp);
3991      if (numberWritten!=1)
3992        return 1;
3993      numberWritten = fwrite(status_,sizeof(char),length, fp);
3994      if (numberWritten!=length)
3995        return 1;
3996    } else {
3997      length = 0;
3998      numberWritten = fwrite(&length,sizeof(int),1,fp);
3999      if (numberWritten!=1)
4000        return 1;
4001    }
4002    if (lengthNames_) {
4003      char * array = 
4004        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4005      char * put = array;
4006      assert (numberRows_ == (int) rowNames_.size());
4007      for (i=0;i<numberRows_;i++) {
4008        assert((int)rowNames_[i].size()<=lengthNames_);
4009        strcpy(put,rowNames_[i].c_str());
4010        put += lengthNames_+1;
4011      }
4012      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
4013      if (numberWritten!=numberRows_)
4014        return 1;
4015      put=array;
4016      assert (numberColumns_ == (int) columnNames_.size());
4017      for (i=0;i<numberColumns_;i++) {
4018        assert((int) columnNames_[i].size()<=lengthNames_);
4019        strcpy(put,columnNames_[i].c_str());
4020        put += lengthNames_+1;
4021      }
4022      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
4023      if (numberWritten!=numberColumns_)
4024        return 1;
4025      delete [] array;
4026    }
4027    // just standard type at present
4028    assert (matrix_->type()==1);
4029    assert (matrix_->getNumCols() == numberColumns_);
4030    assert (matrix_->getNumRows() == numberRows_);
4031    // we are going to save with gaps
4032    length = matrix_->getVectorStarts()[numberColumns_-1]
4033      + matrix_->getVectorLengths()[numberColumns_-1];
4034    numberWritten = fwrite(&length,sizeof(int),1,fp);
4035    if (numberWritten!=1)
4036      return 1;
4037    numberWritten = fwrite(matrix_->getElements(),
4038                           sizeof(double),length,fp);
4039    if (numberWritten!=length)
4040      return 1;
4041    numberWritten = fwrite(matrix_->getIndices(),
4042                           sizeof(int),length,fp);
4043    if (numberWritten!=length)
4044      return 1;
4045    numberWritten = fwrite(matrix_->getVectorStarts(),
4046                           sizeof(int),numberColumns_+1,fp);
4047    if (numberWritten!=numberColumns_+1)
4048      return 1;
4049    numberWritten = fwrite(matrix_->getVectorLengths(),
4050                           sizeof(int),numberColumns_,fp);
4051    if (numberWritten!=numberColumns_)
4052      return 1;
4053    // finished
4054    fclose(fp);
4055    return 0;
4056  } else {
4057    return -1;
4058  }
4059}
4060
4061int inDoubleArray(double * &array, int length, FILE * fp)
4062{
4063  int numberRead;
4064  int length2;
4065  numberRead = fread(&length2,sizeof(int),1,fp);
4066  if (numberRead!=1)
4067    return 1;
4068  if (length2) {
4069    // lengths must match
4070    if (length!=length2)
4071      return 2;
4072    array = new double[length];
4073    numberRead = fread(array,sizeof(double),length,fp);
4074    if (numberRead!=length)
4075      return 1;
4076  } 
4077  return 0;
4078}
4079/* Restore model from file, returns 0 if success,
4080   deletes current model */
4081int 
4082ClpSimplex::restoreModel(const char * fileName)
4083{
4084  FILE * fp = fopen(fileName,"rb");
4085  if (fp) {
4086    // Get rid of current model
4087    // save event handler in case already set
4088    ClpEventHandler * handler = eventHandler_->clone();
4089    ClpModel::gutsOfDelete();
4090    eventHandler_ = handler;
4091    gutsOfDelete(0);
4092    int i;
4093    for (i=0;i<6;i++) {
4094      rowArray_[i]=NULL;
4095      columnArray_[i]=NULL;
4096    }
4097    // get an empty factorization so we can set tolerances etc
4098    factorization_ = new ClpFactorization();
4099    // Say sparse
4100    factorization_->sparseThreshold(1);
4101    Clp_scalars scalars;
4102    CoinBigIndex numberRead;
4103
4104    // get scalars
4105    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
4106    if (numberRead!=1)
4107      return 1;
4108    // Fill in scalars
4109    optimizationDirection_ = scalars.optimizationDirection;
4110    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
4111    objectiveValue_ = scalars.objectiveValue;
4112    dualBound_ = scalars.dualBound;
4113    dualTolerance_ = scalars.dualTolerance;
4114    primalTolerance_ = scalars.primalTolerance;
4115    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
4116    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
4117    infeasibilityCost_ = scalars.infeasibilityCost;
4118    numberRows_ = scalars.numberRows;
4119    numberColumns_ = scalars.numberColumns;
4120    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
4121    numberIterations_ = scalars.numberIterations;
4122    problemStatus_ = scalars.problemStatus;
4123    setMaximumIterations(scalars.maximumIterations);
4124    lengthNames_ = scalars.lengthNames;
4125    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
4126    numberDualInfeasibilitiesWithoutFree_
4127      = scalars.numberDualInfeasibilitiesWithoutFree;
4128    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
4129    numberRefinements_ = scalars.numberRefinements;
4130    scalingFlag_ = scalars.scalingFlag;
4131    algorithm_ = scalars.algorithm;
4132    specialOptions_ = scalars.specialOptions;
4133    // strings
4134    CoinBigIndex length;
4135    for (i=0;i<ClpLastStrParam;i++) {
4136      numberRead = fread(&length,sizeof(int),1,fp);
4137      if (numberRead!=1)
4138        return 1;
4139      if (length) {
4140        char * array = new char[length+1];
4141        numberRead = fread(array,length,1,fp);
4142        if (numberRead!=1)
4143          return 1;
4144        array[length]='\0';
4145        strParam_[i]=array;
4146        delete [] array;
4147      }
4148    }
4149    // arrays - in no particular order
4150    if (inDoubleArray(rowActivity_,numberRows_,fp))
4151        return 1;
4152    if (inDoubleArray(columnActivity_,numberColumns_,fp))
4153        return 1;
4154    if (inDoubleArray(dual_,numberRows_,fp))
4155        return 1;
4156    if (inDoubleArray(reducedCost_,numberColumns_,fp))
4157        return 1;
4158    if (inDoubleArray(rowLower_,numberRows_,fp))
4159        return 1;
4160    if (inDoubleArray(rowUpper_,numberRows_,fp))
4161        return 1;
4162    double * objective;
4163    if (inDoubleArray(objective,numberColumns_,fp))
4164        return 1;
4165    delete objective_;
4166    objective_ = new ClpLinearObjective(objective,numberColumns_);
4167    delete [] objective;
4168    if (inDoubleArray(rowObjective_,numberRows_,fp))
4169        return 1;
4170    if (inDoubleArray(columnLower_,numberColumns_,fp))
4171        return 1;
4172    if (inDoubleArray(columnUpper_,numberColumns_,fp))
4173        return 1;
4174    if (problemStatus_==1) {
4175      if (inDoubleArray(ray_,numberRows_,fp))
4176        return 1;
4177    } else if (problemStatus_==2) {
4178      if (inDoubleArray(ray_,numberColumns_,fp))
4179        return 1;
4180    } else {
4181      // ray should be null
4182      numberRead = fread(&length,sizeof(int),1,fp);
4183      if (numberRead!=1)
4184        return 1;
4185      if (length)
4186        return 2;
4187    }
4188    delete [] status_;
4189    status_=NULL;
4190    // status region
4191    numberRead = fread(&length,sizeof(int),1,fp);
4192    if (numberRead!=1)
4193        return 1;
4194    if (length) {
4195      if (length!=numberRows_+numberColumns_)
4196        return 1;
4197      status_ = new char unsigned[length];
4198      numberRead = fread(status_,sizeof(char),length, fp);
4199      if (numberRead!=length)
4200        return 1;
4201    }
4202    if (lengthNames_) {
4203      char * array = 
4204        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4205      char * get = array;
4206      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
4207      if (numberRead!=numberRows_)
4208        return 1;
4209      rowNames_ = std::vector<std::string> ();
4210      rowNames_.resize(numberRows_);
4211      for (i=0;i<numberRows_;i++) {
4212        rowNames_.push_back(get);
4213        get += lengthNames_+1;
4214      }
4215      get = array;
4216      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
4217      if (numberRead!=numberColumns_)
4218        return 1;
4219      columnNames_ = std::vector<std::string> ();
4220      columnNames_.resize(numberColumns_);
4221      for (i=0;i<numberColumns_;i++) {
4222        columnNames_.push_back(get);
4223        get += lengthNames_+1;
4224      }
4225      delete [] array;
4226    }
4227    // Pivot choices
4228    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
4229    delete dualRowPivot_;
4230    switch ((scalars.dualPivotChoice&63)) {
4231    default:
4232      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
4233    case 1:
4234      // Dantzig
4235      dualRowPivot_ = new ClpDualRowDantzig();
4236      break;
4237    case 2:
4238      // Steepest - use mode
4239      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
4240      break;
4241    }
4242    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
4243    delete primalColumnPivot_;
4244    switch ((scalars.primalPivotChoice&63)) {
4245    default:
4246      printf("Need another primalPivot case %d\n",
4247             scalars.primalPivotChoice&63);
4248    case 1:
4249      // Dantzig
4250      primalColumnPivot_ = new ClpPrimalColumnDantzig();
4251      break;
4252    case 2:
4253      // Steepest - use mode
4254      primalColumnPivot_
4255        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
4256      break;
4257    }
4258    assert(scalars.matrixStorageChoice==1);
4259    delete matrix_;
4260    // get arrays
4261    numberRead = fread(&length,sizeof(int),1,fp);
4262    if (numberRead!=1)
4263      return 1;
4264    double * elements = new double[length];
4265    int * indices = new int[length];
4266    CoinBigIndex * starts = new CoinBigIndex[numberColumns_+1];
4267    int * lengths = new int[numberColumns_];
4268    numberRead = fread(elements, sizeof(double),length,fp);
4269    if (numberRead!=length)
4270      return 1;
4271    numberRead = fread(indices, sizeof(int),length,fp);
4272    if (numberRead!=length)
4273      return 1;
4274    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
4275    if (numberRead!=numberColumns_+1)
4276      return 1;
4277    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
4278    if (numberRead!=numberColumns_)
4279      return 1;
4280    // assign matrix
4281    CoinPackedMatrix * matrix = new CoinPackedMatrix();
4282    matrix->assignMatrix(true, numberRows_, numberColumns_,
4283                         length, elements, indices, starts, lengths);
4284    // and transfer to Clp
4285    matrix_ = new ClpPackedMatrix(matrix);
4286    // finished
4287    fclose(fp);
4288    return 0;
4289  } else {
4290    return -1;
4291  }
4292  return 0;
4293}
4294// value of incoming variable (in Dual)
4295double 
4296ClpSimplex::valueIncomingDual() const
4297{
4298  // Need value of incoming for list of infeasibilities as may be infeasible
4299  double valueIncoming = (dualOut_/alpha_)*directionOut_;
4300  if (directionIn_==-1)
4301    valueIncoming = upperIn_-valueIncoming;
4302  else
4303    valueIncoming = lowerIn_-valueIncoming;
4304  return valueIncoming;
4305}
4306// Sanity check on input data - returns true if okay
4307bool 
4308ClpSimplex::sanityCheck()
4309{
4310  // bad if empty
4311  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
4312    int infeasNumber[2];
4313    double infeasSum[2];
4314    problemStatus_=emptyProblem(infeasNumber,infeasSum,false);
4315    numberDualInfeasibilities_=infeasNumber[0];
4316    sumDualInfeasibilities_=infeasSum[0];
4317    numberPrimalInfeasibilities_=infeasNumber[1];
4318    sumPrimalInfeasibilities_=infeasSum[1];
4319    return false;
4320  }
4321  int numberBad ;
4322  double largestBound, smallestBound, minimumGap;
4323  double smallestObj, largestObj;
4324  int firstBad;
4325  int modifiedBounds=0;
4326  int i;
4327  numberBad=0;
4328  firstBad=-1;
4329  minimumGap=1.0e100;
4330  smallestBound=1.0e100;
4331  largestBound=0.0;
4332  smallestObj=1.0e100;
4333  largestObj=0.0;
4334  // If bounds are too close - fix
4335  double fixTolerance = 10.0*primalTolerance_;
4336  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
4337    double value;
4338    value = fabs(cost_[i]);
4339    if (value>1.0e50) {
4340      numberBad++;
4341      if (firstBad<0)
4342        firstBad=i;
4343    } else if (value) {
4344      if (value>largestObj)
4345        largestObj=value;
4346      if (value<smallestObj)
4347        smallestObj=value;
4348    }
4349    value=upper_[i]-lower_[i];
4350    if (value<-primalTolerance_) {
4351      numberBad++;
4352      if (firstBad<0)
4353        firstBad=i;
4354    } else if (value<=fixTolerance) {
4355      if (value) {
4356        // modify
4357        upper_[i] = lower_[i];
4358        modifiedBounds++;
4359      }
4360    } else {
4361      if (value<minimumGap)
4362        minimumGap=value;
4363    }
4364    if (lower_[i]>-1.0e100&&lower_[i]) {
4365      value = fabs(lower_[i]);
4366      if (value>largestBound)
4367        largestBound=value;
4368      if (value<smallestBound)
4369        smallestBound=value;
4370    }
4371    if (upper_[i]<1.0e100&&upper_[i]) {
4372      value = fabs(upper_[i]);
4373      if (value>largestBound)
4374        largestBound=value;
4375      if (value<smallestBound)
4376        smallestBound=value;
4377    }
4378  }
4379  if (largestBound)
4380    handler_->message(CLP_RIMSTATISTICS3,messages_)
4381      <<smallestBound
4382      <<largestBound
4383      <<minimumGap
4384      <<CoinMessageEol;
4385  minimumGap=1.0e100;
4386  smallestBound=1.0e100;
4387  largestBound=0.0;
4388  for (i=0;i<numberColumns_;i++) {
4389    double value;
4390    value = fabs(cost_[i]);
4391    if (value>1.0e50) {
4392      numberBad++;
4393      if (firstBad<0)
4394        firstBad=i;
4395    } else if (value) {
4396      if (value>largestObj)
4397        largestObj=value;
4398      if (value<smallestObj)
4399        smallestObj=value;
4400    }
4401    value=upper_[i]-lower_[i];
4402    if (value<-primalTolerance_) {
4403      numberBad++;
4404      if (firstBad<0)
4405        firstBad=i;
4406    } else if (value<=fixTolerance) {
4407      if (value) {
4408        // modify
4409        upper_[i] = lower_[i];
4410        modifiedBounds++;
4411      }
4412    } else {
4413      if (value<minimumGap)
4414        minimumGap=value;
4415    }
4416    if (lower_[i]>-1.0e100&&lower_[i]) {
4417      value = fabs(lower_[i]);
4418      if (value>largestBound)
4419        largestBound=value;
4420      if (value<smallestBound)
4421        smallestBound=value;
4422    }
4423    if (upper_[i]<1.0e100&&upper_[i]) {
4424      value = fabs(upper_[i]);
4425      if (value>largestBound)
4426        largestBound=value;
4427      if (value<smallestBound)
4428        smallestBound=value;
4429    }
4430  }
4431  char rowcol[]={'R','C'};
4432  if (numberBad) {
4433    handler_->message(CLP_BAD_BOUNDS,messages_)
4434      <<numberBad
4435      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
4436      <<CoinMessageEol;
4437    problemStatus_=4;
4438    return false;
4439  }
4440  if (modifiedBounds)
4441    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
4442      <<modifiedBounds
4443      <<CoinMessageEol;
4444  handler_->message(CLP_RIMSTATISTICS1,messages_)
4445    <<smallestObj
4446    <<largestObj
4447    <<CoinMessageEol;  if (largestBound)
4448    handler_->message(CLP_RIMSTATISTICS2,messages_)
4449      <<smallestBound
4450      <<largestBound
4451      <<minimumGap
4452      <<CoinMessageEol;
4453  return true;
4454}
4455// Set up status array (for OsiClp)
4456void 
4457ClpSimplex::createStatus() 
4458{
4459  if(!status_)
4460    status_ = new unsigned char [numberColumns_+numberRows_];
4461  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
4462  int i;
4463  // set column status to one nearest zero
4464  for (i=0;i<numberColumns_;i++) {
4465#if 0
4466    if (columnLower_[i]>=0.0) {
4467      setColumnStatus(i,atLowerBound);
4468    } else if (columnUpper_[i]<=0.0) {
4469      setColumnStatus(i,atUpperBound);
4470    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4471      // free
4472      setColumnStatus(i,isFree);
4473    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4474      setColumnStatus(i,atLowerBound);
4475    } else {
4476      setColumnStatus(i,atUpperBound);
4477    }
4478#else
4479    setColumnStatus(i,atLowerBound);
4480#endif
4481  }
4482  for (i=0;i<numberRows_;i++) {
4483    setRowStatus(i,basic);
4484  }
4485}
4486/* Sets up all slack basis and resets solution to
4487   as it was after initial load or readMps */
4488void ClpSimplex::allSlackBasis(bool resetSolution)
4489{
4490  createStatus();
4491  if (resetSolution) {
4492    // put back to as it was originally
4493    int i;
4494    // set column status to one nearest zero
4495    // But set value to zero if lb <0.0 and ub>0.0
4496    for (i=0;i<numberColumns_;i++) {
4497      if (columnLower_[i]>=0.0) {
4498        columnActivity_[i]=columnLower_[i];
4499        setColumnStatus(i,atLowerBound);
4500      } else if (columnUpper_[i]<=0.0) {
4501        columnActivity_[i]=columnUpper_[i];
4502        setColumnStatus(i,atUpperBound);
4503      } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4504        // free
4505        columnActivity_[i]=0.0;
4506        setColumnStatus(i,isFree);
4507      } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4508        columnActivity_[i]=0.0;
4509        setColumnStatus(i,atLowerBound);
4510      } else {
4511        columnActivity_[i]=0.0;
4512        setColumnStatus(i,atUpperBound);
4513      }
4514    }
4515  }
4516}
4517/* Loads a problem (the constraints on the
4518   rows are given by lower and upper bounds). If a pointer is 0 then the
4519   following values are the default:
4520   <ul>
4521   <li> <code>colub</code>: all columns have upper bound infinity
4522   <li> <code>collb</code>: all columns have lower bound 0
4523   <li> <code>rowub</code>: all rows have upper bound infinity
4524   <li> <code>rowlb</code>: all rows have lower bound -infinity
4525   <li> <code>obj</code>: all variables have 0 objective coefficient
4526   </ul>
4527*/
4528void 
4529ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
4530                    const double* collb, const double* colub,   
4531                    const double* obj,
4532                    const double* rowlb, const double* rowub,
4533                    const double * rowObjective)
4534{
4535  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4536                        rowObjective);
4537  createStatus();
4538}
4539void 
4540ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
4541                    const double* collb, const double* colub,   
4542                    const double* obj,
4543                    const double* rowlb, const double* rowub,
4544                    const double * rowObjective)
4545{
4546  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4547                        rowObjective);
4548  createStatus();
4549}
4550
4551/* Just like the other loadProblem() method except that the matrix is
4552   given in a standard column major ordered format (without gaps). */
4553void 
4554ClpSimplex::loadProblem (  const int numcols, const int numrows,
4555                    const CoinBigIndex* start, const int* index,
4556                    const double* value,
4557                    const double* collb, const double* colub,   
4558                    const double* obj,
4559                    const double* rowlb, const double* rowub,
4560                    const double * rowObjective)
4561{
4562  ClpModel::loadProblem(numcols, numrows, start, index, value,
4563                          collb, colub, obj, rowlb, rowub,
4564                          rowObjective);
4565  createStatus();
4566}
4567void 
4568ClpSimplex::loadProblem (  const int numcols, const int numrows,
4569                           const CoinBigIndex* start, const int* index,
4570                           const double* value,const int * length,
4571                           const double* collb, const double* colub,   
4572                           const double* obj,
4573                           const double* rowlb, const double* rowub,
4574                           const double * rowObjective)
4575{
4576  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
4577                          collb, colub, obj, rowlb, rowub,
4578                          rowObjective);
4579  createStatus();
4580}
4581// Read an mps file from the given filename
4582int 
4583ClpSimplex::readMps(const char *filename,
4584            bool keepNames,
4585            bool ignoreErrors)
4586{
4587  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
4588  createStatus();
4589  return status;
4590}
4591// Just check solution (for external use)
4592void 
4593ClpSimplex::checkSolution()
4594{
4595  // put in standard form
4596  createRim(7+8+16);
4597  dualTolerance_=dblParam_[ClpDualTolerance];
4598  primalTolerance_=dblParam_[ClpPrimalTolerance];
4599  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
4600  checkDualSolution();
4601  if (!numberDualInfeasibilities_&&
4602      !numberPrimalInfeasibilities_)
4603    problemStatus_=0;
4604  else
4605    problemStatus_=-1;
4606#ifdef CLP_DEBUG
4607  int i;
4608  double value=0.0;
4609  for (i=0;i<numberRows_+numberColumns_;i++)
4610    value += dj_[i]*solution_[i];
4611  printf("dual value %g, primal %g\n",value,objectiveValue());
4612#endif
4613  // release extra memory
4614  deleteRim(0);
4615}
4616/* Crash - at present just aimed at dual, returns
4617   -2 if dual preferred and crash basis created
4618   -1 if dual preferred and all slack basis preferred
4619   0 if basis going in was not all slack
4620   1 if primal preferred and all slack basis preferred
4621   2 if primal preferred and crash basis created.
4622   
4623   if gap between bounds <="gap" variables can be flipped
4624   
4625   If "pivot" is
4626   0 No pivoting (so will just be choice of algorithm)
4627   1 Simple pivoting e.g. gub
4628   2 Mini iterations
4629*/
4630int 
4631ClpSimplex::crash(double gap,int pivot)
4632{
4633  assert(!rowObjective_); // not coded
4634  int iColumn;
4635  int numberBad=0;
4636  int numberBasic=0;
4637  double dualTolerance=dblParam_[ClpDualTolerance];
4638  //double primalTolerance=dblParam_[ClpPrimalTolerance];
4639  int returnCode=0;
4640  // If no basis then make all slack one
4641  if (!status_)
4642    createStatus();
4643 
4644  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4645    if (getColumnStatus(iColumn)==basic)
4646      numberBasic++;
4647  }
4648  if (!numberBasic) {
4649    // all slack
4650    double * dj = new double [numberColumns_];
4651    double * solution = columnActivity_;
4652    const double * linearObjective = objective();
4653    //double objectiveValue=0.0;
4654    int iColumn;
4655    double direction = optimizationDirection_;
4656    // direction is actually scale out not scale in
4657    if (direction)
4658      direction = 1.0/direction;
4659    for (iColumn=0;iColumn<numberColumns_;iColumn++)
4660      dj[iColumn] = direction*linearObjective[iColumn];
4661    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4662      // assume natural place is closest to zero
4663      double lowerBound = columnLower_[iColumn];
4664      double upperBound = columnUpper_[iColumn];
4665      if (lowerBound>-1.0e20||upperBound<1.0e20) {
4666        bool atLower;
4667        if (fabs(upperBound)<fabs(lowerBound)) {
4668          atLower=false;
4669          setColumnStatus(iColumn,atUpperBound);
4670          solution[iColumn]=upperBound;
4671        } else {
4672          atLower=true;
4673          setColumnStatus(iColumn,atLowerBound);
4674          solution[iColumn]=lowerBound;
4675        }
4676        if (dj[iColumn]<0.0) {
4677          // should be at upper bound
4678          if (atLower) {
4679            // can we flip
4680            if (upperBound-lowerBound<=gap) {
4681              columnActivity_[iColumn]=upperBound;
4682              setColumnStatus(iColumn,atUpperBound);
4683            } else if (dj[iColumn]<-dualTolerance) {
4684              numberBad++;
4685            }
4686          }
4687        } else if (dj[iColumn]>0.0) {
4688          // should be at lower bound
4689          if (!atLower) {
4690            // can we flip
4691            if (upperBound-lowerBound<=gap) {
4692              columnActivity_[iColumn]=lowerBound;
4693              setColumnStatus(iColumn,atLowerBound);
4694            } else if (dj[iColumn]>dualTolerance) {
4695              numberBad++;
4696            }
4697          }
4698        }
4699      } else {
4700        // free
4701        setColumnStatus(iColumn,isFree);
4702        if (fabs(dj[iColumn])>dualTolerance) 
4703          numberBad++;
4704      }
4705    }
4706    if (numberBad||pivot) {
4707      if (!pivot) {
4708        delete [] dj;
4709        returnCode = 1;
4710      } else {
4711        // see if can be made dual feasible with gubs etc
4712        double * pi = new double[numberRows_];
4713        memset (pi,0,numberRows_*sizeof(double));
4714        int * way = new int[numberColumns_];
4715        int numberIn = 0;
4716
4717        // Get column copy
4718        CoinPackedMatrix * columnCopy = matrix();
4719        // Get a row copy in standard format
4720        CoinPackedMatrix copy;
4721        copy.reverseOrderedCopyOf(*columnCopy);
4722        // get matrix data pointers
4723        const int * column = copy.getIndices();
4724        const CoinBigIndex * rowStart = copy.getVectorStarts();
4725        const int * rowLength = copy.getVectorLengths(); 
4726        const double * elementByRow = copy.getElements();
4727        //const int * row = columnCopy->getIndices();
4728        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4729        //const int * columnLength = columnCopy->getVectorLengths();
4730        //const double * element = columnCopy->getElements();
4731
4732
4733        // if equality row and bounds mean artificial in basis bad
4734        // then do anyway
4735
4736        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4737          // - if we want to reduce dj, + if we want to increase
4738          int thisWay = 100;
4739          double lowerBound = columnLower_[iColumn];
4740          double upperBound = columnUpper_[iColumn];
4741          if (upperBound>lowerBound) {
4742            switch(getColumnStatus(iColumn)) {
4743             
4744            case basic:
4745              thisWay=0;
4746            case ClpSimplex::isFixed:
4747              break;
4748            case isFree:
4749            case superBasic:
4750              if (dj[iColumn]<-dualTolerance) 
4751                thisWay = 1;
4752              else if (dj[iColumn]>dualTolerance) 
4753                thisWay = -1;
4754              else
4755                thisWay =0;
4756              break;
4757            case atUpperBound:
4758              if (dj[iColumn]>dualTolerance) 
4759                thisWay = -1;
4760              else if (dj[iColumn]<-dualTolerance) 
4761                thisWay = -3;
4762              else
4763                thisWay = -2;
4764              break;
4765            case atLowerBound:
4766              if (dj[iColumn]<-dualTolerance) 
4767                thisWay = 1;
4768              else if (dj[iColumn]>dualTolerance) 
4769                thisWay = 3;
4770              else
4771                thisWay = 2;
4772              break;
4773            }
4774          }
4775          way[iColumn] = thisWay;
4776        }
4777        /*if (!numberBad)
4778          printf("Was dual feasible before passes - rows %d\n",
4779          numberRows_);*/
4780        int lastNumberIn = -100000;
4781        int numberPasses=5;
4782        while (numberIn>lastNumberIn+numberRows_/100) {
4783          lastNumberIn = numberIn;
4784          // we need to maximize chance of doing good
4785          int iRow;
4786          for (iRow=0;iRow<numberRows_;iRow++) {
4787            double lowerBound = rowLower_[iRow];
4788            double upperBound = rowUpper_[iRow];
4789            if (getRowStatus(iRow)==basic) {
4790              // see if we can find a column to pivot on
4791              int j;
4792              // down is amount pi can go down
4793              double maximumDown = COIN_DBL_MAX;
4794              double maximumUp = COIN_DBL_MAX;
4795              double minimumDown =0.0;
4796              double minimumUp =0.0;
4797              int iUp=-1;
4798              int iDown=-1;
4799              int iUpB=-1;
4800              int iDownB=-1;
4801              if (lowerBound<-1.0e20)
4802                maximumUp = -1.0;
4803              if (upperBound>1.0e20)
4804                maximumDown = -1.0;
4805              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4806                int iColumn = column[j];
4807                double value = elementByRow[j];
4808                double djValue = dj[iColumn];
4809                /* way -
4810                   -3 - okay at upper bound with negative dj
4811                   -2 - marginal at upper bound with zero dj - can only decrease
4812                   -1 - bad at upper bound
4813                   0 - we can never pivot on this row
4814                   1 - bad at lower bound
4815                   2 - marginal at lower bound with zero dj - can only increase
4816                   3 - okay at lower bound with positive dj
4817                   100 - fine we can just ignore
4818                */
4819                if (way[iColumn]!=100) {
4820                  switch(way[iColumn]) {
4821                   
4822                  case -3:
4823                    if (value>0.0) {
4824                      if (maximumDown*value>-djValue) {
4825                        maximumDown = -djValue/value;
4826                        iDown = iColumn;
4827                      }
4828                    } else {
4829                      if (-maximumUp*value>-djValue) {
4830                        maximumUp = djValue/value;
4831                        iUp = iColumn;
4832                      }
4833                    }
4834                    break;
4835                  case -2:
4836                    if (value>0.0) {
4837                      maximumDown = 0.0;
4838                    } else {
4839                      maximumUp = 0.0;
4840                    }
4841                    break;
4842                  case -1:
4843                    // see if could be satisfied
4844                    // dj value > 0
4845                    if (value>0.0) {
4846                      maximumDown=0.0;
4847                      if (maximumUp*value<djValue-dualTolerance) {
4848                        maximumUp = 0.0; // would improve but not enough
4849                      } else {
4850                        if (minimumUp*value<djValue) {
4851                          minimumUp = djValue/value;
4852                          iUpB = iColumn;
4853                        }
4854                      }
4855                    } else {
4856                      maximumUp=0.0;
4857                      if (-maximumDown*value<djValue-dualTolerance) {
4858                        maximumDown = 0.0; // would improve but not enough
4859                      } else {
4860                        if (-minimumDown*value<djValue) {
4861                          minimumDown = -djValue/value;
4862                          iDownB = iColumn;
4863                        }
4864                      }
4865                    }
4866                   
4867                    break;
4868                  case 0:
4869                    maximumDown = -1.0;
4870                    maximumUp=-1.0;
4871                    break;
4872                  case 1:
4873                    // see if could be satisfied
4874                    // dj value < 0
4875                    if (value>0.0) {
4876                      maximumUp=0.0;
4877                      if (maximumDown*value<-djValue-dualTolerance) {
4878                        maximumDown = 0.0; // would improve but not enough
4879                      } else {
4880                        if (minimumDown*value<-djValue) {
4881                          minimumDown = -djValue/value;
4882                          iDownB = iColumn;
4883                        }
4884                      }
4885                    } else {
4886                      maximumDown=0.0;
4887                      if (-maximumUp*value<-djValue-dualTolerance) {
4888                        maximumUp = 0.0; // would improve but not enough
4889                      } else {
4890                        if (-minimumUp*value<-djValue) {
4891                          minimumUp = djValue/value;
4892                          iUpB = iColumn;
4893                        }
4894                      }
4895                    }
4896                   
4897                    break;
4898                  case 2:
4899                    if (value>0.0) {
4900                      maximumUp = 0.0;
4901                    } else {
4902                      maximumDown = 0.0;
4903                    }
4904                   
4905                    break;
4906                  case 3:
4907                    if (value>0.0) {
4908                      if (maximumUp*value>djValue) {
4909                        maximumUp = djValue/value;
4910                        iUp = iColumn;
4911                      }
4912                    } else {
4913                      if (-maximumDown*value>djValue) {
4914                        maximumDown = -djValue/value;
4915                        iDown = iColumn;
4916                      }
4917                    }
4918                   
4919                    break;
4920                  default:
4921                    break;
4922                  }
4923                }
4924              }
4925              if (iUpB>=0)
4926                iUp=iUpB;
4927              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
4928                iUp=-1;
4929              if (iDownB>=0)
4930                iDown=iDownB;
4931              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
4932                iDown=-1;
4933              if (iUp>=0||iDown>=0) {
4934                // do something
4935                if (iUp>=0&&iDown>=0) {
4936                  if (maximumDown>maximumUp)
4937                    iUp=-1;
4938                }
4939                double change;
4940                int kColumn;
4941                if (iUp>=0) {
4942                  kColumn=iUp;
4943                  change=maximumUp;
4944                  // just do minimum if was dual infeasible
4945                  // ? only if maximum large?
4946                  if (minimumUp>0.0)
4947                    change=minimumUp;
4948                  setRowStatus(iRow,atUpperBound);
4949                } else {
4950                  kColumn=iDown;
4951                  change=-maximumDown;
4952                  // just do minimum if was dual infeasible
4953                  // ? only if maximum large?
4954                  if (minimumDown>0.0)
4955                    change=-minimumDown;
4956                  setRowStatus(iRow,atLowerBound);
4957                }
4958                assert (fabs(change)<1.0e20);
4959                setColumnStatus(kColumn,basic);
4960                numberIn++;
4961                pi[iRow]=change;
4962                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4963                  int iColumn = column[j];
4964                  double value = elementByRow[j];
4965                  double djValue = dj[iColumn]-change*value;
4966                  dj[iColumn]=djValue;
4967                  if (abs(way[iColumn])==1) {
4968                    numberBad--;
4969                    /*if (!numberBad)
4970                      printf("Became dual feasible at row %d out of %d\n",
4971                      iRow, numberRows_);*/
4972                    lastNumberIn=-1000000;
4973                  }
4974                  int thisWay = 100;
4975                  double lowerBound = columnLower_[iColumn];
4976                  double upperBound = columnUpper_[iColumn];
4977                  if (upperBound>lowerBound) {
4978                    switch(getColumnStatus(iColumn)) {
4979                     
4980                    case basic:
4981                      thisWay=0;
4982                    case isFixed:
4983                      break;
4984                    case isFree:
4985                    case superBasic:
4986                      if (djValue<-dualTolerance) 
4987                        thisWay = 1;
4988                      else if (djValue>dualTolerance) 
4989                        thisWay = -1;
4990                      else
4991                        { thisWay =0; abort();}
4992                      break;
4993                    case atUpperBound:
4994                      if (djValue>dualTolerance) 
4995                        { thisWay =-1; abort();}
4996                      else if (djValue<-dualTolerance) 
4997                        thisWay = -3;
4998                      else
4999                        thisWay = -2;
5000                      break;
5001                    case atLowerBound:
5002                      if (djValue<-dualTolerance) 
5003                        { thisWay =1; abort();}
5004                      else if (djValue>dualTolerance) 
5005                        thisWay = 3;
5006                      else
5007                        thisWay = 2;
5008                      break;
5009                    }
5010                  }
5011                  way[iColumn] = thisWay;
5012                }
5013              }
5014            }
5015          }
5016          if (numberIn==lastNumberIn||numberBad||pivot<2)
5017            break;
5018          if (!(--numberPasses))
5019            break;
5020          //printf("%d put in so far\n",numberIn);
5021        }
5022        // last attempt to flip
5023        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5024          double lowerBound = columnLower_[iColumn];
5025          double upperBound = columnUpper_[iColumn];
5026          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
5027            double djValue=dj[iColumn];
5028            switch(getColumnStatus(iColumn)) {
5029             
5030            case basic:
5031            case ClpSimplex::isFixed:
5032              break;
5033            case isFree:
5034            case superBasic:
5035              break;
5036            case atUpperBound:
5037              if (djValue>dualTolerance) {
5038                setColumnStatus(iColumn,atUpperBound);
5039                solution[iColumn]=upperBound;
5040              } 
5041              break;
5042            case atLowerBound:
5043              if (djValue<-dualTolerance) {
5044                setColumnStatus(iColumn,atUpperBound);
5045                solution[iColumn]=upperBound;
5046              }
5047              break;
5048            }
5049          }
5050        }
5051        delete [] pi;
5052        delete [] dj;
5053        delete [] way;
5054        handler_->message(CLP_CRASH,messages_)
5055          <<numberIn
5056          <<numberBad
5057          <<CoinMessageEol;
5058        returnCode =  -1;
5059      }
5060    } else {
5061      delete [] dj;
5062      returnCode =  -1;
5063    }
5064    //cleanStatus();
5065  }
5066  return returnCode;
5067}
5068/* Pivot in a variable and out a variable.  Returns 0 if okay,
5069   1 if inaccuracy forced re-factorization, -1 if would be singular.
5070   Also updates primal/dual infeasibilities.
5071   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
5072*/
5073int ClpSimplex::pivot()
5074{
5075  // scaling not allowed
5076  assert (!scalingFlag_);
5077  // assume In_ and Out_ are correct and directionOut_ set
5078  // (or In_ if flip
5079  lowerIn_ = lower_[sequenceIn_];
5080  valueIn_ = solution_[sequenceIn_];
5081  upperIn_ = upper_[sequenceIn_];
5082  dualIn_ = dj_[sequenceIn_];
5083  lowerOut_ = lower_[sequenceOut_];
5084  valueOut_ = solution_[sequenceOut_];
5085  upperOut_ = upper_[sequenceOut_];
5086  // for now assume primal is feasible (or in dual)
5087  dualOut_ = dj_[sequenceOut_];
5088  assert(fabs(dualOut_)<1.0e-6);
5089  bool roundAgain = true;
5090  int returnCode=0;
5091  while (roundAgain) {
5092    roundAgain=false;
5093    unpack(rowArray_[1]);
5094    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
5095    alpha_=0.0;
5096    int i;
5097    int * index = rowArray_[1]->getIndices();
5098    int number = rowArray_[1]->getNumElements();
5099    double * element = rowArray_[1]->denseVector();
5100    for (i=0;i<number;i++) {
5101      int ii = index[i];
5102      if ( pivotVariable_[ii]==sequenceOut_) {
5103        pivotRow_=ii;
5104        alpha_=element[pivotRow_];
5105        break;
5106      }
5107    }
5108    assert (fabs(alpha_)>1.0e-12);
5109    // we are going to subtract movement from current basic
5110    double movement;
5111    // see where incoming will go to
5112    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
5113      // flip so go to bound
5114      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
5115    } else {
5116      // get where outgoing needs to get to
5117      double outValue = (directionOut_<0) ? upperOut_ : lowerOut_;
5118      // solutionOut_ - movement*alpha_ == outValue
5119      movement = (valueOut_-outValue)/alpha_;
5120      // set directionIn_ correctly
5121      directionIn_ = (movement>0) ? 1 :-1;
5122    }
5123    // update primal solution
5124    for (i=0;i<number;i++) {
5125      int ii = index[i];
5126      // get column
5127      int ij = pivotVariable_[ii];
5128      solution_[ij] -= movement*element[ii];
5129      element[ii]=0.0;
5130    }
5131    rowArray_[1]->setNumElements(0);
5132    // see where something went to
5133    if (sequenceOut_<0) {
5134      if (directionIn_<0) {
5135        assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
5136        solution_[sequenceIn_]=upperIn_;
5137      } else {
5138        assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
5139        solution_[sequenceIn_]=lowerIn_;
5140      }
5141    } else {
5142      if (directionOut_<0) {
5143        assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
5144        solution_[sequenceOut_]=upperOut_;
5145      } else {
5146        assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
5147        solution_[sequenceOut_]=lowerOut_;
5148      }
5149      solution_[sequenceIn_]=valueIn_+movement;
5150    }
5151    double objectiveChange = dualIn_*movement;
5152    // update duals
5153    if (pivotRow_>=0) {
5154      assert (fabs(alpha_)>1.0e-8);
5155      double multiplier = dualIn_/alpha_;
5156      rowArray_[0]->insert(pivotRow_,multiplier);
5157      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
5158      // put row of tableau in rowArray[0] and columnArray[0]
5159      matrix_->transposeTimes(this,-1.0,
5160                              rowArray_[0],columnArray_[1],columnArray_[0]);
5161      // update column djs
5162      int i;
5163      int * index = columnArray_[0]->getIndices();
5164      int number = columnArray_[0]->getNumElements();
5165      double * element = columnArray_[0]->denseVector();
5166      for (i=0;i<number;i++) {
5167        int ii = index[i];
5168        dj_[ii] += element[ii];
5169        element[ii]=0.0;
5170      }
5171      columnArray_[0]->setNumElements(0);
5172      // and row djs
5173      index = rowArray_[0]->getIndices();
5174      number = rowArray_[0]->getNumElements();
5175      element = rowArray_[0]->denseVector();
5176      for (i=0;i<number;i++) {
5177        int ii = index[i];
5178        dj_[ii+numberColumns_] += element[ii];
5179        dual_[ii] = dj_[ii+numberColumns_];
5180        element[ii]=0.0;
5181      }
5182      rowArray_[0]->setNumElements(0);
5183      // check incoming
5184      assert (fabs(dj_[sequenceIn_])<1.0e-6);
5185    }
5186   
5187    // if stable replace in basis
5188    int updateStatus = factorization_->replaceColumn(this,
5189                                                     rowArray_[2],
5190                                                     rowArray_[1],
5191                                                     pivotRow_,
5192                                                     alpha_);
5193    bool takePivot=true;
5194    // if no pivots, bad update but reasonable alpha - take and invert
5195    if (updateStatus==2&&
5196        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
5197      updateStatus=4;
5198    if (updateStatus==1||updateStatus==4) {
5199      // slight error
5200      if (factorization_->pivots()>5||updateStatus==4) {
5201        returnCode=-1;
5202      }
5203    } else if (updateStatus==2) {
5204      // major error
5205      rowArray_[1]->clear();
5206      takePivot=false;
5207      if (factorization_->pivots()) {
5208        // refactorize here
5209        statusOfProblem();
5210        roundAgain=true;
5211      } else {
5212        returnCode=1;
5213      }
5214    } else if (updateStatus==3) {
5215      // out of memory
5216      // increase space if not many iterations
5217      if (factorization_->pivots()<
5218          0.5*factorization_->maximumPivots()&&
5219          factorization_->pivots()<200)
5220        factorization_->areaFactor(
5221                                   factorization_->areaFactor() * 1.1);
5222      returnCode =-1; // factorize now
5223    }
5224    if (takePivot) {
5225      int save = algorithm_;
5226      // make simple so always primal
5227      algorithm_=1;
5228      housekeeping(objectiveChange);
5229      algorithm_=save;
5230    }
5231  }
5232  if (returnCode == -1) {
5233    // refactorize here
5234    statusOfProblem();
5235  } else {
5236    // just for now - recompute anyway
5237    gutsOfSolution(NULL,NULL);
5238  }
5239  return returnCode;
5240}
5241
5242/* Pivot in a variable and choose an outgoing one.  Assumes primal
5243   feasible - will not go through a bound.  Returns step length in theta
5244   Returns ray in ray_ (or NULL if no pivot)
5245   Return codes as before but -1 means no acceptable pivot
5246*/
5247int ClpSimplex::primalPivotResult()
5248{
5249  assert (sequenceIn_>=0);
5250  valueIn_=solution_[sequenceIn_];
5251  lowerIn_=lower_[sequenceIn_];
5252  upperIn_=upper_[sequenceIn_];
5253  dualIn_=dj_[sequenceIn_];
5254
5255  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
5256  if (returnCode<0&&returnCode>-4) {
5257    return 0;
5258  } else {
5259    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
5260           returnCode);
5261    return -1;
5262  }
5263}
5264 
5265/* Pivot out a variable and choose an incoing one.  Assumes dual
5266   feasible - will not go through a reduced cost. 
5267   Returns step length in theta
5268   Returns ray in ray_ (or NULL if no pivot)
5269   Return codes as before but -1 means no acceptable pivot
5270*/
5271int 
5272ClpSimplex::dualPivotResult()
5273{
5274  return ((ClpSimplexDual *) this)->pivotResult();
5275}
5276// Factorization frequency
5277int 
5278ClpSimplex::factorizationFrequency() const
5279{
5280  if (factorization_)
5281    return factorization_->maximumPivots();
5282  else 
5283    return -1;
5284}
5285void 
5286ClpSimplex::setFactorizationFrequency(int value)
5287{
5288  if (factorization_)
5289    factorization_->maximumPivots(value);
5290}
5291// Common bits of coding for dual and primal
5292int 
5293ClpSimplex::startup(int ifValuesPass)
5294{
5295  // sanity check
5296  // bad if empty (trap here to avoid using bad matrix_)
5297  if (!matrix_||(!matrix_->getNumElements()&&objective_->type()<2)) {
5298    int infeasNumber[2];
5299    double infeasSum[2];
5300    problemStatus_=emptyProblem(infeasNumber,infeasSum);
5301    numberDualInfeasibilities_=infeasNumber[0];
5302    sumDualInfeasibilities_=infeasSum[0];
5303    numberPrimalInfeasibilities_=infeasNumber[1];
5304    sumPrimalInfeasibilities_=infeasSum[1];
5305    return 2;
5306  }
5307  pivotRow_=-1;
5308  sequenceIn_=-1;
5309  sequenceOut_=-1;
5310  secondaryStatus_=0;
5311
5312  primalTolerance_=dblParam_[ClpPrimalTolerance];
5313  dualTolerance_=dblParam_[ClpDualTolerance];
5314  if (problemStatus_!=10)
5315    numberIterations_=0;
5316
5317  // put in standard form (and make row copy)
5318  // create modifiable copies of model rim and do optional scaling
5319  bool goodMatrix=createRim(7+8+16,true);
5320
5321  if (goodMatrix) {
5322    // Model looks okay
5323    // Do initial factorization
5324    // and set certain stuff
5325    // We can either set increasing rows so ...IsBasic gives pivot row
5326    // or we can just increment iBasic one by one
5327    // for now let ...iBasic give pivot row
5328    factorization_->increasingRows(2);
5329    // row activities have negative sign
5330    factorization_->slackValue(-1.0);
5331    factorization_->zeroTolerance(1.0e-13);
5332    // Switch off dense (unless special option set)
5333    int saveThreshold = factorization_->denseThreshold();
5334    factorization_->setDenseThreshold(0);
5335    // If values pass then perturb (otherwise may be optimal so leave a bit)
5336    if (ifValuesPass) {
5337      // do perturbation if asked for
5338     
5339      if (perturbation_<100) {
5340        if (algorithm_>0&&(objective_->type()<2||!objective_->activated())) {
5341          ((ClpSimplexPrimal *) this)->perturb(0);
5342        } else if (algorithm_<0) {
5343        ((ClpSimplexDual *) this)->perturb();
5344        }
5345      }
5346    }
5347    // for primal we will change bounds using infeasibilityCost_
5348    if (nonLinearCost_==NULL&&algorithm_>0) {
5349      // get a valid nonlinear cost function
5350      delete nonLinearCost_;
5351      nonLinearCost_= new ClpNonLinearCost(this);
5352    }
5353   
5354    // loop round to clean up solution if values pass
5355    int numberThrownOut = -1;
5356    int totalNumberThrownOut=0;
5357    while(numberThrownOut) {
5358      int status = internalFactorize(ifValuesPass ? 10 : 0);
5359      if (status<0)
5360        return 1; // some error
5361      else
5362        numberThrownOut = status;
5363     
5364      // for this we need clean basis so it is after factorize
5365      if (!numberThrownOut)
5366        numberThrownOut = gutsOfSolution(  NULL,NULL,
5367                                         ifValuesPass!=0);
5368      else
5369        matrix_->rhsOffset(this,true); // redo rhs offset
5370      totalNumberThrownOut+= numberThrownOut;
5371     
5372    }
5373   
5374    if (totalNumberThrownOut)
5375      handler_->message(CLP_SINGULARITIES,messages_)
5376        <<totalNumberThrownOut
5377        <<CoinMessageEol;
5378    // Switch back dense
5379    factorization_->setDenseThreshold(saveThreshold);
5380   
5381    problemStatus_ = -1;
5382   
5383    // number of times we have declared optimality
5384    numberTimesOptimal_=0;
5385
5386    return 0;
5387  } else {
5388    // bad matrix
5389    return 2;
5390  }
5391   
5392}
5393
5394
5395void 
5396ClpSimplex::finish()
5397{
5398  // Get rid of some arrays and empty factorization
5399  deleteRim();
5400  // Skip message if changing algorithms
5401  if (problemStatus_!=10) {
5402    if (problemStatus_==-1)
5403      problemStatus_=4;
5404    assert(problemStatus_>=0&&problemStatus_<6);
5405    handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
5406      <<objectiveValue()
5407      <<CoinMessageEol;
5408  }
5409  factorization_->relaxAccuracyCheck(1.0);
5410  // get rid of any network stuff - could do more
5411  factorization_->cleanUp();
5412}
5413// Save data
5414ClpDataSave
5415ClpSimplex::saveData() 
5416{
5417  ClpDataSave saved;
5418  saved.dualBound_ = dualBound_;
5419  saved.infeasibilityCost_ = infeasibilityCost_;
5420  saved.sparseThreshold_ = factorization_->sparseThreshold();
5421  saved.perturbation_ = perturbation_;
5422  // Progress indicator
5423  delete progress_;
5424  progress_ = new ClpSimplexProgress (this);
5425  return saved;
5426}
5427// Restore data
5428void 
5429ClpSimplex::restoreData(ClpDataSave saved)
5430{
5431  factorization_->sparseThreshold(saved.sparseThreshold_);
5432  perturbation_ = saved.perturbation_;
5433  infeasibilityCost_ = saved.infeasibilityCost_;
5434  dualBound_ = saved.dualBound_;
5435  delete progress_;
5436  progress_=NULL;
5437}
5438// To flag a variable (not inline to allow for column generation)
5439void 
5440ClpSimplex::setFlagged( int sequence)
5441{
5442  status_[sequence] |= 64;
5443  matrix_->generalExpanded(this,7,sequence);
5444  lastFlaggedIteration_=numberIterations_;
5445}
5446/* Factorizes and returns true if optimal.  Used by user */
5447bool
5448ClpSimplex::statusOfProblem(bool initial)
5449{
5450  // is factorization okay?
5451  if (initial) {
5452    // First time - allow singularities
5453    int numberThrownOut = -1;
5454    int totalNumberThrownOut=0;
5455    while(numberThrownOut) {
5456      int status = internalFactorize(0);
5457      if (status<0)
5458        return false; // some error
5459      else
5460        numberThrownOut = status;
5461     
5462      // for this we need clean basis so it is after factorize
5463      //if (!numberThrownOut)
5464      //numberThrownOut = gutsOfSolution(  NULL,NULL,
5465      //                                   false);
5466      //else
5467      //matrix_->rhsOffset(this,true); // redo rhs offset
5468      totalNumberThrownOut+= numberThrownOut;
5469     
5470    }
5471   
5472    if (totalNumberThrownOut)
5473      handler_->message(CLP_SINGULARITIES,messages_)
5474        <<totalNumberThrownOut
5475        <<CoinMessageEol;
5476  } else {
5477#ifndef NDEBUG   
5478    int returnCode=internalFactorize(1);
5479    assert (!returnCode);
5480#else
5481    internalFactorize(1);
5482#endif
5483  }
5484  // put back original costs and then check
5485  // also move to work arrays
5486  createRim(4+32);
5487  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
5488  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
5489  gutsOfSolution(NULL,NULL);
5490  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
5491  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
5492  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
5493  deleteRim(-1);
5494  return (primalFeasible()&&dualFeasible());
5495}
5496/* Return model - updates any scalars */
5497void 
5498ClpSimplex::returnModel(ClpSimplex & otherModel)
5499{
5500  ClpModel::returnModel(otherModel);
5501  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
5502  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
5503  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
5504  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
5505  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
5506  otherModel.columnDualSequence_ = columnDualSequence_;
5507  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
5508  otherModel.rowDualSequence_ = rowDualSequence_;
5509  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5510  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
5511  otherModel.largestPrimalError_ = largestPrimalError_;
5512  otherModel.largestDualError_ = largestDualError_;
5513  otherModel.largestSolutionError_ = largestSolutionError_;
5514  otherModel.alpha_ = alpha_;
5515  otherModel.theta_ = theta_;
5516  otherModel.lowerIn_ = lowerIn_;
5517  otherModel.valueIn_ = valueIn_;
5518  otherModel.upperIn_ = upperIn_;
5519  otherModel.dualIn_ = dualIn_;
5520  otherModel.sequenceIn_ = sequenceIn_;
5521  otherModel.directionIn_ = directionIn_;
5522  otherModel.lowerOut_ = lowerOut_;
5523  otherModel.valueOut_ = valueOut_;
5524  otherModel.upperOut_ = upperOut_;
5525  otherModel.dualOut_ = dualOut_;
5526  otherModel.sequenceOut_ = sequenceOut_;
5527  otherModel.directionOut_ = directionOut_;
5528  otherModel.pivotRow_ = pivotRow_;
5529  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
5530  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
5531  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
5532    numberDualInfeasibilitiesWithoutFree_;
5533  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
5534  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
5535  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
5536  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
5537  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
5538}
5539/* Constructs a non linear cost from list of non-linearities (columns only)
5540   First lower of each column is taken as real lower
5541   Last lower is taken as real upper and cost ignored
5542   
5543   Returns nonzero if bad data e.g. lowers not monotonic
5544*/
5545int 
5546ClpSimplex::createPiecewiseLinearCosts(const int * starts,
5547                                       const double * lower, const double * gradient)
5548{
5549  delete nonLinearCost_;
5550  // Set up feasible bounds and check monotonicity
5551  int iColumn;
5552  int returnCode=0;
5553
5554  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5555    int iIndex = starts[iColumn];
5556    int end = starts[iColumn+1]-1;
5557    columnLower_[iColumn] = lower[iIndex];
5558    columnUpper_[iColumn] = lower[end];
5559    double value = columnLower_[iColumn];
5560    iIndex++;
5561    for (;iIndex<end;iIndex++) {
5562      if (lower[iIndex]<value)
5563        returnCode++; // not monotonic
5564      value=lower[iIndex];
5565    }
5566  }
5567  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
5568  specialOptions_ |= 2; // say keep
5569  return returnCode;
5570}
5571/* For advanced use.  When doing iterative solves things can get
5572   nasty so on values pass if incoming solution has largest
5573   infeasibility < incomingInfeasibility throw out variables
5574   from basis until largest infeasibility < allowedInfeasibility
5575   or incoming largest infeasibility.
5576   If allowedInfeasibility>= incomingInfeasibility this is
5577   always possible altough you may end up with an all slack basis.
5578   
5579   Defaults are 1.0,10.0
5580*/
5581void 
5582ClpSimplex::setValuesPassAction(float incomingInfeasibility,
5583                                float allowedInfeasibility)
5584{
5585  incomingInfeasibility_=incomingInfeasibility;
5586  allowedInfeasibility_=allowedInfeasibility;
5587  assert(incomingInfeasibility_>=0.0);
5588  assert(allowedInfeasibility_>=incomingInfeasibility_);
5589}
5590//#############################################################################
5591
5592ClpSimplexProgress::ClpSimplexProgress () 
5593{
5594  int i;
5595  for (i=0;i<CLP_PROGRESS;i++) {
5596    objective_[i] = COIN_DBL_MAX;
5597    infeasibility_[i] = -1.0; // set to an impossible value
5598    realInfeasibility_[i] = COIN_DBL_MAX;
5599    numberInfeasibilities_[i]=-1; 
5600    iterationNumber_[i]=-1;
5601  }
5602  for (i=0;i<CLP_CYCLE;i++) {
5603    //obj_[i]=COIN_DBL_MAX;
5604    in_[i]=-1;
5605    out_[i]=-1;
5606    way_[i]=0;
5607  }
5608  numberTimes_ = 0;
5609  numberBadTimes_ = 0;
5610  model_ = NULL;
5611  oddState_=0;
5612}
5613
5614
5615//-----------------------------------------------------------------------------
5616
5617ClpSimplexProgress::~ClpSimplexProgress ()
5618{
5619}
5620// Copy constructor.
5621ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
5622{
5623  int i;
5624  for (i=0;i<CLP_PROGRESS;i++) {
5625    objective_[i] = rhs.objective_[i];
5626    infeasibility_[i] = rhs.infeasibility_[i];
5627    realInfeasibility_[i] = rhs.realInfeasibility_[i];
5628    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
5629    iterationNumber_[i]=rhs.iterationNumber_[i];
5630  }
5631  for (i=0;i<CLP_CYCLE;i++) {
5632    //obj_[i]=rhs.obj_[i];
5633    in_[i]=rhs.in_[i];
5634    out_[i]=rhs.out_[i];
5635    way_[i]=rhs.way_[i];
5636  }
5637  numberTimes_ = rhs.numberTimes_;
5638  numberBadTimes_ = rhs.numberBadTimes_;
5639  model_ = rhs.model_;
5640  oddState_=rhs.oddState_;
5641}
5642// Copy constructor.from model
5643ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
5644{
5645  model_ = model;
5646  int i;
5647  for (i=0;i<CLP_PROGRESS;i++) {
5648    if (model_->algorithm()>=0)
5649      objective_[i] = COIN_DBL_MAX;
5650    else
5651      objective_[i] = -COIN_DBL_MAX;
5652    infeasibility_[i] = -1.0; // set to an impossible value
5653    realInfeasibility_[i] = COIN_DBL_MAX;
5654    numberInfeasibilities_[i]=-1; 
5655    iterationNumber_[i]=-1;
5656  }
5657  for (i=0;i<CLP_CYCLE;i++) {
5658    //obj_[i]=COIN_DBL_MAX;
5659    in_[i]=-1;
5660    out_[i]=-1;
5661    way_[i]=0;
5662  }
5663  numberTimes_ = 0;
5664  numberBadTimes_ = 0;
5665  oddState_=0;
5666}
5667// Assignment operator. This copies the data
5668ClpSimplexProgress & 
5669ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
5670{
5671  if (this != &rhs) {
5672    int i;
5673    for (i=0;i<CLP_PROGRESS;i++) {
5674      objective_[i] = rhs.objective_[i];
5675      infeasibility_[i] = rhs.infeasibility_[i];
5676      realInfeasibility_[i] = rhs.realInfeasibility_[i];
5677      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
5678      iterationNumber_[i]=rhs.iterationNumber_[i];
5679    }
5680    for (i=0;i<CLP_CYCLE;i++) {
5681      //obj_[i]=rhs.obj_[i];
5682      in_[i]=rhs.in_[i