source: trunk/ClpSimplex.cpp @ 423

Last change on this file since 423 was 423, checked in by forrest, 15 years ago

various Cholesky

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 189.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/* Loads a problem (the constraints on the
4487   rows are given by lower and upper bounds). If a pointer is 0 then the
4488   following values are the default:
4489   <ul>
4490   <li> <code>colub</code>: all columns have upper bound infinity
4491   <li> <code>collb</code>: all columns have lower bound 0
4492   <li> <code>rowub</code>: all rows have upper bound infinity
4493   <li> <code>rowlb</code>: all rows have lower bound -infinity
4494   <li> <code>obj</code>: all variables have 0 objective coefficient
4495   </ul>
4496*/
4497void 
4498ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
4499                    const double* collb, const double* colub,   
4500                    const double* obj,
4501                    const double* rowlb, const double* rowub,
4502                    const double * rowObjective)
4503{
4504  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4505                        rowObjective);
4506  createStatus();
4507}
4508void 
4509ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
4510                    const double* collb, const double* colub,   
4511                    const double* obj,
4512                    const double* rowlb, const double* rowub,
4513                    const double * rowObjective)
4514{
4515  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4516                        rowObjective);
4517  createStatus();
4518}
4519
4520/* Just like the other loadProblem() method except that the matrix is
4521   given in a standard column major ordered format (without gaps). */
4522void 
4523ClpSimplex::loadProblem (  const int numcols, const int numrows,
4524                    const CoinBigIndex* start, const int* index,
4525                    const double* value,
4526                    const double* collb, const double* colub,   
4527                    const double* obj,
4528                    const double* rowlb, const double* rowub,
4529                    const double * rowObjective)
4530{
4531  ClpModel::loadProblem(numcols, numrows, start, index, value,
4532                          collb, colub, obj, rowlb, rowub,
4533                          rowObjective);
4534  createStatus();
4535}
4536void 
4537ClpSimplex::loadProblem (  const int numcols, const int numrows,
4538                           const CoinBigIndex* start, const int* index,
4539                           const double* value,const int * length,
4540                           const double* collb, const double* colub,   
4541                           const double* obj,
4542                           const double* rowlb, const double* rowub,
4543                           const double * rowObjective)
4544{
4545  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
4546                          collb, colub, obj, rowlb, rowub,
4547                          rowObjective);
4548  createStatus();
4549}
4550// Read an mps file from the given filename
4551int 
4552ClpSimplex::readMps(const char *filename,
4553            bool keepNames,
4554            bool ignoreErrors)
4555{
4556  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
4557  createStatus();
4558  return status;
4559}
4560// Just check solution (for external use)
4561void 
4562ClpSimplex::checkSolution()
4563{
4564  // put in standard form
4565  createRim(7+8+16);
4566  dualTolerance_=dblParam_[ClpDualTolerance];
4567  primalTolerance_=dblParam_[ClpPrimalTolerance];
4568  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
4569  checkDualSolution();
4570  if (!numberDualInfeasibilities_&&
4571      !numberPrimalInfeasibilities_)
4572    problemStatus_=0;
4573  else
4574    problemStatus_=-1;
4575#ifdef CLP_DEBUG
4576  int i;
4577  double value=0.0;
4578  for (i=0;i<numberRows_+numberColumns_;i++)
4579    value += dj_[i]*solution_[i];
4580  printf("dual value %g, primal %g\n",value,objectiveValue());
4581#endif
4582  // release extra memory
4583  deleteRim(0);
4584}
4585/* Crash - at present just aimed at dual, returns
4586   -2 if dual preferred and crash basis created
4587   -1 if dual preferred and all slack basis preferred
4588   0 if basis going in was not all slack
4589   1 if primal preferred and all slack basis preferred
4590   2 if primal preferred and crash basis created.
4591   
4592   if gap between bounds <="gap" variables can be flipped
4593   
4594   If "pivot" is
4595   0 No pivoting (so will just be choice of algorithm)
4596   1 Simple pivoting e.g. gub
4597   2 Mini iterations
4598*/
4599int 
4600ClpSimplex::crash(double gap,int pivot)
4601{
4602  assert(!rowObjective_); // not coded
4603  int iColumn;
4604  int numberBad=0;
4605  int numberBasic=0;
4606  double dualTolerance=dblParam_[ClpDualTolerance];
4607  //double primalTolerance=dblParam_[ClpPrimalTolerance];
4608  int returnCode=0;
4609  // If no basis then make all slack one
4610  if (!status_)
4611    createStatus();
4612 
4613  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4614    if (getColumnStatus(iColumn)==basic)
4615      numberBasic++;
4616  }
4617  if (!numberBasic) {
4618    // all slack
4619    double * dj = new double [numberColumns_];
4620    double * solution = columnActivity_;
4621    const double * linearObjective = objective();
4622    //double objectiveValue=0.0;
4623    int iColumn;
4624    double direction = optimizationDirection_;
4625    // direction is actually scale out not scale in
4626    if (direction)
4627      direction = 1.0/direction;
4628    for (iColumn=0;iColumn<numberColumns_;iColumn++)
4629      dj[iColumn] = direction*linearObjective[iColumn];
4630    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4631      // assume natural place is closest to zero
4632      double lowerBound = columnLower_[iColumn];
4633      double upperBound = columnUpper_[iColumn];
4634      if (lowerBound>-1.0e20||upperBound<1.0e20) {
4635        bool atLower;
4636        if (fabs(upperBound)<fabs(lowerBound)) {
4637          atLower=false;
4638          setColumnStatus(iColumn,atUpperBound);
4639          solution[iColumn]=upperBound;
4640        } else {
4641          atLower=true;
4642          setColumnStatus(iColumn,atLowerBound);
4643          solution[iColumn]=lowerBound;
4644        }
4645        if (dj[iColumn]<0.0) {
4646          // should be at upper bound
4647          if (atLower) {
4648            // can we flip
4649            if (upperBound-lowerBound<=gap) {
4650              columnActivity_[iColumn]=upperBound;
4651              setColumnStatus(iColumn,atUpperBound);
4652            } else if (dj[iColumn]<-dualTolerance) {
4653              numberBad++;
4654            }
4655          }
4656        } else if (dj[iColumn]>0.0) {
4657          // should be at lower bound
4658          if (!atLower) {
4659            // can we flip
4660            if (upperBound-lowerBound<=gap) {
4661              columnActivity_[iColumn]=lowerBound;
4662              setColumnStatus(iColumn,atLowerBound);
4663            } else if (dj[iColumn]>dualTolerance) {
4664              numberBad++;
4665            }
4666          }
4667        }
4668      } else {
4669        // free
4670        setColumnStatus(iColumn,isFree);
4671        if (fabs(dj[iColumn])>dualTolerance) 
4672          numberBad++;
4673      }
4674    }
4675    if (numberBad||pivot) {
4676      if (!pivot) {
4677        delete [] dj;
4678        returnCode = 1;
4679      } else {
4680        // see if can be made dual feasible with gubs etc
4681        double * pi = new double[numberRows_];
4682        memset (pi,0,numberRows_*sizeof(double));
4683        int * way = new int[numberColumns_];
4684        int numberIn = 0;
4685
4686        // Get column copy
4687        CoinPackedMatrix * columnCopy = matrix();
4688        // Get a row copy in standard format
4689        CoinPackedMatrix copy;
4690        copy.reverseOrderedCopyOf(*columnCopy);
4691        // get matrix data pointers
4692        const int * column = copy.getIndices();
4693        const CoinBigIndex * rowStart = copy.getVectorStarts();
4694        const int * rowLength = copy.getVectorLengths(); 
4695        const double * elementByRow = copy.getElements();
4696        //const int * row = columnCopy->getIndices();
4697        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4698        //const int * columnLength = columnCopy->getVectorLengths();
4699        //const double * element = columnCopy->getElements();
4700
4701
4702        // if equality row and bounds mean artificial in basis bad
4703        // then do anyway
4704
4705        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4706          // - if we want to reduce dj, + if we want to increase
4707          int thisWay = 100;
4708          double lowerBound = columnLower_[iColumn];
4709          double upperBound = columnUpper_[iColumn];
4710          if (upperBound>lowerBound) {
4711            switch(getColumnStatus(iColumn)) {
4712             
4713            case basic:
4714              thisWay=0;
4715            case ClpSimplex::isFixed:
4716              break;
4717            case isFree:
4718            case superBasic:
4719              if (dj[iColumn]<-dualTolerance) 
4720                thisWay = 1;
4721              else if (dj[iColumn]>dualTolerance) 
4722                thisWay = -1;
4723              else
4724                thisWay =0;
4725              break;
4726            case atUpperBound:
4727              if (dj[iColumn]>dualTolerance) 
4728                thisWay = -1;
4729              else if (dj[iColumn]<-dualTolerance) 
4730                thisWay = -3;
4731              else
4732                thisWay = -2;
4733              break;
4734            case atLowerBound:
4735              if (dj[iColumn]<-dualTolerance) 
4736                thisWay = 1;
4737              else if (dj[iColumn]>dualTolerance) 
4738                thisWay = 3;
4739              else
4740                thisWay = 2;
4741              break;
4742            }
4743          }
4744          way[iColumn] = thisWay;
4745        }
4746        /*if (!numberBad)
4747          printf("Was dual feasible before passes - rows %d\n",
4748          numberRows_);*/
4749        int lastNumberIn = -100000;
4750        int numberPasses=5;
4751        while (numberIn>lastNumberIn+numberRows_/100) {
4752          lastNumberIn = numberIn;
4753          // we need to maximize chance of doing good
4754          int iRow;
4755          for (iRow=0;iRow<numberRows_;iRow++) {
4756            double lowerBound = rowLower_[iRow];
4757            double upperBound = rowUpper_[iRow];
4758            if (getRowStatus(iRow)==basic) {
4759              // see if we can find a column to pivot on
4760              int j;
4761              // down is amount pi can go down
4762              double maximumDown = COIN_DBL_MAX;
4763              double maximumUp = COIN_DBL_MAX;
4764              double minimumDown =0.0;
4765              double minimumUp =0.0;
4766              int iUp=-1;
4767              int iDown=-1;
4768              int iUpB=-1;
4769              int iDownB=-1;
4770              if (lowerBound<-1.0e20)
4771                maximumUp = -1.0;
4772              if (upperBound>1.0e20)
4773                maximumDown = -1.0;
4774              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4775                int iColumn = column[j];
4776                double value = elementByRow[j];
4777                double djValue = dj[iColumn];
4778                /* way -
4779                   -3 - okay at upper bound with negative dj
4780                   -2 - marginal at upper bound with zero dj - can only decrease
4781                   -1 - bad at upper bound
4782                   0 - we can never pivot on this row
4783                   1 - bad at lower bound
4784                   2 - marginal at lower bound with zero dj - can only increase
4785                   3 - okay at lower bound with positive dj
4786                   100 - fine we can just ignore
4787                */
4788                if (way[iColumn]!=100) {
4789                  switch(way[iColumn]) {
4790                   
4791                  case -3:
4792                    if (value>0.0) {
4793                      if (maximumDown*value>-djValue) {
4794                        maximumDown = -djValue/value;
4795                        iDown = iColumn;
4796                      }
4797                    } else {
4798                      if (-maximumUp*value>-djValue) {
4799                        maximumUp = djValue/value;
4800                        iUp = iColumn;
4801                      }
4802                    }
4803                    break;
4804                  case -2:
4805                    if (value>0.0) {
4806                      maximumDown = 0.0;
4807                    } else {
4808                      maximumUp = 0.0;
4809                    }
4810                    break;
4811                  case -1:
4812                    // see if could be satisfied
4813                    // dj value > 0
4814                    if (value>0.0) {
4815                      maximumDown=0.0;
4816                      if (maximumUp*value<djValue-dualTolerance) {
4817                        maximumUp = 0.0; // would improve but not enough
4818                      } else {
4819                        if (minimumUp*value<djValue) {
4820                          minimumUp = djValue/value;
4821                          iUpB = iColumn;
4822                        }
4823                      }
4824                    } else {
4825                      maximumUp=0.0;
4826                      if (-maximumDown*value<djValue-dualTolerance) {
4827                        maximumDown = 0.0; // would improve but not enough
4828                      } else {
4829                        if (-minimumDown*value<djValue) {
4830                          minimumDown = -djValue/value;
4831                          iDownB = iColumn;
4832                        }
4833                      }
4834                    }
4835                   
4836                    break;
4837                  case 0:
4838                    maximumDown = -1.0;
4839                    maximumUp=-1.0;
4840                    break;
4841                  case 1:
4842                    // see if could be satisfied
4843                    // dj value < 0
4844                    if (value>0.0) {
4845                      maximumUp=0.0;
4846                      if (maximumDown*value<-djValue-dualTolerance) {
4847                        maximumDown = 0.0; // would improve but not enough
4848                      } else {
4849                        if (minimumDown*value<-djValue) {
4850                          minimumDown = -djValue/value;
4851                          iDownB = iColumn;
4852                        }
4853                      }
4854                    } else {
4855                      maximumDown=0.0;
4856                      if (-maximumUp*value<-djValue-dualTolerance) {
4857                        maximumUp = 0.0; // would improve but not enough
4858                      } else {
4859                        if (-minimumUp*value<-djValue) {
4860                          minimumUp = djValue/value;
4861                          iUpB = iColumn;
4862                        }
4863                      }
4864                    }
4865                   
4866                    break;
4867                  case 2:
4868                    if (value>0.0) {
4869                      maximumUp = 0.0;
4870                    } else {
4871                      maximumDown = 0.0;
4872                    }
4873                   
4874                    break;
4875                  case 3:
4876                    if (value>0.0) {
4877                      if (maximumUp*value>djValue) {
4878                        maximumUp = djValue/value;
4879                        iUp = iColumn;
4880                      }
4881                    } else {
4882                      if (-maximumDown*value>djValue) {
4883                        maximumDown = -djValue/value;
4884                        iDown = iColumn;
4885                      }
4886                    }
4887                   
4888                    break;
4889                  default:
4890                    break;
4891                  }
4892                }
4893              }
4894              if (iUpB>=0)
4895                iUp=iUpB;
4896              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
4897                iUp=-1;
4898              if (iDownB>=0)
4899                iDown=iDownB;
4900              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
4901                iDown=-1;
4902              if (iUp>=0||iDown>=0) {
4903                // do something
4904                if (iUp>=0&&iDown>=0) {
4905                  if (maximumDown>maximumUp)
4906                    iUp=-1;
4907                }
4908                double change;
4909                int kColumn;
4910                if (iUp>=0) {
4911                  kColumn=iUp;
4912                  change=maximumUp;
4913                  // just do minimum if was dual infeasible
4914                  // ? only if maximum large?
4915                  if (minimumUp>0.0)
4916                    change=minimumUp;
4917                  setRowStatus(iRow,atUpperBound);
4918                } else {
4919                  kColumn=iDown;
4920                  change=-maximumDown;
4921                  // just do minimum if was dual infeasible
4922                  // ? only if maximum large?
4923                  if (minimumDown>0.0)
4924                    change=-minimumDown;
4925                  setRowStatus(iRow,atLowerBound);
4926                }
4927                assert (fabs(change)<1.0e20);
4928                setColumnStatus(kColumn,basic);
4929                numberIn++;
4930                pi[iRow]=change;
4931                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4932                  int iColumn = column[j];
4933                  double value = elementByRow[j];
4934                  double djValue = dj[iColumn]-change*value;
4935                  dj[iColumn]=djValue;
4936                  if (abs(way[iColumn])==1) {
4937                    numberBad--;
4938                    /*if (!numberBad)
4939                      printf("Became dual feasible at row %d out of %d\n",
4940                      iRow, numberRows_);*/
4941                    lastNumberIn=-1000000;
4942                  }
4943                  int thisWay = 100;
4944                  double lowerBound = columnLower_[iColumn];
4945                  double upperBound = columnUpper_[iColumn];
4946                  if (upperBound>lowerBound) {
4947                    switch(getColumnStatus(iColumn)) {
4948                     
4949                    case basic:
4950                      thisWay=0;
4951                    case isFixed:
4952                      break;
4953                    case isFree:
4954                    case superBasic:
4955                      if (djValue<-dualTolerance) 
4956                        thisWay = 1;
4957                      else if (djValue>dualTolerance) 
4958                        thisWay = -1;
4959                      else
4960                        { thisWay =0; abort();}
4961                      break;
4962                    case atUpperBound:
4963                      if (djValue>dualTolerance) 
4964                        { thisWay =-1; abort();}
4965                      else if (djValue<-dualTolerance) 
4966                        thisWay = -3;
4967                      else
4968                        thisWay = -2;
4969                      break;
4970                    case atLowerBound:
4971                      if (djValue<-dualTolerance) 
4972                        { thisWay =1; abort();}
4973                      else if (djValue>dualTolerance) 
4974                        thisWay = 3;
4975                      else
4976                        thisWay = 2;
4977                      break;
4978                    }
4979                  }
4980                  way[iColumn] = thisWay;
4981                }
4982              }
4983            }
4984          }
4985          if (numberIn==lastNumberIn||numberBad||pivot<2)
4986            break;
4987          if (!(--numberPasses))
4988            break;
4989          //printf("%d put in so far\n",numberIn);
4990        }
4991        // last attempt to flip
4992        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4993          double lowerBound = columnLower_[iColumn];
4994          double upperBound = columnUpper_[iColumn];
4995          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
4996            double djValue=dj[iColumn];
4997            switch(getColumnStatus(iColumn)) {
4998             
4999            case basic:
5000            case ClpSimplex::isFixed:
5001              break;
5002            case isFree:
5003            case superBasic:
5004              break;
5005            case atUpperBound:
5006              if (djValue>dualTolerance) {
5007                setColumnStatus(iColumn,atUpperBound);
5008                solution[iColumn]=upperBound;
5009              } 
5010              break;
5011            case atLowerBound:
5012              if (djValue<-dualTolerance) {
5013                setColumnStatus(iColumn,atUpperBound);
5014                solution[iColumn]=upperBound;
5015              }
5016              break;
5017            }
5018          }
5019        }
5020        delete [] pi;
5021        delete [] dj;
5022        delete [] way;
5023        handler_->message(CLP_CRASH,messages_)
5024          <<numberIn
5025          <<numberBad
5026          <<CoinMessageEol;
5027        returnCode =  -1;
5028      }
5029    } else {
5030      delete [] dj;
5031      returnCode =  -1;
5032    }
5033    //cleanStatus();
5034  }
5035  return returnCode;
5036}
5037/* Pivot in a variable and out a variable.  Returns 0 if okay,
5038   1 if inaccuracy forced re-factorization, -1 if would be singular.
5039   Also updates primal/dual infeasibilities.
5040   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
5041*/
5042int ClpSimplex::pivot()
5043{
5044  // scaling not allowed
5045  assert (!scalingFlag_);
5046  // assume In_ and Out_ are correct and directionOut_ set
5047  // (or In_ if flip
5048  lowerIn_ = lower_[sequenceIn_];
5049  valueIn_ = solution_[sequenceIn_];
5050  upperIn_ = upper_[sequenceIn_];
5051  dualIn_ = dj_[sequenceIn_];
5052  lowerOut_ = lower_[sequenceOut_];
5053  valueOut_ = solution_[sequenceOut_];
5054  upperOut_ = upper_[sequenceOut_];
5055  // for now assume primal is feasible (or in dual)
5056  dualOut_ = dj_[sequenceOut_];
5057  assert(fabs(dualOut_)<1.0e-6);
5058  bool roundAgain = true;
5059  int returnCode=0;
5060  while (roundAgain) {
5061    roundAgain=false;
5062    unpack(rowArray_[1]);
5063    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
5064    alpha_=0.0;
5065    int i;
5066    int * index = rowArray_[1]->getIndices();
5067    int number = rowArray_[1]->getNumElements();
5068    double * element = rowArray_[1]->denseVector();
5069    for (i=0;i<number;i++) {
5070      int ii = index[i];
5071      if ( pivotVariable_[ii]==sequenceOut_) {
5072        pivotRow_=ii;
5073        alpha_=element[pivotRow_];
5074        break;
5075      }
5076    }
5077    assert (fabs(alpha_)>1.0e-12);
5078    // we are going to subtract movement from current basic
5079    double movement;
5080    // see where incoming will go to
5081    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
5082      // flip so go to bound
5083      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
5084    } else {
5085      // get where outgoing needs to get to
5086      double outValue = (directionOut_<0) ? upperOut_ : lowerOut_;
5087      // solutionOut_ - movement*alpha_ == outValue
5088      movement = (valueOut_-outValue)/alpha_;
5089      // set directionIn_ correctly
5090      directionIn_ = (movement>0) ? 1 :-1;
5091    }
5092    // update primal solution
5093    for (i=0;i<number;i++) {
5094      int ii = index[i];
5095      // get column
5096      int ij = pivotVariable_[ii];
5097      solution_[ij] -= movement*element[ii];
5098      element[ii]=0.0;
5099    }
5100    rowArray_[1]->setNumElements(0);
5101    // see where something went to
5102    if (sequenceOut_<0) {
5103      if (directionIn_<0) {
5104        assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
5105        solution_[sequenceIn_]=upperIn_;
5106      } else {
5107        assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
5108        solution_[sequenceIn_]=lowerIn_;
5109      }
5110    } else {
5111      if (directionOut_<0) {
5112        assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
5113        solution_[sequenceOut_]=upperOut_;
5114      } else {
5115        assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
5116        solution_[sequenceOut_]=lowerOut_;
5117      }
5118      solution_[sequenceIn_]=valueIn_+movement;
5119    }
5120    double objectiveChange = dualIn_*movement;
5121    // update duals
5122    if (pivotRow_>=0) {
5123      assert (fabs(alpha_)>1.0e-8);
5124      double multiplier = dualIn_/alpha_;
5125      rowArray_[0]->insert(pivotRow_,multiplier);
5126      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
5127      // put row of tableau in rowArray[0] and columnArray[0]
5128      matrix_->transposeTimes(this,-1.0,
5129                              rowArray_[0],columnArray_[1],columnArray_[0]);
5130      // update column djs
5131      int i;
5132      int * index = columnArray_[0]->getIndices();
5133      int number = columnArray_[0]->getNumElements();
5134      double * element = columnArray_[0]->denseVector();
5135      for (i=0;i<number;i++) {
5136        int ii = index[i];
5137        dj_[ii] += element[ii];
5138        element[ii]=0.0;
5139      }
5140      columnArray_[0]->setNumElements(0);
5141      // and row djs
5142      index = rowArray_[0]->getIndices();
5143      number = rowArray_[0]->getNumElements();
5144      element = rowArray_[0]->denseVector();
5145      for (i=0;i<number;i++) {
5146        int ii = index[i];
5147        dj_[ii+numberColumns_] += element[ii];
5148        dual_[ii] = dj_[ii+numberColumns_];
5149        element[ii]=0.0;
5150      }
5151      rowArray_[0]->setNumElements(0);
5152      // check incoming
5153      assert (fabs(dj_[sequenceIn_])<1.0e-6);
5154    }
5155   
5156    // if stable replace in basis
5157    int updateStatus = factorization_->replaceColumn(this,
5158                                                     rowArray_[2],
5159                                                     rowArray_[1],
5160                                                     pivotRow_,
5161                                                     alpha_);
5162    bool takePivot=true;
5163    // if no pivots, bad update but reasonable alpha - take and invert
5164    if (updateStatus==2&&
5165        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
5166      updateStatus=4;
5167    if (updateStatus==1||updateStatus==4) {
5168      // slight error
5169      if (factorization_->pivots()>5||updateStatus==4) {
5170        returnCode=-1;
5171      }
5172    } else if (updateStatus==2) {
5173      // major error
5174      rowArray_[1]->clear();
5175      takePivot=false;
5176      if (factorization_->pivots()) {
5177        // refactorize here
5178        statusOfProblem();
5179        roundAgain=true;
5180      } else {
5181        returnCode=1;
5182      }
5183    } else if (updateStatus==3) {
5184      // out of memory
5185      // increase space if not many iterations
5186      if (factorization_->pivots()<
5187          0.5*factorization_->maximumPivots()&&
5188          factorization_->pivots()<200)
5189        factorization_->areaFactor(
5190                                   factorization_->areaFactor() * 1.1);
5191      returnCode =-1; // factorize now
5192    }
5193    if (takePivot) {
5194      int save = algorithm_;
5195      // make simple so always primal
5196      algorithm_=1;
5197      housekeeping(objectiveChange);
5198      algorithm_=save;
5199    }
5200  }
5201  if (returnCode == -1) {
5202    // refactorize here
5203    statusOfProblem();
5204  } else {
5205    // just for now - recompute anyway
5206    gutsOfSolution(NULL,NULL);
5207  }
5208  return returnCode;
5209}
5210
5211/* Pivot in a variable and choose an outgoing one.  Assumes primal
5212   feasible - will not go through a bound.  Returns step length in theta
5213   Returns ray in ray_ (or NULL if no pivot)
5214   Return codes as before but -1 means no acceptable pivot
5215*/
5216int ClpSimplex::primalPivotResult()
5217{
5218  assert (sequenceIn_>=0);
5219  valueIn_=solution_[sequenceIn_];
5220  lowerIn_=lower_[sequenceIn_];
5221  upperIn_=upper_[sequenceIn_];
5222  dualIn_=dj_[sequenceIn_];
5223
5224  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
5225  if (returnCode<0&&returnCode>-4) {
5226    return 0;
5227  } else {
5228    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
5229           returnCode);
5230    return -1;
5231  }
5232}
5233 
5234/* Pivot out a variable and choose an incoing one.  Assumes dual
5235   feasible - will not go through a reduced cost. 
5236   Returns step length in theta
5237   Returns ray in ray_ (or NULL if no pivot)
5238   Return codes as before but -1 means no acceptable pivot
5239*/
5240int 
5241ClpSimplex::dualPivotResult()
5242{
5243  return ((ClpSimplexDual *) this)->pivotResult();
5244}
5245// Factorization frequency
5246int 
5247ClpSimplex::factorizationFrequency() const
5248{
5249  if (factorization_)
5250    return factorization_->maximumPivots();
5251  else 
5252    return -1;
5253}
5254void 
5255ClpSimplex::setFactorizationFrequency(int value)
5256{
5257  if (factorization_)
5258    factorization_->maximumPivots(value);
5259}
5260// Common bits of coding for dual and primal
5261int 
5262ClpSimplex::startup(int ifValuesPass)
5263{
5264  // sanity check
5265  // bad if empty (trap here to avoid using bad matrix_)
5266  if (!matrix_||(!matrix_->getNumElements()&&objective_->type()<2)) {
5267    int infeasNumber[2];
5268    double infeasSum[2];
5269    problemStatus_=emptyProblem(infeasNumber,infeasSum);
5270    numberDualInfeasibilities_=infeasNumber[0];
5271    sumDualInfeasibilities_=infeasSum[0];
5272    numberPrimalInfeasibilities_=infeasNumber[1];
5273    sumPrimalInfeasibilities_=infeasSum[1];
5274    return 2;
5275  }
5276  pivotRow_=-1;
5277  sequenceIn_=-1;
5278  sequenceOut_=-1;
5279  secondaryStatus_=0;
5280
5281  primalTolerance_=dblParam_[ClpPrimalTolerance];
5282  dualTolerance_=dblParam_[ClpDualTolerance];
5283  if (problemStatus_!=10)
5284    numberIterations_=0;
5285
5286  // put in standard form (and make row copy)
5287  // create modifiable copies of model rim and do optional scaling
5288  bool goodMatrix=createRim(7+8+16,true);
5289
5290  if (goodMatrix) {
5291    // Model looks okay
5292    // Do initial factorization
5293    // and set certain stuff
5294    // We can either set increasing rows so ...IsBasic gives pivot row
5295    // or we can just increment iBasic one by one
5296    // for now let ...iBasic give pivot row
5297    factorization_->increasingRows(2);
5298    // row activities have negative sign
5299    factorization_->slackValue(-1.0);
5300    factorization_->zeroTolerance(1.0e-13);
5301    // Switch off dense (unless special option set)
5302    int saveThreshold = factorization_->denseThreshold();
5303    factorization_->setDenseThreshold(0);
5304    // If values pass then perturb (otherwise may be optimal so leave a bit)
5305    if (ifValuesPass) {
5306      // do perturbation if asked for
5307     
5308      if (perturbation_<100) {
5309        if (algorithm_>0&&(objective_->type()<2||!objective_->activated())) {
5310          ((ClpSimplexPrimal *) this)->perturb(0);
5311        } else if (algorithm_<0) {
5312        ((ClpSimplexDual *) this)->perturb();
5313        }
5314      }
5315    }
5316    // for primal we will change bounds using infeasibilityCost_
5317    if (nonLinearCost_==NULL&&algorithm_>0) {
5318      // get a valid nonlinear cost function
5319      delete nonLinearCost_;
5320      nonLinearCost_= new ClpNonLinearCost(this);
5321    }
5322   
5323    // loop round to clean up solution if values pass
5324    int numberThrownOut = -1;
5325    int totalNumberThrownOut=0;
5326    while(numberThrownOut) {
5327      int status = internalFactorize(ifValuesPass ? 10 : 0);
5328      if (status<0)
5329        return 1; // some error
5330      else
5331        numberThrownOut = status;
5332     
5333      // for this we need clean basis so it is after factorize
5334      if (!numberThrownOut)
5335        numberThrownOut = gutsOfSolution(  NULL,NULL,
5336                                         ifValuesPass!=0);
5337      else
5338        matrix_->rhsOffset(this,true); // redo rhs offset
5339      totalNumberThrownOut+= numberThrownOut;
5340     
5341    }
5342   
5343    if (totalNumberThrownOut)
5344      handler_->message(CLP_SINGULARITIES,messages_)
5345        <<totalNumberThrownOut
5346        <<CoinMessageEol;
5347    // Switch back dense
5348    factorization_->setDenseThreshold(saveThreshold);
5349   
5350    problemStatus_ = -1;
5351   
5352    // number of times we have declared optimality
5353    numberTimesOptimal_=0;
5354
5355    return 0;
5356  } else {
5357    // bad matrix
5358    return 2;
5359  }
5360   
5361}
5362
5363
5364void 
5365ClpSimplex::finish()
5366{
5367  // Get rid of some arrays and empty factorization
5368  deleteRim();
5369  // Skip message if changing algorithms
5370  if (problemStatus_!=10) {
5371    if (problemStatus_==-1)
5372      problemStatus_=4;
5373    assert(problemStatus_>=0&&problemStatus_<6);
5374    handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
5375      <<objectiveValue()
5376      <<CoinMessageEol;
5377  }
5378  factorization_->relaxAccuracyCheck(1.0);
5379  // get rid of any network stuff - could do more
5380  factorization_->cleanUp();
5381}
5382// Save data
5383ClpDataSave
5384ClpSimplex::saveData() 
5385{
5386  ClpDataSave saved;
5387  saved.dualBound_ = dualBound_;
5388  saved.infeasibilityCost_ = infeasibilityCost_;
5389  saved.sparseThreshold_ = factorization_->sparseThreshold();
5390  saved.perturbation_ = perturbation_;
5391  // Progress indicator
5392  delete progress_;
5393  progress_ = new ClpSimplexProgress (this);
5394  return saved;
5395}
5396// Restore data
5397void 
5398ClpSimplex::restoreData(ClpDataSave saved)
5399{
5400  factorization_->sparseThreshold(saved.sparseThreshold_);
5401  perturbation_ = saved.perturbation_;
5402  infeasibilityCost_ = saved.infeasibilityCost_;
5403  dualBound_ = saved.dualBound_;
5404  delete progress_;
5405  progress_=NULL;
5406}
5407// To flag a variable (not inline to allow for column generation)
5408void 
5409ClpSimplex::setFlagged( int sequence)
5410{
5411  status_[sequence] |= 64;
5412  matrix_->generalExpanded(this,7,sequence);
5413  lastFlaggedIteration_=numberIterations_;
5414}
5415/* Factorizes and returns true if optimal.  Used by user */
5416bool
5417ClpSimplex::statusOfProblem(bool initial)
5418{
5419  // is factorization okay?
5420  if (initial) {
5421    // First time - allow singularities
5422    int numberThrownOut = -1;
5423    int totalNumberThrownOut=0;
5424    while(numberThrownOut) {
5425      int status = internalFactorize(0);
5426      if (status<0)
5427        return false; // some error
5428      else
5429        numberThrownOut = status;
5430     
5431      // for this we need clean basis so it is after factorize
5432      //if (!numberThrownOut)
5433      //numberThrownOut = gutsOfSolution(  NULL,NULL,
5434      //                                   false);
5435      //else
5436      //matrix_->rhsOffset(this,true); // redo rhs offset
5437      totalNumberThrownOut+= numberThrownOut;
5438     
5439    }
5440   
5441    if (totalNumberThrownOut)
5442      handler_->message(CLP_SINGULARITIES,messages_)
5443        <<totalNumberThrownOut
5444        <<CoinMessageEol;
5445  } else {
5446#ifndef NDEBUG   
5447    int returnCode=internalFactorize(1);
5448    assert (!returnCode);
5449#else
5450    internalFactorize(1);
5451#endif
5452  }
5453  // put back original costs and then check
5454  // also move to work arrays
5455  createRim(4+32);
5456  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
5457  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
5458  gutsOfSolution(NULL,NULL);
5459  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
5460  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
5461  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
5462  deleteRim(-1);
5463  return (primalFeasible()&&dualFeasible());
5464}
5465/* Return model - updates any scalars */
5466void 
5467ClpSimplex::returnModel(ClpSimplex & otherModel)
5468{
5469  ClpModel::returnModel(otherModel);
5470  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
5471  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
5472  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
5473  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
5474  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
5475  otherModel.columnDualSequence_ = columnDualSequence_;
5476  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
5477  otherModel.rowDualSequence_ = rowDualSequence_;
5478  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5479  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
5480  otherModel.largestPrimalError_ = largestPrimalError_;
5481  otherModel.largestDualError_ = largestDualError_;
5482  otherModel.largestSolutionError_ = largestSolutionError_;
5483  otherModel.alpha_ = alpha_;
5484  otherModel.theta_ = theta_;
5485  otherModel.lowerIn_ = lowerIn_;
5486  otherModel.valueIn_ = valueIn_;
5487  otherModel.upperIn_ = upperIn_;
5488  otherModel.dualIn_ = dualIn_;
5489  otherModel.sequenceIn_ = sequenceIn_;
5490  otherModel.directionIn_ = directionIn_;
5491  otherModel.lowerOut_ = lowerOut_;
5492  otherModel.valueOut_ = valueOut_;
5493  otherModel.upperOut_ = upperOut_;
5494  otherModel.dualOut_ = dualOut_;
5495  otherModel.sequenceOut_ = sequenceOut_;
5496  otherModel.directionOut_ = directionOut_;
5497  otherModel.pivotRow_ = pivotRow_;
5498  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
5499  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
5500  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
5501    numberDualInfeasibilitiesWithoutFree_;
5502  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
5503  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
5504  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
5505  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
5506  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
5507}
5508/* Constructs a non linear cost from list of non-linearities (columns only)
5509   First lower of each column is taken as real lower
5510   Last lower is taken as real upper and cost ignored
5511   
5512   Returns nonzero if bad data e.g. lowers not monotonic
5513*/
5514int 
5515ClpSimplex::createPiecewiseLinearCosts(const int * starts,
5516                                       const double * lower, const double * gradient)
5517{
5518  delete nonLinearCost_;
5519  // Set up feasible bounds and check monotonicity
5520  int iColumn;
5521  int returnCode=0;
5522
5523  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5524    int iIndex = starts[iColumn];
5525    int end = starts[iColumn+1]-1;
5526    columnLower_[iColumn] = lower[iIndex];
5527    columnUpper_[iColumn] = lower[end];
5528    double value = columnLower_[iColumn];
5529    iIndex++;
5530    for (;iIndex<end;iIndex++) {
5531      if (lower[iIndex]<value)
5532        returnCode++; // not monotonic
5533      value=lower[iIndex];
5534    }
5535  }
5536  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
5537  specialOptions_ |= 2; // say keep
5538  return returnCode;
5539}
5540/* For advanced use.  When doing iterative solves things can get
5541   nasty so on values pass if incoming solution has largest
5542   infeasibility < incomingInfeasibility throw out variables
5543   from basis until largest infeasibility < allowedInfeasibility
5544   or incoming largest infeasibility.
5545   If allowedInfeasibility>= incomingInfeasibility this is
5546   always possible altough you may end up with an all slack basis.
5547   
5548   Defaults are 1.0,10.0
5549*/
5550void 
5551ClpSimplex::setValuesPassAction(float incomingInfeasibility,
5552                                float allowedInfeasibility)
5553{
5554  incomingInfeasibility_=incomingInfeasibility;
5555  allowedInfeasibility_=allowedInfeasibility;
5556  assert(incomingInfeasibility_>=0.0);
5557  assert(allowedInfeasibility_>=incomingInfeasibility_);
5558}
5559//#############################################################################
5560
5561ClpSimplexProgress::ClpSimplexProgress () 
5562{
5563  int i;
5564  for (i=0;i<CLP_PROGRESS;i++) {
5565    objective_[i] = COIN_DBL_MAX;
5566    infeasibility_[i] = -1.0; // set to an impossible value
5567    realInfeasibility_[i] = COIN_DBL_MAX;
5568    numberInfeasibilities_[i]=-1; 
5569    iterationNumber_[i]=-1;
5570  }
5571  for (i=0;i<CLP_CYCLE;i++) {
5572    //obj_[i]=COIN_DBL_MAX;
5573    in_[i]=-1;
5574    out_[i]=-1;
5575    way_[i]=0;
5576  }
5577  numberTimes_ = 0;
5578  numberBadTimes_ = 0;
5579  model_ = NULL;
5580  oddState_=0;
5581}
5582
5583
5584//-----------------------------------------------------------------------------
5585
5586ClpSimplexProgress::~ClpSimplexProgress ()
5587{
5588}
5589// Copy constructor.
5590ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
5591{
5592  int i;
5593  for (i=0;i<CLP_PROGRESS;i++) {
5594    objective_[i] = rhs.objective_[i];
5595    infeasibility_[i] = rhs.infeasibility_[i];
5596    realInfeasibility_[i] = rhs.realInfeasibility_[i];
5597    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
5598    iterationNumber_[i]=rhs.iterationNumber_[i];
5599  }
5600  for (i=0;i<CLP_CYCLE;i++) {
5601    //obj_[i]=rhs.obj_[i];
5602    in_[i]=rhs.in_[i];
5603    out_[i]=rhs.out_[i];
5604    way_[i]=rhs.way_[i];
5605  }
5606  numberTimes_ = rhs.numberTimes_;
5607  numberBadTimes_ = rhs.numberBadTimes_;
5608  model_ = rhs.model_;
5609  oddState_=rhs.oddState_;
5610}
5611// Copy constructor.from model
5612ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
5613{
5614  model_ = model;
5615  int i;
5616  for (i=0;i<CLP_PROGRESS;i++) {
5617    if (model_->algorithm()>=0)
5618      objective_[i] = COIN_DBL_MAX;
5619    else
5620      objective_[i] = -COIN_DBL_MAX;
5621    infeasibility_[i] = -1.0; // set to an impossible value
5622    realInfeasibility_[i] = COIN_DBL_MAX;
5623    numberInfeasibilities_[i]=-1; 
5624    iterationNumber_[i]=-1;
5625  }
5626  for (i=0;i<CLP_CYCLE;i++) {
5627    //obj_[i]=COIN_DBL_MAX;
5628    in_[i]=-1;
5629    out_[i]=-1;
5630    way_[i]=0;
5631  }
5632  numberTimes_ = 0;
5633  numberBadTimes_ = 0;
5634  oddState_=0;
5635}
5636// Assignment operator. This copies the data
5637ClpSimplexProgress & 
5638ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
5639{
5640  if (this != &rhs) {
5641    int i;
5642    for (i=0;i<CLP_PROGRESS;i++) {
5643      objective_[i] = rhs.objective_[i];
5644      infeasibility_[i] = rhs.infeasibility_[i];
5645      realInfeasibility_[i] = rhs.realInfeasibility_[i];
5646      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
5647      iterationNumber_[i]=rhs.iterationNumber_[i];
5648    }
5649    for (i=0;i<CLP_CYCLE;i++) {
5650      //obj_[i]=rhs.obj_[i];
5651      in_[i]=rhs.in_[i];
5652      out_[i]=rhs.out_[i];
5653      way_[i]=rhs.way_[i];
5654    }
5655    numberTimes_ = rhs.numberTimes_;
5656    numberBadTimes_ = rhs.numberBadTimes_;
5657    model_ = rhs.model_;
5658    oddState_=rhs.oddState_;
5659  }
5660  return *this;
5661}
5662// Seems to be something odd about exact comparison of doubles on linux
5663static bool equalDouble(double value1, double value2)
5664{
5665  unsigned int *i1 = (unsigned int *) &value1;
5666  unsigned int *i2 = (unsigned int *) &value2;
5667  if (sizeof(unsigned int)*2==sizeof(double)) 
5668    return (i1[0]==i2[0]&&i1[1]==i2[1]);
5669  else
5670    return (i1[0]==i2[0]);
5671}
5672int
5673ClpSimplexProgress::looping()
5674{
5675  if (!model_)
5676    return -1;
5677  double objective = model_->objectiveValue();
5678  double infeasibility;
5679  double realInfeasibility=0.0;
5680  int numberInfeasibilities;
5681  int iterationNumber = model_->numberIterations();
5682  if (model_->algorithm()