source: trunk/ClpSimplex.cpp @ 468

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

to make sbb faster

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