source: trunk/ClpSimplex.cpp @ 372

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

Moving scaling to model plus robustness fixes in simplex
still messing with interior point

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