source: trunk/ClpSimplex.cpp @ 348

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

New status 5

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