source: trunk/ClpSimplex.cpp @ 460

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

a bit faster plus compiler bug in ClpsimplexPrimal?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 193.1 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  if ((what&(16+32))!=0) {
2252    // move information to work arrays
2253    double direction = optimizationDirection_;
2254    // direction is actually scale out not scale in
2255    if (direction)
2256      direction = 1.0/direction;
2257    if (direction!=1.0) {
2258      // reverse all dual signs
2259      for (i=0;i<numberColumns_;i++) 
2260        reducedCost_[i] *= direction;
2261      for (i=0;i<numberRows_;i++) 
2262        dual_[i] *= direction;
2263    }
2264    // row reduced costs
2265    if (!dj_) {
2266      dj_ = new double[numberRows2+numberColumns_];
2267      reducedCostWork_ = dj_;
2268      rowReducedCost_ = dj_+numberColumns_;
2269      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
2270      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
2271    }
2272    if (!solution_||(what&32)!=0) {
2273      if (!solution_)
2274        solution_ = new double[numberRows2+numberColumns_];
2275      columnActivityWork_ = solution_;
2276      rowActivityWork_ = solution_+numberColumns_;
2277      if (status_) {
2278        for (i=0;i<numberColumns_;i++) {
2279          columnActivityWork_[i] = columnActivity_[i];
2280          setFakeBound(i,noFake);
2281        }
2282        for (i=0;i<numberRows_;i++) {
2283          rowActivityWork_[i] = rowActivity_[i];
2284          setFakeBound(i+numberColumns_,noFake);
2285        }
2286      } else {
2287        for (i=0;i<numberColumns_;i++) {
2288          columnActivityWork_[i] = columnActivity_[i];
2289        }
2290        for (i=0;i<numberRows_;i++) {
2291          rowActivityWork_[i] = rowActivity_[i];
2292        }
2293      }
2294      for (i=0;i<numberColumns_+numberRows_;i++) {
2295        if (fabs(solution_[i])>1.0e20) {
2296          //printf("bad value of %g for column %d\n",solution_[i],i);
2297          setColumnStatus(i,superBasic);
2298          if (upper_[i]<0.0) {
2299            solution_[i]=upper_[i];
2300          } else if (lower_[i]>0.0) {
2301            solution_[i]=lower_[i];
2302          } else {
2303            solution_[i]=0.0;
2304          }
2305        }
2306      }
2307    }
2308  }
2309  if ((what&16)!=0) {
2310    //check matrix
2311    if (!matrix_)
2312      matrix_=new ClpPackedMatrix();
2313    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20)) {
2314      problemStatus_=4;
2315      goodMatrix= false;
2316    }
2317    if (makeRowCopy) {
2318      delete rowCopy_;
2319      // may return NULL if can't give row copy
2320      rowCopy_ = matrix_->reverseOrderedCopy();
2321#ifdef TAKEOUT
2322      {
2323
2324        ClpPackedMatrix* rowCopy =
2325          dynamic_cast< ClpPackedMatrix*>(rowCopy_);
2326        const int * column = rowCopy->getIndices();
2327        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2328        const double * element = rowCopy->getElements();
2329        int i;
2330        for (i=133;i<numberRows_;i++) {
2331          if (rowStart[i+1]-rowStart[i]==10||rowStart[i+1]-rowStart[i]==15)
2332            printf("Row %d has %d elements\n",i,rowStart[i+1]-rowStart[i]);
2333        }
2334      } 
2335#endif
2336    }
2337  }
2338  if ((what&4)!=0) {
2339    delete [] cost_;
2340    // extra copy with original costs
2341    int nTotal = numberRows2+numberColumns_;
2342    //cost_ = new double[2*nTotal];
2343    cost_ = new double[nTotal];
2344    objectiveWork_ = cost_;
2345    rowObjectiveWork_ = cost_+numberColumns_;
2346    memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double));
2347    if (rowObjective_)
2348      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
2349    else
2350      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2351    // and initialize changes to zero
2352    //memset(cost_+nTotal,0,nTotal*sizeof(double));
2353  }
2354  // do scaling if needed
2355  if (scalingFlag_>0&&!rowScale_) {
2356    if (matrix_->scale(this))
2357      scalingFlag_=-scalingFlag_; // not scaled after all
2358    if (rowScale_&&automaticScale_) {
2359      // try automatic scaling
2360      double smallestObj=1.0e100;
2361      double largestObj=0.0;
2362      double largestRhs=0.0;
2363      for (i=0;i<numberColumns_+numberRows_;i++) {
2364        double value = fabs(cost_[i]);
2365        if (i<numberColumns_)
2366          value *= columnScale_[i];
2367        if (value&&lower_[i]!=upper_[i]) {
2368          smallestObj = CoinMin(smallestObj,value);
2369          largestObj = CoinMax(largestObj,value);
2370        }
2371        if (lower_[i]>0.0||upper_[i]<0.0) {
2372          double scale;
2373          if (i<numberColumns_)
2374            scale = 1.0/columnScale_[i];
2375          else
2376            scale = rowScale_[i-numberColumns_];
2377          //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2378          if (lower_[i]>0)
2379            largestRhs=CoinMax(largestRhs,lower_[i]*scale);
2380          if (upper_[i]<0.0)
2381            largestRhs=CoinMax(largestRhs,-upper_[i]*scale);
2382        }
2383      }
2384      printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2385      bool scalingDone=false;
2386      // look at element range
2387      double smallestNegative;
2388      double largestNegative;
2389      double smallestPositive;
2390      double largestPositive;
2391      matrix_->rangeOfElements(smallestNegative, largestNegative,
2392                               smallestPositive, largestPositive);
2393      smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2394      largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2395      if (largestObj) {
2396        double ratio = largestObj/smallestObj;
2397        double scale=1.0;
2398        if (ratio<1.0e8) {
2399          // reasonable
2400          if (smallestObj<1.0e-4) {
2401            // may as well scale up
2402            scalingDone=true;
2403            scale=1.0e-3/smallestObj;
2404          } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2405            //done=true;
2406          } else {
2407            scalingDone=true;
2408            if (algorithm_<0) {
2409              scale = 1.0e6/largestObj;
2410            } else {
2411              scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2412            }
2413          }
2414        } else if (ratio<1.0e12) {
2415          // not so good
2416          if (smallestObj<1.0e-7) {
2417            // may as well scale up
2418            scalingDone=true;
2419            scale=1.0e-6/smallestObj;
2420          } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2421            //done=true;
2422          } else {
2423            scalingDone=true;
2424            if (algorithm_<0) {
2425              scale = 1.0e7/largestObj;
2426            } else {
2427              scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2428            }
2429          }
2430        } else {
2431          // Really nasty problem
2432          if (smallestObj<1.0e-8) {
2433            // may as well scale up
2434            scalingDone=true;
2435            scale=1.0e-7/smallestObj;
2436            largestObj *= scale;
2437          } 
2438          if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2439            //done=true;
2440          } else {
2441            scalingDone=true;
2442            if (algorithm_<0) {
2443              scale = 1.0e7/largestObj;
2444            } else {
2445              scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2446            }
2447          }
2448        }
2449        objectiveScale_=scale;
2450      }
2451      if (largestRhs>1.0e12) {
2452        scalingDone=true;
2453        rhsScale_=1.0e9/largestRhs;
2454      } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
2455        scalingDone=true;
2456        rhsScale_=1.0e6/largestRhs;
2457      } else {
2458        rhsScale_=1.0;
2459      }
2460      if (scalingDone) {
2461        handler_->message(CLP_RIM_SCALE,messages_)
2462          <<objectiveScale_<<rhsScale_
2463          <<CoinMessageEol;
2464      }
2465    }
2466  } else if (makeRowCopy&&(what&16)!=0&&scalingFlag_>0) {
2467    matrix_->scaleRowCopy(this);
2468  }
2469  if ((what&4)!=0) {
2470    double direction = optimizationDirection_*objectiveScale_;
2471    // but also scale by scale factors
2472    // not really sure about row scaling
2473    if (!rowScale_) {
2474      if (direction!=1.0) {
2475        for (i=0;i<numberRows_;i++)
2476          rowObjectiveWork_[i] *= direction;
2477        for (i=0;i<numberColumns_;i++)
2478          objectiveWork_[i] *= direction;
2479      }
2480    } else {
2481      for (i=0;i<numberRows_;i++)
2482        rowObjectiveWork_[i] *= direction/rowScale_[i];
2483      for (i=0;i<numberColumns_;i++)
2484        objectiveWork_[i] *= direction*columnScale_[i];
2485    }
2486  }
2487  if ((what&(1+32))!=0) {
2488    if(rowScale_) {
2489      for (i=0;i<numberColumns_;i++) {
2490        double multiplier = rhsScale_/columnScale_[i];
2491        if (columnLowerWork_[i]>-1.0e50)
2492          columnLowerWork_[i] *= multiplier;
2493        if (columnUpperWork_[i]<1.0e50)
2494          columnUpperWork_[i] *= multiplier;
2495       
2496      }
2497      for (i=0;i<numberRows_;i++) {
2498        double multiplier = rhsScale_*rowScale_[i];
2499        if (rowLowerWork_[i]>-1.0e50)
2500          rowLowerWork_[i] *= multiplier;
2501        if (rowUpperWork_[i]<1.0e50)
2502          rowUpperWork_[i] *= multiplier;
2503      }
2504    } else if (rhsScale_!=1.0) {
2505      for (i=0;i<numberColumns_+numberRows_;i++) {
2506        if (lower_[i]>-1.0e50)
2507          lower_[i] *= rhsScale_;
2508        if (upper_[i]<1.0e50)
2509          upper_[i] *= rhsScale_;
2510       
2511      }
2512    }
2513  }
2514  if ((what&(8+32))!=0) {
2515    if (rowScale_) {
2516      // on entry
2517      for (i=0;i<numberColumns_;i++) {
2518        columnActivityWork_[i] /= columnScale_[i];
2519        columnActivityWork_[i] *= rhsScale_;
2520        reducedCostWork_[i] *= objectiveScale_*columnScale_[i];
2521      }
2522      for (i=0;i<numberRows_;i++) {
2523        rowActivityWork_[i] *= rowScale_[i]*rhsScale_;
2524        dual_[i] /= rowScale_[i];
2525        dual_[i] *= objectiveScale_;
2526        rowReducedCost_[i] = dual_[i];
2527      }
2528    } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
2529      // on entry
2530      for (i=0;i<numberColumns_;i++) {
2531        columnActivityWork_[i] *= rhsScale_;
2532        reducedCostWork_[i] *= objectiveScale_;
2533      }
2534      for (i=0;i<numberRows_;i++) {
2535        rowActivityWork_[i] *= rhsScale_;
2536        dual_[i] *= objectiveScale_;
2537        rowReducedCost_[i] = dual_[i];
2538      }
2539    }
2540  }
2541 
2542  if ((what&16)!=0) {
2543    // check rim of problem okay
2544    if (!sanityCheck())
2545      goodMatrix=false;
2546  } 
2547  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
2548  // maybe we need to move scales to SimplexModel for factorization?
2549  if ((what&8)!=0&&!pivotVariable_) {
2550    pivotVariable_=new int[numberRows2];
2551    for (int i=0;i<numberRows2;i++)
2552      pivotVariable_[i]=-1;
2553  }
2554
2555  if ((what&16)!=0&&!rowArray_[2]) {
2556    // get some arrays
2557    int iRow,iColumn;
2558    // these are "indexed" arrays so we always know where nonzeros are
2559    /**********************************************************
2560      rowArray_[3] is long enough for rows+columns
2561    *********************************************************/
2562    for (iRow=0;iRow<4;iRow++) {
2563      if (!rowArray_[iRow]) {
2564        rowArray_[iRow]=new CoinIndexedVector();
2565        int length =numberRows2+factorization_->maximumPivots();
2566        if (iRow==3||objective_->type()>1)
2567          length += numberColumns_;
2568        rowArray_[iRow]->reserve(length);
2569      }
2570    }
2571   
2572    for (iColumn=0;iColumn<2;iColumn++) {
2573      if (!columnArray_[iColumn]) {
2574        columnArray_[iColumn]=new CoinIndexedVector();
2575        if (!iColumn)
2576          columnArray_[iColumn]->reserve(numberColumns_);
2577        else
2578          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
2579      }
2580    }   
2581  }
2582  double primalTolerance=dblParam_[ClpPrimalTolerance];
2583  if ((what&1)!=0) {
2584    // fix any variables with tiny gaps
2585    for (i=0;i<numberColumns_;i++) {
2586      if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2587        if (columnLowerWork_[i]>=0.0) {
2588          columnUpperWork_[i] = columnLowerWork_[i];
2589        } else if (columnUpperWork_[i]<=0.0) {
2590          columnLowerWork_[i] = columnUpperWork_[i];
2591        } else {
2592          columnUpperWork_[i] = 0.0;
2593          columnLowerWork_[i] = 0.0;
2594        }
2595      }
2596    }
2597    for (i=0;i<numberRows_;i++) {
2598      if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2599        if (rowLowerWork_[i]>=0.0) {
2600          rowUpperWork_[i] = rowLowerWork_[i];
2601        } else if (rowUpperWork_[i]<=0.0) {
2602          rowLowerWork_[i] = rowUpperWork_[i];
2603        } else {
2604          rowUpperWork_[i] = 0.0;
2605          rowLowerWork_[i] = 0.0;
2606        }
2607      }
2608    }
2609  }
2610  if (problemStatus_==10) {
2611    problemStatus_=-1;
2612    handler_->setLogLevel(saveLevel); // switch back messages
2613  }
2614  if ((what&5)!=0) 
2615    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
2616  return goodMatrix;
2617}
2618void
2619ClpSimplex::deleteRim(int getRidOfFactorizationData)
2620{
2621  int i;
2622  if (problemStatus_!=1&&problemStatus_!=2) {
2623    delete [] ray_;
2624    ray_=NULL;
2625  }
2626  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
2627  upperOut_=1.0;
2628  // ray may be null if in branch and bound
2629  if (rowScale_) {
2630    // Collect infeasibilities
2631    int numberPrimalScaled=0;
2632    int numberPrimalUnscaled=0;
2633    int numberDualScaled=0;
2634    int numberDualUnscaled=0;
2635    double scaleC = 1.0/objectiveScale_;
2636    double scaleR = 1.0/rhsScale_;
2637    for (i=0;i<numberColumns_;i++) {
2638      double scaleFactor = columnScale_[i];
2639      double valueScaled = columnActivityWork_[i];
2640      double lowerScaled = columnLowerWork_[i];
2641      double upperScaled = columnUpperWork_[i];
2642      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2643        if (valueScaled<lowerScaled-primalTolerance_||
2644            valueScaled>upperScaled+primalTolerance_)
2645          numberPrimalScaled++;
2646        else
2647          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2648      }
2649      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
2650      double value = columnActivity_[i];
2651      if (value<columnLower_[i]-primalTolerance_)
2652        numberPrimalUnscaled++;
2653      else if (value>columnUpper_[i]+primalTolerance_)
2654        numberPrimalUnscaled++;
2655      double valueScaledDual = reducedCostWork_[i];
2656      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2657        numberDualScaled++;
2658      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2659        numberDualScaled++;
2660      reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
2661      double valueDual = reducedCost_[i];
2662      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2663        numberDualUnscaled++;
2664      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2665        numberDualUnscaled++;
2666    }
2667    for (i=0;i<numberRows_;i++) {
2668      double scaleFactor = rowScale_[i];
2669      double valueScaled = rowActivityWork_[i];
2670      double lowerScaled = rowLowerWork_[i];
2671      double upperScaled = rowUpperWork_[i];
2672      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2673        if (valueScaled<lowerScaled-primalTolerance_||
2674            valueScaled>upperScaled+primalTolerance_)
2675          numberPrimalScaled++;
2676        else
2677          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2678      }
2679      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
2680      double value = rowActivity_[i];
2681      if (value<rowLower_[i]-primalTolerance_)
2682        numberPrimalUnscaled++;
2683      else if (value>rowUpper_[i]+primalTolerance_)
2684        numberPrimalUnscaled++;
2685      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2686      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2687        numberDualScaled++;
2688      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2689        numberDualScaled++;
2690      dual_[i] = valueScaledDual*scaleFactor*scaleC;
2691      double valueDual = dual_[i]; 
2692      if (rowObjective_)
2693        valueDual += rowObjective_[i];
2694      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2695        numberDualUnscaled++;
2696      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2697        numberDualUnscaled++;
2698    }
2699    if (!problemStatus_&&!secondaryStatus_) {
2700      // See if we need to set secondary status
2701      if (numberPrimalUnscaled) {
2702        if (numberDualUnscaled) 
2703          secondaryStatus_=4;
2704        else
2705          secondaryStatus_=2;
2706      } else {
2707        if (numberDualUnscaled) 
2708          secondaryStatus_=3;
2709      }
2710    }
2711    if (problemStatus_==2) {
2712      for (i=0;i<numberColumns_;i++) {
2713        ray_[i] *= columnScale_[i];
2714      }
2715    } else if (problemStatus_==1&&ray_) {
2716      for (i=0;i<numberRows_;i++) {
2717        ray_[i] *= rowScale_[i];
2718      }
2719    }
2720  } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
2721    // Collect infeasibilities
2722    int numberPrimalScaled=0;
2723    int numberPrimalUnscaled=0;
2724    int numberDualScaled=0;
2725    int numberDualUnscaled=0;
2726    double scaleC = 1.0/objectiveScale_;
2727    double scaleR = 1.0/rhsScale_;
2728    for (i=0;i<numberColumns_;i++) {
2729      double valueScaled = columnActivityWork_[i];
2730      double lowerScaled = columnLowerWork_[i];
2731      double upperScaled = columnUpperWork_[i];
2732      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2733        if (valueScaled<lowerScaled-primalTolerance_||
2734            valueScaled>upperScaled+primalTolerance_)
2735          numberPrimalScaled++;
2736        else
2737          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2738      }
2739      columnActivity_[i] = valueScaled*scaleR;
2740      double value = columnActivity_[i];
2741      if (value<columnLower_[i]-primalTolerance_)
2742        numberPrimalUnscaled++;
2743      else if (value>columnUpper_[i]+primalTolerance_)
2744        numberPrimalUnscaled++;
2745      double valueScaledDual = reducedCostWork_[i];
2746      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2747        numberDualScaled++;
2748      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2749        numberDualScaled++;
2750      reducedCost_[i] = valueScaledDual*scaleC;
2751      double valueDual = reducedCost_[i];
2752      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2753        numberDualUnscaled++;
2754      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2755        numberDualUnscaled++;
2756    }
2757    for (i=0;i<numberRows_;i++) {
2758      double valueScaled = rowActivityWork_[i];
2759      double lowerScaled = rowLowerWork_[i];
2760      double upperScaled = rowUpperWork_[i];
2761      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2762        if (valueScaled<lowerScaled-primalTolerance_||
2763            valueScaled>upperScaled+primalTolerance_)
2764          numberPrimalScaled++;
2765        else
2766          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2767      }
2768      rowActivity_[i] = valueScaled*scaleR;
2769      double value = rowActivity_[i];
2770      if (value<rowLower_[i]-primalTolerance_)
2771        numberPrimalUnscaled++;
2772      else if (value>rowUpper_[i]+primalTolerance_)
2773        numberPrimalUnscaled++;
2774      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2775      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2776        numberDualScaled++;
2777      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2778        numberDualScaled++;
2779      dual_[i] = valueScaledDual*scaleC;
2780      double valueDual = dual_[i]; 
2781      if (rowObjective_)
2782        valueDual += rowObjective_[i];
2783      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2784        numberDualUnscaled++;
2785      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2786        numberDualUnscaled++;
2787    }
2788    if (!problemStatus_&&!secondaryStatus_) {
2789      // See if we need to set secondary status
2790      if (numberPrimalUnscaled) {
2791        if (numberDualUnscaled) 
2792          secondaryStatus_=4;
2793        else
2794          secondaryStatus_=2;
2795      } else {
2796        if (numberDualUnscaled) 
2797          secondaryStatus_=3;
2798      }
2799    }
2800  } else {
2801    if (columnActivityWork_) {
2802      for (i=0;i<numberColumns_;i++) {
2803        double value = columnActivityWork_[i];
2804        double lower = columnLowerWork_[i];
2805        double upper = columnUpperWork_[i];
2806        if (lower>-1.0e20||upper<1.0e20) {
2807          if (value>lower&&value<upper)
2808            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2809        }
2810        columnActivity_[i] = columnActivityWork_[i];
2811        reducedCost_[i] = reducedCostWork_[i];
2812      }
2813      for (i=0;i<numberRows_;i++) {
2814        double value = rowActivityWork_[i];
2815        double lower = rowLowerWork_[i];
2816        double upper = rowUpperWork_[i];
2817        if (lower>-1.0e20||upper<1.0e20) {
2818          if (value>lower&&value<upper)
2819            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2820        }
2821        rowActivity_[i] = rowActivityWork_[i];
2822      }
2823    }
2824  }
2825  // switch off scalefactor if auto
2826  if (automaticScale_) {
2827    rhsScale_=1.0;
2828    objectiveScale_=1.0;
2829  }
2830  if (optimizationDirection_!=1.0) {
2831    // and modify all dual signs
2832    for (i=0;i<numberColumns_;i++) 
2833      reducedCost_[i] *= optimizationDirection_;
2834    for (i=0;i<numberRows_;i++) 
2835      dual_[i] *= optimizationDirection_;
2836  }
2837  // scaling may have been turned off
2838  scalingFlag_ = abs(scalingFlag_);
2839  if(getRidOfFactorizationData>=0)
2840    gutsOfDelete(getRidOfFactorizationData+1);
2841  // get rid of data
2842  matrix_->generalExpanded(this,13,scalingFlag_);
2843}
2844void 
2845ClpSimplex::setDualBound(double value)
2846{
2847  if (value>0.0)
2848    dualBound_=value;
2849}
2850void 
2851ClpSimplex::setInfeasibilityCost(double value)
2852{
2853  if (value>0.0)
2854    infeasibilityCost_=value;
2855}
2856void ClpSimplex::setNumberRefinements( int value) 
2857{
2858  if (value>=0&&value<10)
2859    numberRefinements_=value;
2860}
2861// Sets row pivot choice algorithm in dual
2862void 
2863ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
2864{
2865  delete dualRowPivot_;
2866  dualRowPivot_ = choice.clone(true);
2867}
2868// Sets row pivot choice algorithm in dual
2869void 
2870ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
2871{
2872  delete primalColumnPivot_;
2873  primalColumnPivot_ = choice.clone(true);
2874}
2875// Passes in factorization
2876void 
2877ClpSimplex::setFactorization( ClpFactorization & factorization)
2878{
2879  delete factorization_;
2880  factorization_= new ClpFactorization(factorization);
2881}
2882/* Perturbation:
2883   -50 to +50 - perturb by this power of ten (-6 sounds good)
2884   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
2885   101 - we are perturbed
2886   102 - don't try perturbing again
2887   default is 100
2888*/
2889void 
2890ClpSimplex::setPerturbation(int value)
2891{
2892  if(value<=100&&value >=-1000) {
2893    perturbation_=value;
2894  } 
2895}
2896// Sparsity on or off
2897bool 
2898ClpSimplex::sparseFactorization() const
2899{
2900  return factorization_->sparseThreshold()!=0;
2901}
2902void 
2903ClpSimplex::setSparseFactorization(bool value)
2904{
2905  if (value) {
2906    if (!factorization_->sparseThreshold())
2907      factorization_->goSparse();
2908  } else {
2909    factorization_->sparseThreshold(0);
2910  }
2911}
2912void checkCorrect(ClpSimplex * model,int iRow,
2913                  const double * element,const int * rowStart,const int * rowLength,
2914                  const int * column,
2915                  const double * columnLower_, const double * columnUpper_,
2916                  int infiniteUpperC,
2917                  int infiniteLowerC,
2918                  double &maximumUpC,
2919                  double &maximumDownC)
2920{
2921  int infiniteUpper = 0;
2922  int infiniteLower = 0;
2923  double maximumUp = 0.0;
2924  double maximumDown = 0.0;
2925  CoinBigIndex rStart = rowStart[iRow];
2926  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2927  CoinBigIndex j;
2928  double large=1.0e15;
2929  int iColumn;
2930  // Compute possible lower and upper ranges
2931 
2932  for (j = rStart; j < rEnd; ++j) {
2933    double value=element[j];
2934    iColumn = column[j];
2935    if (value > 0.0) {
2936      if (columnUpper_[iColumn] >= large) {
2937        ++infiniteUpper;
2938      } else {
2939        maximumUp += columnUpper_[iColumn] * value;
2940      }
2941      if (columnLower_[iColumn] <= -large) {
2942        ++infiniteLower;
2943      } else {
2944        maximumDown += columnLower_[iColumn] * value;
2945      }
2946    } else if (value<0.0) {
2947      if (columnUpper_[iColumn] >= large) {
2948        ++infiniteLower;
2949      } else {
2950        maximumDown += columnUpper_[iColumn] * value;
2951      }
2952      if (columnLower_[iColumn] <= -large) {
2953        ++infiniteUpper;
2954      } else {
2955        maximumUp += columnLower_[iColumn] * value;
2956      }
2957    }
2958  }
2959  assert (infiniteLowerC==infiniteLower);
2960  assert (infiniteUpperC==infiniteUpper);
2961  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
2962    printf("row %d comp up %g, true up %g\n",iRow,
2963           maximumUpC,maximumUp);
2964  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
2965    printf("row %d comp down %g, true down %g\n",iRow,
2966           maximumDownC,maximumDown);
2967  maximumUpC=maximumUp;
2968  maximumDownC=maximumDown;
2969}
2970
2971/* Tightens primal bounds to make dual faster.  Unless
2972   fixed, bounds are slightly looser than they could be.
2973   This is to make dual go faster and is probably not needed
2974   with a presolve.  Returns non-zero if problem infeasible
2975
2976   Fudge for branch and bound - put bounds on columns of factor *
2977   largest value (at continuous) - should improve stability
2978   in branch and bound on infeasible branches (0.0 is off)
2979*/
2980int 
2981ClpSimplex::tightenPrimalBounds(double factor)
2982{
2983 
2984  // Get a row copy in standard format
2985  CoinPackedMatrix copy;
2986  copy.reverseOrderedCopyOf(*matrix());
2987  // get matrix data pointers
2988  const int * column = copy.getIndices();
2989  const CoinBigIndex * rowStart = copy.getVectorStarts();
2990  const int * rowLength = copy.getVectorLengths(); 
2991  const double * element = copy.getElements();
2992  int numberChanged=1,iPass=0;
2993  double large = largeValue(); // treat bounds > this as infinite
2994#ifndef NDEBUG
2995  double large2= 1.0e10*large;
2996#endif
2997  int numberInfeasible=0;
2998  int totalTightened = 0;
2999
3000  double tolerance = primalTolerance();
3001
3002
3003  // Save column bounds
3004  double * saveLower = new double [numberColumns_];
3005  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
3006  double * saveUpper = new double [numberColumns_];
3007  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
3008
3009  int iRow, iColumn;
3010
3011  // If wanted - tighten column bounds using solution
3012  if (factor) {
3013    assert (factor>1.0);
3014    double largest=0.0;
3015    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3016      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3017        largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
3018      }
3019    }
3020    largest *= factor;
3021    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3022      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3023        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
3024        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
3025      }
3026    }
3027  }
3028#define MAXPASS 10
3029
3030  // Loop round seeing if we can tighten bounds
3031  // Would be faster to have a stack of possible rows
3032  // and we put altered rows back on stack
3033  int numberCheck=-1;
3034  while(numberChanged>numberCheck) {
3035
3036    numberChanged = 0; // Bounds tightened this pass
3037   
3038    if (iPass==MAXPASS) break;
3039    iPass++;
3040   
3041    for (iRow = 0; iRow < numberRows_; iRow++) {
3042
3043      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
3044
3045        // possible row
3046        int infiniteUpper = 0;
3047        int infiniteLower = 0;
3048        double maximumUp = 0.0;
3049        double maximumDown = 0.0;
3050        double newBound;
3051        CoinBigIndex rStart = rowStart[iRow];
3052        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
3053        CoinBigIndex j;
3054        // Compute possible lower and upper ranges
3055     
3056        for (j = rStart; j < rEnd; ++j) {
3057          double value=element[j];
3058          iColumn = column[j];
3059          if (value > 0.0) {
3060            if (columnUpper_[iColumn] >= large) {
3061              ++infiniteUpper;
3062            } else {
3063              maximumUp += columnUpper_[iColumn] * value;
3064            }
3065            if (columnLower_[iColumn] <= -large) {
3066              ++infiniteLower;
3067            } else {
3068              maximumDown += columnLower_[iColumn] * value;
3069            }
3070          } else if (value<0.0) {
3071            if (columnUpper_[iColumn] >= large) {
3072              ++infiniteLower;
3073            } else {
3074              maximumDown += columnUpper_[iColumn] * value;
3075            }
3076            if (columnLower_[iColumn] <= -large) {
3077              ++infiniteUpper;
3078            } else {
3079              maximumUp += columnLower_[iColumn] * value;
3080            }
3081          }
3082        }
3083        // Build in a margin of error
3084        maximumUp += 1.0e-8*fabs(maximumUp);
3085        maximumDown -= 1.0e-8*fabs(maximumDown);
3086        double maxUp = maximumUp+infiniteUpper*1.0e31;
3087        double maxDown = maximumDown-infiniteLower*1.0e31;
3088        if (maxUp <= rowUpper_[iRow] + tolerance && 
3089            maxDown >= rowLower_[iRow] - tolerance) {
3090         
3091          // Row is redundant - make totally free
3092          // NO - can't do this for postsolve
3093          // rowLower_[iRow]=-COIN_DBL_MAX;
3094          // rowUpper_[iRow]=COIN_DBL_MAX;
3095          //printf("Redundant row in presolveX %d\n",iRow);
3096
3097        } else {
3098          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
3099              maxDown > rowUpper_[iRow]+100.0*tolerance) {
3100            // problem is infeasible - exit at once
3101            numberInfeasible++;
3102            break;
3103          }
3104          double lower = rowLower_[iRow];
3105          double upper = rowUpper_[iRow];
3106          for (j = rStart; j < rEnd; ++j) {
3107            double value=element[j];
3108            iColumn = column[j];
3109            double nowLower = columnLower_[iColumn];
3110            double nowUpper = columnUpper_[iColumn];
3111            if (value > 0.0) {
3112              // positive value
3113              if (lower>-large) {
3114                if (!infiniteUpper) {
3115                  assert(nowUpper < large2);
3116                  newBound = nowUpper + 
3117                    (lower - maximumUp) / value;
3118                  // relax if original was large
3119                  if (fabs(maximumUp)>1.0e8)
3120                    newBound -= 1.0e-12*fabs(maximumUp);
3121                } else if (infiniteUpper==1&&nowUpper>large) {
3122                  newBound = (lower -maximumUp) / value;
3123                  // relax if original was large
3124                  if (fabs(maximumUp)>1.0e8)
3125                    newBound -= 1.0e-12*fabs(maximumUp);
3126                } else {
3127                  newBound = -COIN_DBL_MAX;
3128                }
3129                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3130                  // Tighten the lower bound
3131                  columnLower_[iColumn] = newBound;
3132                  numberChanged++;
3133                  // check infeasible (relaxed)
3134                  if (nowUpper - newBound < 
3135                      -100.0*tolerance) {
3136                    numberInfeasible++;
3137                  }
3138                  // adjust
3139                  double now;
3140                  if (nowLower<-large) {
3141                    now=0.0;
3142                    infiniteLower--;
3143                  } else {
3144                    now = nowLower;
3145                  }
3146                  maximumDown += (newBound-now) * value;
3147                  nowLower = newBound;
3148#ifdef DEBUG
3149                  checkCorrect(this,iRow,
3150                               element, rowStart, rowLength,
3151                               column,
3152                               columnLower_,  columnUpper_,
3153                               infiniteUpper,
3154                               infiniteLower,
3155                               maximumUp,
3156                               maximumDown);
3157#endif
3158                }
3159              } 
3160              if (upper <large) {
3161                if (!infiniteLower) {
3162                  assert(nowLower >- large2);
3163                  newBound = nowLower + 
3164                    (upper - maximumDown) / value;
3165                  // relax if original was large
3166                  if (fabs(maximumDown)>1.0e8)
3167                    newBound += 1.0e-12*fabs(maximumDown);
3168                } else if (infiniteLower==1&&nowLower<-large) {
3169                  newBound =   (upper - maximumDown) / value;
3170                  // relax if original was large
3171                  if (fabs(maximumDown)>1.0e8)
3172                    newBound += 1.0e-12*fabs(maximumDown);
3173                } else {
3174                  newBound = COIN_DBL_MAX;
3175                }
3176                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3177                  // Tighten the upper bound
3178                  columnUpper_[iColumn] = newBound;
3179                  numberChanged++;
3180                  // check infeasible (relaxed)
3181                  if (newBound - nowLower < 
3182                      -100.0*tolerance) {
3183                    numberInfeasible++;
3184                  }
3185                  // adjust
3186                  double now;
3187                  if (nowUpper>large) {
3188                    now=0.0;
3189                    infiniteUpper--;
3190                  } else {
3191                    now = nowUpper;
3192                  }
3193                  maximumUp += (newBound-now) * value;
3194                  nowUpper = newBound;
3195#ifdef DEBUG
3196                  checkCorrect(this,iRow,
3197                               element, rowStart, rowLength,
3198                               column,
3199                               columnLower_,  columnUpper_,
3200                               infiniteUpper,
3201                               infiniteLower,
3202                               maximumUp,
3203                               maximumDown);
3204#endif
3205                }
3206              }
3207            } else {
3208              // negative value
3209              if (lower>-large) {
3210                if (!infiniteUpper) {
3211                  assert(nowLower < large2);
3212                  newBound = nowLower + 
3213                    (lower - maximumUp) / value;
3214                  // relax if original was large
3215                  if (fabs(maximumUp)>1.0e8)
3216                    newBound += 1.0e-12*fabs(maximumUp);
3217                } else if (infiniteUpper==1&&nowLower<-large) {
3218                  newBound = (lower -maximumUp) / value;
3219                  // relax if original was large
3220                  if (fabs(maximumUp)>1.0e8)
3221                    newBound += 1.0e-12*fabs(maximumUp);
3222                } else {
3223                  newBound = COIN_DBL_MAX;
3224                }
3225                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3226                  // Tighten the upper bound
3227                  columnUpper_[iColumn] = newBound;
3228                  numberChanged++;
3229                  // check infeasible (relaxed)
3230                  if (newBound - nowLower < 
3231                      -100.0*tolerance) {
3232                    numberInfeasible++;
3233                  }
3234                  // adjust
3235                  double now;
3236                  if (nowUpper>large) {
3237                    now=0.0;
3238                    infiniteLower--;
3239                  } else {
3240                    now = nowUpper;
3241                  }
3242                  maximumDown += (newBound-now) * value;
3243                  nowUpper = newBound;
3244#ifdef DEBUG
3245                  checkCorrect(this,iRow,
3246                               element, rowStart, rowLength,
3247                               column,
3248                               columnLower_,  columnUpper_,
3249                               infiniteUpper,
3250                               infiniteLower,
3251                               maximumUp,
3252                               maximumDown);
3253#endif
3254                }
3255              }
3256              if (upper <large) {
3257                if (!infiniteLower) {
3258                  assert(nowUpper < large2);
3259                  newBound = nowUpper + 
3260                    (upper - maximumDown) / value;
3261                  // relax if original was large
3262                  if (fabs(maximumDown)>1.0e8)
3263                    newBound -= 1.0e-12*fabs(maximumDown);
3264                } else if (infiniteLower==1&&nowUpper>large) {
3265                  newBound =   (upper - maximumDown) / value;
3266                  // relax if original was large
3267                  if (fabs(maximumDown)>1.0e8)
3268                    newBound -= 1.0e-12*fabs(maximumDown);
3269                } else {
3270                  newBound = -COIN_DBL_MAX;
3271                }
3272                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3273                  // Tighten the lower bound
3274                  columnLower_[iColumn] = newBound;
3275                  numberChanged++;
3276                  // check infeasible (relaxed)
3277                  if (nowUpper - newBound < 
3278                      -100.0*tolerance) {
3279                    numberInfeasible++;
3280                  }
3281                  // adjust
3282                  double now;
3283                  if (nowLower<-large) {
3284                    now=0.0;
3285                    infiniteUpper--;
3286                  } else {
3287                    now = nowLower;
3288                  }
3289                  maximumUp += (newBound-now) * value;
3290                  nowLower = newBound;
3291#ifdef DEBUG
3292                  checkCorrect(this,iRow,
3293                               element, rowStart, rowLength,
3294                               column,
3295                               columnLower_,  columnUpper_,
3296                               infiniteUpper,
3297                               infiniteLower,
3298                               maximumUp,
3299                               maximumDown);
3300#endif
3301                }
3302              }
3303            }
3304          }
3305        }
3306      }
3307    }
3308    totalTightened += numberChanged;
3309    if (iPass==1)
3310      numberCheck=numberChanged>>4;
3311    if (numberInfeasible) break;
3312  }
3313  if (!numberInfeasible) {
3314    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
3315      <<totalTightened
3316      <<CoinMessageEol;
3317    // Set bounds slightly loose
3318    double useTolerance = 1.0e-3;
3319    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3320      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
3321        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) {
3322          // relax enough so will have correct dj
3323#if 1
3324          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3325                                    columnLower_[iColumn]-100.0*useTolerance);
3326          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3327                                    columnUpper_[iColumn]+100.0*useTolerance);
3328#else
3329          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
3330            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
3331              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
3332            } else {
3333              columnLower_[iColumn]=saveLower[iColumn];
3334              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3335                                        saveLower[iColumn]+100.0*useTolerance);
3336            }
3337          } else {
3338            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
3339              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
3340            } else {
3341              columnUpper_[iColumn]=saveUpper[iColumn];
3342              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3343                                        saveUpper[iColumn]-100.0*useTolerance);
3344            }
3345          }
3346#endif
3347        } else {
3348          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
3349            // relax a bit
3350            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
3351                                        saveUpper[iColumn]);
3352          }
3353          if (columnLower_[iColumn]>saveLower[iColumn]) {
3354            // relax a bit
3355            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
3356                                        saveLower[iColumn]);
3357          }
3358        }
3359      }
3360    }
3361  } else {
3362    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
3363      <<numberInfeasible
3364      <<CoinMessageEol;
3365    // restore column bounds
3366    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
3367    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
3368  }
3369  delete [] saveLower;
3370  delete [] saveUpper;
3371  return (numberInfeasible);
3372}
3373// dual
3374#include "ClpSimplexDual.hpp"
3375#include "ClpSimplexPrimal.hpp"
3376int ClpSimplex::dual (int ifValuesPass )
3377{
3378  int saveQuadraticActivated = objective_->activated();
3379  objective_->setActivated(0);
3380  assert (ifValuesPass>=0&&ifValuesPass<3);
3381  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3382      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3383      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3384      additional data and have no destructor or (non-default) constructor.
3385
3386      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3387      in primal and dual.
3388
3389      As far as I can see this is perfectly safe.
3390  */
3391  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass);
3392  if (problemStatus_==10) {
3393    //printf("Cleaning up with primal\n");
3394    int savePerturbation = perturbation_;
3395    perturbation_=100;
3396    bool denseFactorization = initialDenseFactorization();
3397    // It will be safe to allow dense
3398    setInitialDenseFactorization(true);
3399    // Allow for catastrophe
3400    int saveMax = intParam_[ClpMaxNumIteration];
3401    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
3402      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
3403    // check which algorithms allowed
3404    int dummy;
3405    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
3406      returnCode = ((ClpSimplexPrimal *) this)->primal(1);
3407    else
3408      returnCode = ((ClpSimplexDual *) this)->dual(0);
3409    if (problemStatus_==3&&numberIterations_<saveMax) {
3410      printf("looks like trouble - too many iterations in clean up - trying again\n");
3411      // flatten solution and try again
3412      int iRow,iColumn;
3413      for (iRow=0;iRow<numberRows_;iRow++) {
3414        if (getRowStatus(iRow)!=basic) {
3415          setRowStatus(iRow,superBasic);
3416          // but put to bound if close
3417          if (fabs(rowActivity_[iRow]-rowLower_[iRow])
3418              <=primalTolerance_) {
3419            rowActivity_[iRow]=rowLower_[iRow];
3420            setRowStatus(iRow,atLowerBound);
3421          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])
3422                     <=primalTolerance_) {
3423            rowActivity_[iRow]=rowUpper_[iRow];
3424            setRowStatus(iRow,atUpperBound);
3425          }
3426        }
3427      }
3428      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3429        if (getColumnStatus(iColumn)!=basic) {
3430          setColumnStatus(iColumn,superBasic);
3431          // but put to bound if close
3432          if (fabs(columnActivity_[iColumn]-columnLower_[iColumn])
3433              <=primalTolerance_) {
3434            columnActivity_[iColumn]=columnLower_[iColumn];
3435            setColumnStatus(iColumn,atLowerBound);
3436          } else if (fabs(columnActivity_[iColumn]
3437                          -columnUpper_[iColumn])
3438                     <=primalTolerance_) {
3439            columnActivity_[iColumn]=columnUpper_[iColumn];
3440            setColumnStatus(iColumn,atUpperBound);
3441          }
3442        }
3443      }
3444      problemStatus_=-1;
3445      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 
3446                                          2*numberRows_+numberColumns_,saveMax);
3447      perturbation_=savePerturbation;
3448      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
3449      if (problemStatus_==3&&numberIterations_<saveMax) 
3450        printf("looks like real trouble - too many iterations in second clean up - giving up\n");
3451    }
3452    intParam_[ClpMaxNumIteration] = saveMax;
3453
3454    setInitialDenseFactorization(denseFactorization);
3455    perturbation_=savePerturbation;
3456    if (problemStatus_==10) 
3457      problemStatus_=0;
3458  }
3459  objective_->setActivated(saveQuadraticActivated);
3460  return returnCode;
3461}
3462#include "ClpQuadraticObjective.hpp"
3463// primal
3464int ClpSimplex::primal (int ifValuesPass )
3465{
3466  // See if nonlinear
3467  if (objective_->type()>1&&objective_->activated()) 
3468    return reducedGradient();
3469 
3470  assert (ifValuesPass>=0&&ifValuesPass<3);
3471  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3472      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3473      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3474      additional data and have no destructor or (non-default) constructor.
3475
3476      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3477      in primal and dual.
3478
3479      As far as I can see this is perfectly safe.
3480  */
3481  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
3482  if (problemStatus_==10) {
3483    //printf("Cleaning up with dual\n");
3484    int savePerturbation = perturbation_;
3485    perturbation_=100;
3486    bool denseFactorization = initialDenseFactorization();
3487    // It will be safe to allow dense
3488    setInitialDenseFactorization(true);
3489    // check which algorithms allowed
3490    int dummy;
3491    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3492      double saveBound = dualBound_;
3493      // upperOut_ has largest away from bound
3494      dualBound_=CoinMin(2.0*upperOut_,dualBound_);
3495      returnCode = ((ClpSimplexDual *) this)->dual(0);
3496      dualBound_=saveBound;
3497    } else {
3498      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
3499    }
3500    setInitialDenseFactorization(denseFactorization);
3501    perturbation_=savePerturbation;
3502    if (problemStatus_==10) 
3503      problemStatus_=0;
3504  }
3505  return returnCode;
3506}
3507/* Dual ranging.
3508   This computes increase/decrease in cost for each given variable and corresponding
3509   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3510   and numberColumns.. for artificials/slacks.
3511   For non-basic variables the sequence number will be that of the non-basic variables.
3512   
3513   Up to user to provide correct length arrays.
3514   
3515   Returns non-zero if infeasible unbounded etc
3516*/
3517#include "ClpSimplexOther.hpp"
3518int ClpSimplex::dualRanging(int numberCheck,const int * which,
3519                            double * costIncrease, int * sequenceIncrease,
3520                            double * costDecrease, int * sequenceDecrease)
3521{
3522  int savePerturbation = perturbation_;
3523  perturbation_=100;
3524  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3525  if (problemStatus_==10) {
3526    //printf("Cleaning up with dual\n");
3527    bool denseFactorization = initialDenseFactorization();
3528    // It will be safe to allow dense
3529    setInitialDenseFactorization(true);
3530    // check which algorithms allowed
3531    int dummy;
3532    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3533      // upperOut_ has largest away from bound
3534      double saveBound=dualBound_;
3535      if (upperOut_>0.0)
3536        dualBound_=2.0*upperOut_;
3537      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3538      dualBound_=saveBound;
3539    } else {
3540      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3541    }
3542    setInitialDenseFactorization(denseFactorization);
3543    if (problemStatus_==10) 
3544      problemStatus_=0;
3545  }
3546  perturbation_=savePerturbation;
3547  if (problemStatus_) {
3548    finish(); // get rid of arrays
3549    return 1; // odd status
3550  }
3551  ((ClpSimplexOther *) this)->dualRanging(numberCheck,which,
3552                                          costIncrease,sequenceIncrease,
3553                                          costDecrease,sequenceDecrease);
3554  finish(); // get rid of arrays
3555  return 0;
3556}
3557/* Primal ranging.
3558   This computes increase/decrease in value for each given variable and corresponding
3559   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3560   and numberColumns.. for artificials/slacks.
3561   For basic variables the sequence number will be that of the basic variables.
3562   
3563   Up to user to providen correct length arrays.
3564   
3565   Returns non-zero if infeasible unbounded etc
3566*/
3567int ClpSimplex::primalRanging(int numberCheck,const int * which,
3568                  double * valueIncrease, int * sequenceIncrease,
3569                  double * valueDecrease, int * sequenceDecrease)
3570{
3571  int savePerturbation = perturbation_;
3572  perturbation_=100;
3573  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3574  if (problemStatus_==10) {
3575    //printf("Cleaning up with dual\n");
3576    bool denseFactorization = initialDenseFactorization();
3577    // It will be safe to allow dense
3578    setInitialDenseFactorization(true);
3579    // check which algorithms allowed
3580    int dummy;
3581    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3582      // upperOut_ has largest away from bound
3583      double saveBound=dualBound_;
3584      if (upperOut_>0.0)
3585        dualBound_=2.0*upperOut_;
3586      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3587      dualBound_=saveBound;
3588    } else {
3589      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3590    }
3591    setInitialDenseFactorization(denseFactorization);
3592    if (problemStatus_==10) 
3593      problemStatus_=0;
3594  }
3595  perturbation_=savePerturbation;
3596  if (problemStatus_) {
3597    finish(); // get rid of arrays
3598    return 1; // odd status
3599  }
3600  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
3601                                          valueIncrease,sequenceIncrease,
3602                                          valueDecrease,sequenceDecrease);
3603  finish(); // get rid of arrays
3604  return 0;
3605}
3606#include "ClpSimplexNonlinear.hpp"
3607/* Solves nonlinear problem using SLP - may be used as crash
3608   for other algorithms when number of iterations small
3609*/
3610int 
3611ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
3612{
3613  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
3614}
3615// Solves non-linear using reduced gradient
3616int ClpSimplex::reducedGradient(int phase)
3617{
3618  if (objective_->type()<2||!objective_->activated()) {
3619    // no quadratic part
3620    return primal(0);
3621  }
3622  // get feasible
3623  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
3624    objective_->setActivated(0);
3625    double saveDirection = optimizationDirection();
3626    setOptimizationDirection(0.0);
3627    primal(1);
3628    setOptimizationDirection(saveDirection);
3629    objective_->setActivated(1);
3630    // still infeasible
3631    if (numberPrimalInfeasibilities())
3632      return 0;
3633  }
3634  // Now enter method
3635  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
3636  return returnCode;
3637}
3638#include "ClpPredictorCorrector.hpp"
3639#include "ClpCholeskyBase.hpp"
3640// Preference is WSSMP, TAUCS, UFL (just ordering) then base
3641#if WSSMP_BARRIER
3642#include "ClpCholeskyWssmp.hpp"
3643#include "ClpCholeskyWssmpKKT.hpp"
3644#elif UFL_BARRIER
3645#include "ClpCholeskyUfl.hpp"
3646#elif TAUCS_BARRIER
3647#include "ClpCholeskyTaucs.hpp"
3648#endif
3649#include "ClpPresolve.hpp"
3650/* Solves using barrier (assumes you have good cholesky factor code).
3651   Does crossover to simplex if asked*/
3652int 
3653ClpSimplex::barrier(bool crossover)
3654{
3655  ClpSimplex * model2 = this;
3656  int savePerturbation=perturbation_;
3657  ClpInterior barrier;
3658  barrier.borrowModel(*model2);
3659// Preference is WSSMP, UFL, TAUCS then base
3660#ifdef WSSMP_BARRIER
3661  ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
3662  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
3663  //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
3664#elif UFL_BARRIER
3665  ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
3666#elif TAUCS_BARRIER
3667  ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
3668#else
3669  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
3670#endif
3671  barrier.setCholesky(cholesky);
3672  int numberRows = model2->numberRows();
3673  int numberColumns = model2->numberColumns();
3674  int saveMaxIts = model2->maximumIterations();
3675  if (saveMaxIts<1000) {
3676    barrier.setMaximumBarrierIterations(saveMaxIts);
3677    model2->setMaximumIterations(1000000);
3678  }
3679  barrier.primalDual();
3680  int barrierStatus=barrier.status();
3681  double gap = barrier.complementarityGap();
3682  // get which variables are fixed
3683  double * saveLower=NULL;
3684  double * saveUpper=NULL;
3685  ClpPresolve pinfo2;
3686  ClpSimplex * saveModel2=NULL;
3687  int numberFixed = barrier.numberFixed();
3688  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
3689    // may as well do presolve
3690    int numberRows = barrier.numberRows();
3691    int numberColumns = barrier.numberColumns();
3692    int numberTotal = numberRows+numberColumns;
3693    saveLower = new double [numberTotal];
3694    saveUpper = new double [numberTotal];
3695    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
3696    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
3697    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
3698    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
3699    barrier.fixFixed();
3700    saveModel2=model2;
3701  }
3702  barrier.returnModel(*model2);
3703  double * rowPrimal = new double [numberRows];
3704  double * columnPrimal = new double [numberColumns];
3705  double * rowDual = new double [numberRows];
3706  double * columnDual = new double [numberColumns];
3707  // move solutions other way
3708  CoinMemcpyN(model2->primalRowSolution(),
3709              numberRows,rowPrimal);
3710  CoinMemcpyN(model2->dualRowSolution(),
3711              numberRows,rowDual);
3712  CoinMemcpyN(model2->primalColumnSolution(),
3713              numberColumns,columnPrimal);
3714  CoinMemcpyN(model2->dualColumnSolution(),
3715              numberColumns,columnDual);
3716  if (saveModel2) {
3717    // do presolve
3718    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
3719                                   false,5,true);
3720  }
3721  if (barrierStatus<4&&crossover) {
3722    // make sure no status left
3723    model2->createStatus();
3724    // solve
3725    model2->setPerturbation(100);
3726    // throw some into basis
3727    {
3728      int numberRows = model2->numberRows();
3729      int numberColumns = model2->numberColumns();
3730      double * dsort = new double[numberColumns];
3731      int * sort = new int[numberColumns];
3732      int n=0;
3733      const double * columnLower = model2->columnLower();
3734      const double * columnUpper = model2->columnUpper();
3735      const double * primalSolution = model2->primalColumnSolution();
3736      double tolerance = 10.0*primalTolerance_;
3737      int i;
3738      for ( i=0;i<numberRows;i++) 
3739        model2->setRowStatus(i,superBasic);
3740      for ( i=0;i<numberColumns;i++) {
3741        double distance = CoinMin(columnUpper[i]-primalSolution[i],
3742                              primalSolution[i]-columnLower[i]);
3743        if (distance>tolerance) {
3744          dsort[n]=-distance;
3745          sort[n++]=i;
3746          model2->setStatus(i,superBasic);
3747        } else if (distance>primalTolerance_) {
3748          model2->setStatus(i,superBasic);
3749        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
3750          model2->setStatus(i,atLowerBound);
3751        } else {
3752          model2->setStatus(i,atUpperBound);
3753        }
3754      }
3755      CoinSort_2(dsort,dsort+n,sort);
3756      n = CoinMin(numberRows,n);
3757      for ( i=0;i<n;i++) {
3758        int iColumn = sort[i];
3759        model2->setStatus(iColumn,basic);
3760      }
3761      delete [] sort;
3762      delete [] dsort;
3763    }
3764    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
3765      int numberRows = model2->numberRows();
3766      int numberColumns = model2->numberColumns();
3767      // just primal values pass
3768      double saveScale = model2->objectiveScale();
3769      model2->setObjectiveScale(1.0e-3);
3770      model2->primal(2);
3771      model2->setObjectiveScale(saveScale);
3772      // save primal solution and copy back dual
3773      CoinMemcpyN(model2->primalRowSolution(),
3774                  numberRows,rowPrimal);
3775      CoinMemcpyN(rowDual,
3776                  numberRows,model2->dualRowSolution());
3777      CoinMemcpyN(model2->primalColumnSolution(),
3778                  numberColumns,columnPrimal);
3779      CoinMemcpyN(columnDual,
3780                  numberColumns,model2->dualColumnSolution());
3781      //model2->primal(1);
3782      // clean up reduced costs and flag variables
3783      {
3784        double * dj = model2->dualColumnSolution();
3785        double * cost = model2->objective();
3786        double * saveCost = new double[numberColumns];
3787        memcpy(saveCost,cost,numberColumns*sizeof(double));
3788        double * saveLower = new double[numberColumns];
3789        double * lower = model2->columnLower();
3790        memcpy(saveLower,lower,numberColumns*sizeof(double));
3791        double * saveUpper = new double[numberColumns];
3792        double * upper = model2->columnUpper();
3793        memcpy(saveUpper,upper,numberColumns*sizeof(double));
3794        int i;
3795        double tolerance = 10.0*dualTolerance_;
3796        for ( i=0;i<numberColumns;i++) {
3797          if (model2->getStatus(i)==basic) {
3798            dj[i]=0.0;
3799          } else if (model2->getStatus(i)==atLowerBound) {
3800            if (optimizationDirection_*dj[i]<tolerance) {
3801              if (optimizationDirection_*dj[i]<0.0) {
3802                //if (dj[i]<-1.0e-3)
3803                //printf("bad dj at lb %d %g\n",i,dj[i]);
3804                cost[i] -= dj[i];
3805                dj[i]=0.0;
3806              }
3807            } else {
3808              upper[i]=lower[i];
3809            }
3810          } else if (model2->getStatus(i)==atUpperBound) {
3811            if (optimizationDirection_*dj[i]>tolerance) {
3812              if (optimizationDirection_*dj[i]>0.0) {
3813                //if (dj[i]>1.0e-3)
3814                //printf("bad dj at ub %d %g\n",i,dj[i]);
3815                cost[i] -= dj[i];
3816                dj[i]=0.0;
3817              }
3818            } else {
3819              lower[i]=upper[i];
3820            }
3821          }
3822        }
3823        // just dual values pass
3824        //model2->setLogLevel(63);
3825        //model2->setFactorizationFrequency(1);
3826        model2->dual(2);
3827        memcpy(cost,saveCost,numberColumns*sizeof(double));
3828        delete [] saveCost;
3829        memcpy(lower,saveLower,numberColumns*sizeof(double));
3830        delete [] saveLower;
3831        memcpy(upper,saveUpper,numberColumns*sizeof(double));
3832        delete [] saveUpper;
3833      }
3834      // and finish
3835      // move solutions
3836      CoinMemcpyN(rowPrimal,
3837                  numberRows,model2->primalRowSolution());
3838      CoinMemcpyN(columnPrimal,
3839                  numberColumns,model2->primalColumnSolution());
3840    }
3841//     double saveScale = model2->objectiveScale();
3842//     model2->setObjectiveScale(1.0e-3);
3843//     model2->primal(2);
3844//    model2->setObjectiveScale(saveScale);
3845    model2->primal(1);
3846  } else if (barrierStatus==4&&crossover) {
3847    // memory problems
3848    model2->setPerturbation(savePerturbation);
3849    model2->createStatus();
3850    model2->dual();
3851  }
3852  model2->setMaximumIterations(saveMaxIts);
3853  delete [] rowPrimal;
3854  delete [] columnPrimal;
3855  delete [] rowDual;
3856  delete [] columnDual;
3857  if (saveLower) {
3858    pinfo2.postsolve(true);
3859    delete model2;
3860    model2=saveModel2;
3861    int numberRows = model2->numberRows();
3862    int numberColumns = model2->numberColumns();
3863    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
3864    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
3865    delete [] saveLower;
3866    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
3867    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
3868    delete [] saveUpper;
3869    model2->primal(1);
3870  }
3871  model2->setPerturbation(savePerturbation);
3872  return model2->status();
3873}
3874/* For strong branching.  On input lower and upper are new bounds
3875   while on output they are objective function values (>1.0e50 infeasible).
3876   Return code is 0 if nothing interesting, -1 if infeasible both
3877   ways and +1 if infeasible one way (check values to see which one(s))
3878*/
3879int ClpSimplex::strongBranching(int numberVariables,const int * variables,
3880                                double * newLower, double * newUpper,
3881                                double ** outputSolution,
3882                                int * outputStatus, int * outputIterations,
3883                                bool stopOnFirstInfeasible,
3884                                bool alwaysFinish)
3885{
3886  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
3887                                                    newLower,  newUpper,outputSolution,
3888                                                    outputStatus, outputIterations,
3889                                                    stopOnFirstInfeasible,
3890                                                    alwaysFinish);
3891}
3892/* Borrow model.  This is so we dont have to copy large amounts
3893   of data around.  It assumes a derived class wants to overwrite
3894   an empty model with a real one - while it does an algorithm.
3895   This is same as ClpModel one, but sets scaling on etc. */
3896void 
3897ClpSimplex::borrowModel(ClpModel & otherModel) 
3898{
3899  ClpModel::borrowModel(otherModel);
3900  createStatus();
3901  //ClpDualRowSteepest steep1;
3902  //setDualRowPivotAlgorithm(steep1);
3903  //ClpPrimalColumnSteepest steep2;
3904  //setPrimalColumnPivotAlgorithm(steep2);
3905}
3906void 
3907ClpSimplex::borrowModel(ClpSimplex & otherModel) 
3908{
3909  ClpModel::borrowModel(otherModel);
3910  createStatus();
3911  dualBound_ = otherModel.dualBound_;
3912  dualTolerance_ = otherModel.dualTolerance_;
3913  primalTolerance_ = otherModel.primalTolerance_;
3914  delete dualRowPivot_;
3915  dualRowPivot_ = otherModel.dualRowPivot_->clone(true);
3916  delete primalColumnPivot_;
3917  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(true);
3918  perturbation_ = otherModel.perturbation_;
3919  specialOptions_ = otherModel.specialOptions_;
3920  automaticScale_ = otherModel.automaticScale_;
3921}
3922/// Saves scalars for ClpSimplex
3923typedef struct {
3924  double optimizationDirection;
3925  double dblParam[ClpLastDblParam];
3926  double objectiveValue;
3927  double dualBound;
3928  double dualTolerance;
3929  double primalTolerance;
3930  double sumDualInfeasibilities;
3931  double sumPrimalInfeasibilities;
3932  double infeasibilityCost;
3933  int numberRows;
3934  int numberColumns;
3935  int intParam[ClpLastIntParam];
3936  int numberIterations;
3937  int problemStatus;
3938  int maximumIterations;
3939  int lengthNames;
3940  int numberDualInfeasibilities;
3941  int numberDualInfeasibilitiesWithoutFree;
3942  int numberPrimalInfeasibilities;
3943  int numberRefinements;
3944  int scalingFlag;
3945  int algorithm;
3946  unsigned int specialOptions;
3947  int dualPivotChoice;
3948  int primalPivotChoice;
3949  int matrixStorageChoice;
3950} Clp_scalars;
3951
3952int outDoubleArray(double * array, int length, FILE * fp)
3953{
3954  int numberWritten;
3955  if (array&&length) {
3956    numberWritten = fwrite(&length,sizeof(int),1,fp);
3957    if (numberWritten!=1)
3958      return 1;
3959    numberWritten = fwrite(array,sizeof(double),length,fp);
3960    if (numberWritten!=length)
3961      return 1;
3962  } else {
3963    length = 0;
3964    numberWritten = fwrite(&length,sizeof(int),1,fp);
3965    if (numberWritten!=1)
3966      return 1;
3967  }
3968  return 0;
3969}
3970// Save model to file, returns 0 if success
3971int
3972ClpSimplex::saveModel(const char * fileName)
3973{
3974  FILE * fp = fopen(fileName,"wb");
3975  if (fp) {
3976    Clp_scalars scalars;
3977    int i;
3978    CoinBigIndex numberWritten;
3979    // Fill in scalars
3980    scalars.optimizationDirection = optimizationDirection_;
3981    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
3982    scalars.objectiveValue = objectiveValue_;
3983    scalars.dualBound = dualBound_;
3984    scalars.dualTolerance = dualTolerance_;
3985    scalars.primalTolerance = primalTolerance_;
3986    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
3987    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
3988    scalars.infeasibilityCost = infeasibilityCost_;
3989    scalars.numberRows = numberRows_;
3990    scalars.numberColumns = numberColumns_;
3991    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
3992    scalars.numberIterations = numberIterations_;
3993    scalars.problemStatus = problemStatus_;
3994    scalars.maximumIterations = maximumIterations();
3995    scalars.lengthNames = lengthNames_;
3996    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
3997    scalars.numberDualInfeasibilitiesWithoutFree
3998      = numberDualInfeasibilitiesWithoutFree_;
3999    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
4000    scalars.numberRefinements = numberRefinements_;
4001    scalars.scalingFlag = scalingFlag_;
4002    scalars.algorithm = algorithm_;
4003    scalars.specialOptions = specialOptions_;
4004    scalars.dualPivotChoice = dualRowPivot_->type();
4005    scalars.primalPivotChoice = primalColumnPivot_->type();
4006    scalars.matrixStorageChoice = matrix_->type();
4007
4008    // put out scalars
4009    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
4010    if (numberWritten!=1)
4011      return 1;
4012    // strings
4013    CoinBigIndex length;
4014    for (i=0;i<ClpLastStrParam;i++) {
4015      length = strParam_[i].size();
4016      numberWritten = fwrite(&length,sizeof(int),1,fp);
4017      if (numberWritten!=1)
4018        return 1;
4019      if (length) {
4020        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
4021        if (numberWritten!=1)
4022          return 1;
4023      }
4024    }
4025    // arrays - in no particular order
4026    if (outDoubleArray(rowActivity_,numberRows_,fp))
4027        return 1;
4028    if (outDoubleArray(columnActivity_,numberColumns_,fp))
4029        return 1;
4030    if (outDoubleArray(dual_,numberRows_,fp))
4031        return 1;
4032    if (outDoubleArray(reducedCost_,numberColumns_,fp))
4033        return 1;
4034    if (outDoubleArray(rowLower_,numberRows_,fp))
4035        return 1;
4036    if (outDoubleArray(rowUpper_,numberRows_,fp))
4037        return 1;
4038    if (outDoubleArray(objective(),numberColumns_,fp))
4039        return 1;
4040    if (outDoubleArray(rowObjective_,numberRows_,fp))
4041        return 1;
4042    if (outDoubleArray(columnLower_,numberColumns_,fp))
4043        return 1;
4044    if (outDoubleArray(columnUpper_,numberColumns_,fp))
4045        return 1;
4046    if (ray_) {
4047      if (problemStatus_==1)
4048        if (outDoubleArray(ray_,numberRows_,fp))
4049          return 1;
4050      else if (problemStatus_==2)
4051        if (outDoubleArray(ray_,numberColumns_,fp))
4052          return 1;
4053      else
4054        if (outDoubleArray(NULL,0,fp))
4055          return 1;
4056    } else {
4057      if (outDoubleArray(NULL,0,fp))
4058        return 1;
4059    }
4060    if (status_&&(numberRows_+numberColumns_)>0) {
4061      length = numberRows_+numberColumns_;
4062      numberWritten = fwrite(&length,sizeof(int),1,fp);
4063      if (numberWritten!=1)
4064        return 1;
4065      numberWritten = fwrite(status_,sizeof(char),length, fp);
4066      if (numberWritten!=length)
4067        return 1;
4068    } else {
4069      length = 0;
4070      numberWritten = fwrite(&length,sizeof(int),1,fp);
4071      if (numberWritten!=1)
4072        return 1;
4073    }
4074    if (lengthNames_) {
4075      char * array = 
4076        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4077      char * put = array;
4078      assert (numberRows_ == (int) rowNames_.size());
4079      for (i=0;i<numberRows_;i++) {
4080        assert((int)rowNames_[i].size()<=lengthNames_);
4081        strcpy(put,rowNames_[i].c_str());
4082        put += lengthNames_+1;
4083      }
4084      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
4085      if (numberWritten!=numberRows_)
4086        return 1;
4087      put=array;
4088      assert (numberColumns_ == (int) columnNames_.size());
4089      for (i=0;i<numberColumns_;i++) {
4090        assert((int) columnNames_[i].size()<=lengthNames_);
4091        strcpy(put,columnNames_[i].c_str());
4092        put += lengthNames_+1;
4093      }
4094      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
4095      if (numberWritten!=numberColumns_)
4096        return 1;
4097      delete [] array;
4098    }
4099    // just standard type at present
4100    assert (matrix_->type()==1);
4101    assert (matrix_->getNumCols() == numberColumns_);
4102    assert (matrix_->getNumRows() == numberRows_);
4103    // we are going to save with gaps
4104    length = matrix_->getVectorStarts()[numberColumns_-1]
4105      + matrix_->getVectorLengths()[numberColumns_-1];
4106    numberWritten = fwrite(&length,sizeof(int),1,fp);
4107    if (numberWritten!=1)
4108      return 1;
4109    numberWritten = fwrite(matrix_->getElements(),
4110                           sizeof(double),length,fp);
4111    if (numberWritten!=length)
4112      return 1;
4113    numberWritten = fwrite(matrix_->getIndices(),
4114                           sizeof(int),length,fp);
4115    if (numberWritten!=length)
4116      return 1;
4117    numberWritten = fwrite(matrix_->getVectorStarts(),
4118                           sizeof(int),numberColumns_+1,fp);
4119    if (numberWritten!=numberColumns_+1)
4120      return 1;
4121    numberWritten = fwrite(matrix_->getVectorLengths(),
4122                           sizeof(int),numberColumns_,fp);
4123    if (numberWritten!=numberColumns_)
4124      return 1;
4125    // finished
4126    fclose(fp);
4127    return 0;
4128  } else {
4129    return -1;
4130  }
4131}
4132
4133int inDoubleArray(double * &array, int length, FILE * fp)
4134{
4135  int numberRead;
4136  int length2;
4137  numberRead = fread(&length2,sizeof(int),1,fp);
4138  if (numberRead!=1)
4139    return 1;
4140  if (length2) {
4141    // lengths must match
4142    if (length!=length2)
4143      return 2;
4144    array = new double[length];
4145    numberRead = fread(array,sizeof(double),length,fp);
4146    if (numberRead!=length)
4147      return 1;
4148  } 
4149  return 0;
4150}
4151/* Restore model from file, returns 0 if success,
4152   deletes current model */
4153int 
4154ClpSimplex::restoreModel(const char * fileName)
4155{
4156  FILE * fp = fopen(fileName,"rb");
4157  if (fp) {
4158    // Get rid of current model
4159    // save event handler in case already set
4160    ClpEventHandler * handler = eventHandler_->clone();
4161    ClpModel::gutsOfDelete();
4162    eventHandler_ = handler;
4163    gutsOfDelete(0);
4164    int i;
4165    for (i=0;i<6;i++) {
4166      rowArray_[i]=NULL;
4167      columnArray_[i]=NULL;
4168    }
4169    // get an empty factorization so we can set tolerances etc
4170    factorization_ = new ClpFactorization();
4171    // Say sparse
4172    factorization_->sparseThreshold(1);
4173    Clp_scalars scalars;
4174    CoinBigIndex numberRead;
4175
4176    // get scalars
4177    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
4178    if (numberRead!=1)
4179      return 1;
4180    // Fill in scalars
4181    optimizationDirection_ = scalars.optimizationDirection;
4182    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
4183    objectiveValue_ = scalars.objectiveValue;
4184    dualBound_ = scalars.dualBound;
4185    dualTolerance_ = scalars.dualTolerance;
4186    primalTolerance_ = scalars.primalTolerance;
4187    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
4188    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
4189    infeasibilityCost_ = scalars.infeasibilityCost;
4190    numberRows_ = scalars.numberRows;
4191    numberColumns_ = scalars.numberColumns;
4192    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
4193    numberIterations_ = scalars.numberIterations;
4194    problemStatus_ = scalars.problemStatus;
4195    setMaximumIterations(scalars.maximumIterations);
4196    lengthNames_ = scalars.lengthNames;
4197    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
4198    numberDualInfeasibilitiesWithoutFree_
4199      = scalars.numberDualInfeasibilitiesWithoutFree;
4200    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
4201    numberRefinements_ = scalars.numberRefinements;
4202    scalingFlag_ = scalars.scalingFlag;
4203    algorithm_ = scalars.algorithm;
4204    specialOptions_ = scalars.specialOptions;
4205    // strings
4206    CoinBigIndex length;
4207    for (i=0;i<ClpLastStrParam;i++) {
4208      numberRead = fread(&length,sizeof(int),1,fp);
4209      if (numberRead!=1)
4210        return 1;
4211      if (length) {
4212        char * array = new char[length+1];
4213        numberRead = fread(array,length,1,fp);
4214        if (numberRead!=1)
4215          return 1;
4216        array[length]='\0';
4217        strParam_[i]=array;
4218        delete [] array;
4219      }
4220    }
4221    // arrays - in no particular order
4222    if (inDoubleArray(rowActivity_,numberRows_,fp))
4223        return 1;
4224    if (inDoubleArray(columnActivity_,numberColumns_,fp))
4225        return 1;
4226    if (inDoubleArray(dual_,numberRows_,fp))
4227        return 1;
4228    if (inDoubleArray(reducedCost_,numberColumns_,fp))
4229        return 1;
4230    if (inDoubleArray(rowLower_,numberRows_,fp))
4231        return 1;
4232    if (inDoubleArray(rowUpper_,numberRows_,fp))
4233        return 1;
4234    double * objective;
4235    if (inDoubleArray(objective,numberColumns_,fp))
4236        return 1;
4237    delete objective_;
4238    objective_ = new ClpLinearObjective(objective,numberColumns_);
4239    delete [] objective;
4240    if (inDoubleArray(rowObjective_,numberRows_,fp))
4241        return 1;
4242    if (inDoubleArray(columnLower_,numberColumns_,fp))
4243        return 1;
4244    if (inDoubleArray(columnUpper_,numberColumns_,fp))
4245        return 1;
4246    if (problemStatus_==1) {
4247      if (inDoubleArray(ray_,numberRows_,fp))
4248        return 1;
4249    } else if (problemStatus_==2) {
4250      if (inDoubleArray(ray_,numberColumns_,fp))
4251        return 1;
4252    } else {
4253      // ray should be null
4254      numberRead = fread(&length,sizeof(int),1,fp);
4255      if (numberRead!=1)
4256        return 1;
4257      if (length)
4258        return 2;
4259    }
4260    delete [] status_;
4261    status_=NULL;
4262    // status region
4263    numberRead = fread(&length,sizeof(int),1,fp);
4264    if (numberRead!=1)
4265        return 1;
4266    if (length) {
4267      if (length!=numberRows_+numberColumns_)
4268        return 1;
4269      status_ = new char unsigned[length];
4270      numberRead = fread(status_,sizeof(char),length, fp);
4271      if (numberRead!=length)
4272        return 1;
4273    }
4274    if (lengthNames_) {
4275      char * array = 
4276        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4277      char * get = array;
4278      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
4279      if (numberRead!=numberRows_)
4280        return 1;
4281      rowNames_ = std::vector<std::string> ();
4282      rowNames_.resize(numberRows_);
4283      for (i=0;i<numberRows_;i++) {
4284        rowNames_.push_back(get);
4285        get += lengthNames_+1;
4286      }
4287      get = array;
4288      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
4289      if (numberRead!=numberColumns_)
4290        return 1;
4291      columnNames_ = std::vector<std::string> ();
4292      columnNames_.resize(numberColumns_);
4293      for (i=0;i<numberColumns_;i++) {
4294        columnNames_.push_back(get);
4295        get += lengthNames_+1;
4296      }
4297      delete [] array;
4298    }
4299    // Pivot choices
4300    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
4301    delete dualRowPivot_;
4302    switch ((scalars.dualPivotChoice&63)) {
4303    default:
4304      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
4305    case 1:
4306      // Dantzig
4307      dualRowPivot_ = new ClpDualRowDantzig();
4308      break;
4309    case 2:
4310      // Steepest - use mode
4311      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
4312      break;
4313    }
4314    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
4315    delete primalColumnPivot_;
4316    switch ((scalars.primalPivotChoice&63)) {
4317    default:
4318      printf("Need another primalPivot case %d\n",
4319             scalars.primalPivotChoice&63);
4320    case 1:
4321      // Dantzig
4322      primalColumnPivot_ = new ClpPrimalColumnDantzig();
4323      break;
4324    case 2:
4325      // Steepest - use mode
4326      primalColumnPivot_
4327        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
4328      break;
4329    }
4330    assert(scalars.matrixStorageChoice==1);
4331    delete matrix_;
4332    // get arrays
4333    numberRead = fread(&length,sizeof(int),1,fp);
4334    if (numberRead!=1)
4335      return 1;
4336    double * elements = new double[length];
4337    int * indices = new int[length];
4338    CoinBigIndex * starts = new CoinBigIndex[numberColumns_+1];
4339    int * lengths = new int[numberColumns_];
4340    numberRead = fread(elements, sizeof(double),length,fp);
4341    if (numberRead!=length)
4342      return 1;
4343    numberRead = fread(indices, sizeof(int),length,fp);
4344    if (numberRead!=length)
4345      return 1;
4346    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
4347    if (numberRead!=numberColumns_+1)
4348      return 1;
4349    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
4350    if (numberRead!=numberColumns_)
4351      return 1;
4352    // assign matrix
4353    CoinPackedMatrix * matrix = new CoinPackedMatrix();
4354    matrix->assignMatrix(true, numberRows_, numberColumns_,
4355                         length, elements, indices, starts, lengths);
4356    // and transfer to Clp
4357    matrix_ = new ClpPackedMatrix(matrix);
4358    // finished
4359    fclose(fp);
4360    return 0;
4361  } else {
4362    return -1;
4363  }
4364  return 0;
4365}
4366// value of incoming variable (in Dual)
4367double 
4368ClpSimplex::valueIncomingDual() const
4369{
4370  // Need value of incoming for list of infeasibilities as may be infeasible
4371  double valueIncoming = (dualOut_/alpha_)*directionOut_;
4372  if (directionIn_==-1)
4373    valueIncoming = upperIn_-valueIncoming;
4374  else
4375    valueIncoming = lowerIn_-valueIncoming;
4376  return valueIncoming;
4377}
4378// Sanity check on input data - returns true if okay
4379bool 
4380ClpSimplex::sanityCheck()
4381{
4382  // bad if empty
4383  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
4384    int infeasNumber[2];
4385    double infeasSum[2];
4386    problemStatus_=emptyProblem(infeasNumber,infeasSum,false);
4387    numberDualInfeasibilities_=infeasNumber[0];
4388    sumDualInfeasibilities_=infeasSum[0];
4389    numberPrimalInfeasibilities_=infeasNumber[1];
4390    sumPrimalInfeasibilities_=infeasSum[1];
4391    return false;
4392  }
4393  int numberBad ;
4394  double largestBound, smallestBound, minimumGap;
4395  double smallestObj, largestObj;
4396  int firstBad;
4397  int modifiedBounds=0;
4398  int i;
4399  numberBad=0;
4400  firstBad=-1;
4401  minimumGap=1.0e100;
4402  smallestBound=1.0e100;
4403  largestBound=0.0;
4404  smallestObj=1.0e100;
4405  largestObj=0.0;
4406  // If bounds are too close - fix
4407  double fixTolerance = 10.0*primalTolerance_;
4408  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
4409    double value;
4410    value = fabs(cost_[i]);
4411    if (value>1.0e50) {
4412      numberBad++;
4413      if (firstBad<0)
4414        firstBad=i;
4415    } else if (value) {
4416      if (value>largestObj)
4417        largestObj=value;
4418      if (value<smallestObj)
4419        smallestObj=value;
4420    }
4421    value=upper_[i]-lower_[i];
4422    if (value<-primalTolerance_) {
4423      numberBad++;
4424      if (firstBad<0)
4425        firstBad=i;
4426    } else if (value<=fixTolerance) {
4427      if (value) {
4428        // modify
4429        upper_[i] = lower_[i];
4430        modifiedBounds++;
4431      }
4432    } else {
4433      if (value<minimumGap)
4434        minimumGap=value;
4435    }
4436    if (lower_[i]>-1.0e100&&lower_[i]) {
4437      value = fabs(lower_[i]);
4438      if (value>largestBound)
4439        largestBound=value;
4440      if (value<smallestBound)
4441        smallestBound=value;
4442    }
4443    if (upper_[i]<1.0e100&&upper_[i]) {
4444      value = fabs(upper_[i]);
4445      if (value>largestBound)
4446        largestBound=value;
4447      if (value<smallestBound)
4448        smallestBound=value;
4449    }
4450  }
4451  if (largestBound)
4452    handler_->message(CLP_RIMSTATISTICS3,messages_)
4453      <<smallestBound
4454      <<largestBound
4455      <<minimumGap
4456      <<CoinMessageEol;
4457  minimumGap=1.0e100;
4458  smallestBound=1.0e100;
4459  largestBound=0.0;
4460  for (i=0;i<numberColumns_;i++) {
4461    double value;
4462    value = fabs(cost_[i]);
4463    if (value>1.0e50) {
4464      numberBad++;
4465      if (firstBad<0)
4466        firstBad=i;
4467    } else if (value) {
4468      if (value>largestObj)
4469        largestObj=value;
4470      if (value<smallestObj)
4471        smallestObj=value;
4472    }
4473    value=upper_[i]-lower_[i];
4474    if (value<-primalTolerance_) {
4475      numberBad++;
4476      if (firstBad<0)
4477        firstBad=i;
4478    } else if (value<=fixTolerance) {
4479      if (value) {
4480        // modify
4481        upper_[i] = lower_[i];
4482        modifiedBounds++;
4483      }
4484    } else {
4485      if (value<minimumGap)
4486        minimumGap=value;
4487    }
4488    if (lower_[i]>-1.0e100&&lower_[i]) {
4489      value = fabs(lower_[i]);
4490      if (value>largestBound)
4491        largestBound=value;
4492      if (value<smallestBound)
4493        smallestBound=value;
4494    }
4495    if (upper_[i]<1.0e100&&upper_[i]) {
4496      value = fabs(upper_[i]);
4497      if (value>largestBound)
4498        largestBound=value;
4499      if (value<smallestBound)
4500        smallestBound=value;
4501    }
4502  }
4503  char rowcol[]={'R','C'};
4504  if (numberBad) {
4505    handler_->message(CLP_BAD_BOUNDS,messages_)
4506      <<numberBad
4507      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
4508      <<CoinMessageEol;
4509    problemStatus_=4;
4510    return false;
4511  }
4512  if (modifiedBounds)
4513    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
4514      <<modifiedBounds
4515      <<CoinMessageEol;
4516  handler_->message(CLP_RIMSTATISTICS1,messages_)
4517    <<smallestObj
4518    <<largestObj
4519    <<CoinMessageEol;  if (largestBound)
4520    handler_->message(CLP_RIMSTATISTICS2,messages_)
4521      <<smallestBound
4522      <<largestBound
4523      <<minimumGap
4524      <<CoinMessageEol;
4525  return true;
4526}
4527// Set up status array (for OsiClp)
4528void 
4529ClpSimplex::createStatus() 
4530{
4531  if(!status_)
4532    status_ = new unsigned char [numberColumns_+numberRows_];
4533  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
4534  int i;
4535  // set column status to one nearest zero
4536  for (i=0;i<numberColumns_;i++) {
4537#if 0
4538    if (columnLower_[i]>=0.0) {
4539      setColumnStatus(i,atLowerBound);
4540    } else if (columnUpper_[i]<=0.0) {
4541      setColumnStatus(i,atUpperBound);
4542    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4543      // free
4544      setColumnStatus(i,isFree);
4545    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4546      setColumnStatus(i,atLowerBound);
4547    } else {
4548      setColumnStatus(i,atUpperBound);
4549    }
4550#else
4551    setColumnStatus(i,atLowerBound);
4552#endif
4553  }
4554  for (i=0;i<numberRows_;i++) {
4555    setRowStatus(i,basic);
4556  }
4557}
4558/* Sets up all slack basis and resets solution to
4559   as it was after initial load or readMps */
4560void ClpSimplex::allSlackBasis(bool resetSolution)
4561{
4562  createStatus();
4563  if (resetSolution) {
4564    // put back to as it was originally
4565    int i;
4566    // set column status to one nearest zero
4567    // But set value to zero if lb <0.0 and ub>0.0
4568    for (i=0;i<numberColumns_;i++) {
4569      if (columnLower_[i]>=0.0) {
4570        columnActivity_[i]=columnLower_[i];
4571        setColumnStatus(i,atLowerBound);
4572      } else if (columnUpper_[i]<=0.0) {
4573        columnActivity_[i]=columnUpper_[i];
4574        setColumnStatus(i,atUpperBound);
4575      } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4576        // free
4577        columnActivity_[i]=0.0;
4578        setColumnStatus(i,isFree);
4579      } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4580        columnActivity_[i]=0.0;
4581        setColumnStatus(i,atLowerBound);
4582      } else {
4583        columnActivity_[i]=0.0;
4584        setColumnStatus(i,atUpperBound);
4585      }
4586    }
4587  }
4588}
4589/* Loads a problem (the constraints on the
4590   rows are given by lower and upper bounds). If a pointer is 0 then the
4591   following values are the default:
4592   <ul>
4593   <li> <code>colub</code>: all columns have upper bound infinity
4594   <li> <code>collb</code>: all columns have lower bound 0
4595   <li> <code>rowub</code>: all rows have upper bound infinity
4596   <li> <code>rowlb</code>: all rows have lower bound -infinity
4597   <li> <code>obj</code>: all variables have 0 objective coefficient
4598   </ul>
4599*/
4600void 
4601ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
4602                    const double* collb, const double* colub,   
4603                    const double* obj,
4604                    const double* rowlb, const double* rowub,
4605                    const double * rowObjective)
4606{
4607  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4608                        rowObjective);
4609  createStatus();
4610}
4611void 
4612ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
4613                    const double* collb, const double* colub,   
4614                    const double* obj,
4615                    const double* rowlb, const double* rowub,
4616                    const double * rowObjective)
4617{
4618  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4619                        rowObjective);
4620  createStatus();
4621}
4622
4623/* Just like the other loadProblem() method except that the matrix is
4624   given in a standard column major ordered format (without gaps). */
4625void 
4626ClpSimplex::loadProblem (  const int numcols, const int numrows,
4627                    const CoinBigIndex* start, const int* index,
4628                    const double* value,
4629                    const double* collb, const double* colub,   
4630                    const double* obj,
4631                    const double* rowlb, const double* rowub,
4632                    const double * rowObjective)
4633{
4634  ClpModel::loadProblem(numcols, numrows, start, index, value,
4635                          collb, colub, obj, rowlb, rowub,
4636                          rowObjective);
4637  createStatus();
4638}
4639void 
4640ClpSimplex::loadProblem (  const int numcols, const int numrows,
4641                           const CoinBigIndex* start, const int* index,
4642                           const double* value,const int * length,
4643                           const double* collb, const double* colub,   
4644                           const double* obj,
4645                           const double* rowlb, const double* rowub,
4646                           const double * rowObjective)
4647{
4648  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
4649                          collb, colub, obj, rowlb, rowub,
4650                          rowObjective);
4651  createStatus();
4652}
4653// Read an mps file from the given filename
4654int 
4655ClpSimplex::readMps(const char *filename,
4656            bool keepNames,
4657            bool ignoreErrors)
4658{
4659  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
4660  createStatus();
4661  return status;
4662}
4663// Just check solution (for external use)
4664void 
4665ClpSimplex::checkSolution()
4666{
4667  // put in standard form
4668  createRim(7+8+16);
4669  dualTolerance_=dblParam_[ClpDualTolerance];
4670  primalTolerance_=dblParam_[ClpPrimalTolerance];
4671  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
4672  checkDualSolution();
4673  if (!numberDualInfeasibilities_&&
4674      !numberPrimalInfeasibilities_)
4675    problemStatus_=0;
4676  else
4677    problemStatus_=-1;
4678#ifdef CLP_DEBUG
4679  int i;
4680  double value=0.0;
4681  for (i=0;i<numberRows_+numberColumns_;i++)
4682    value += dj_[i]*solution_[i];
4683  printf("dual value %g, primal %g\n",value,objectiveValue());
4684#endif
4685  // release extra memory
4686  deleteRim(0);
4687}
4688/* Crash - at present just aimed at dual, returns
4689   -2 if dual preferred and crash basis created
4690   -1 if dual preferred and all slack basis preferred
4691   0 if basis going in was not all slack
4692   1 if primal preferred and all slack basis preferred
4693   2 if primal preferred and crash basis created.
4694   
4695   if gap between bounds <="gap" variables can be flipped
4696   ( If pivot -1 then can be made super basic!)
4697   
4698   If "pivot" is
4699   -1 No pivoting - always primal
4700   0 No pivoting (so will just be choice of algorithm)
4701   1 Simple pivoting e.g. gub
4702   2 Mini iterations
4703*/
4704int 
4705ClpSimplex::crash(double gap,int pivot)
4706{
4707  assert(!rowObjective_); // not coded
4708  int iColumn;
4709  int numberBad=0;
4710  int numberBasic=0;
4711  double dualTolerance=dblParam_[ClpDualTolerance];
4712  //double primalTolerance=dblParam_[ClpPrimalTolerance];
4713  int returnCode=0;
4714  // If no basis then make all slack one
4715  if (!status_)
4716    createStatus();
4717 
4718  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4719    if (getColumnStatus(iColumn)==basic)
4720      numberBasic++;
4721  }
4722  if (!numberBasic) {
4723    // all slack
4724    double * dj = new double [numberColumns_];
4725    double * solution = columnActivity_;
4726    const double * linearObjective = objective();
4727    //double objectiveValue=0.0;
4728    int iColumn;
4729    double direction = optimizationDirection_;
4730    // direction is actually scale out not scale in
4731    if (direction)
4732      direction = 1.0/direction;
4733    for (iColumn=0;iColumn<numberColumns_;iColumn++)
4734      dj[iColumn] = direction*linearObjective[iColumn];
4735    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4736      // assume natural place is closest to zero
4737      double lowerBound = columnLower_[iColumn];
4738      double upperBound = columnUpper_[iColumn];
4739      if (lowerBound>-1.0e20||upperBound<1.0e20) {
4740        bool atLower;
4741        if (fabs(upperBound)<fabs(lowerBound)) {
4742          atLower=false;
4743          setColumnStatus(iColumn,atUpperBound);
4744          solution[iColumn]=upperBound;
4745        } else {
4746          atLower=true;
4747          setColumnStatus(iColumn,atLowerBound);
4748          solution[iColumn]=lowerBound;
4749        }
4750        if (dj[iColumn]<-dualTolerance_) {
4751          // should be at upper bound
4752          if (atLower) {
4753            // can we flip
4754            if (upperBound-lowerBound<=gap) {
4755              columnActivity_[iColumn]=upperBound;
4756              setColumnStatus(iColumn,atUpperBound);
4757            } else if (pivot<0) {
4758              // set superbasic
4759              columnActivity_[iColumn]=lowerBound+gap;
4760              setColumnStatus(iColumn,superBasic);
4761            } else if (dj[iColumn]<-dualTolerance) {
4762              numberBad++;
4763            }
4764          }
4765        } else if (dj[iColumn]>dualTolerance_) {
4766          // should be at lower bound
4767          if (!atLower) {
4768            // can we flip
4769            if (upperBound-lowerBound<=gap) {
4770              columnActivity_[iColumn]=lowerBound;
4771              setColumnStatus(iColumn,atLowerBound);
4772            } else if (pivot<0) {
4773              // set superbasic
4774              columnActivity_[iColumn]=upperBound-gap;
4775              setColumnStatus(iColumn,superBasic);
4776            } else if (dj[iColumn]>dualTolerance) {
4777              numberBad++;
4778            }
4779          }
4780        }
4781      } else {
4782        // free
4783        setColumnStatus(iColumn,isFree);
4784        if (fabs(dj[iColumn])>dualTolerance) 
4785          numberBad++;
4786      }
4787    }
4788    if (numberBad||pivot) {
4789      if (pivot<=0) {
4790        delete [] dj;
4791        returnCode = 1;
4792      } else {
4793        // see if can be made dual feasible with gubs etc
4794        double * pi = new double[numberRows_];
4795        memset (pi,0,numberRows_*sizeof(double));
4796        int * way = new int[numberColumns_];
4797        int numberIn = 0;
4798
4799        // Get column copy
4800        CoinPackedMatrix * columnCopy = matrix();
4801        // Get a row copy in standard format
4802        CoinPackedMatrix copy;
4803        copy.reverseOrderedCopyOf(*columnCopy);
4804        // get matrix data pointers
4805        const int * column = copy.getIndices();
4806        const CoinBigIndex * rowStart = copy.getVectorStarts();
4807        const int * rowLength = copy.getVectorLengths(); 
4808        const double * elementByRow = copy.getElements();
4809        //const int * row = columnCopy->getIndices();
4810        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
4811        //const int * columnLength = columnCopy->getVectorLengths();
4812        //const double * element = columnCopy->getElements();
4813
4814
4815        // if equality row and bounds mean artificial in basis bad
4816        // then do anyway
4817
4818        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4819          // - if we want to reduce dj, + if we want to increase
4820          int thisWay = 100;
4821          double lowerBound = columnLower_[iColumn];
4822          double upperBound = columnUpper_[iColumn];
4823          if (upperBound>lowerBound) {
4824            switch(getColumnStatus(iColumn)) {
4825             
4826            case basic:
4827              thisWay=0;
4828            case ClpSimplex::isFixed:
4829              break;
4830            case isFree:
4831            case superBasic:
4832              if (dj[iColumn]<-dualTolerance) 
4833                thisWay = 1;
4834              else if (dj[iColumn]>dualTolerance) 
4835                thisWay = -1;
4836              else
4837                thisWay =0;
4838              break;
4839            case atUpperBound:
4840              if (dj[iColumn]>dualTolerance) 
4841                thisWay = -1;
4842              else if (dj[iColumn]<-dualTolerance) 
4843                thisWay = -3;
4844              else
4845                thisWay = -2;
4846              break;
4847            case atLowerBound:
4848              if (dj[iColumn]<-dualTolerance) 
4849                thisWay = 1;
4850              else if (dj[iColumn]>dualTolerance) 
4851                thisWay = 3;
4852              else
4853                thisWay = 2;
4854              break;
4855            }
4856          }
4857          way[iColumn] = thisWay;
4858        }
4859        /*if (!numberBad)
4860          printf("Was dual feasible before passes - rows %d\n",
4861          numberRows_);*/
4862        int lastNumberIn = -100000;
4863        int numberPasses=5;
4864        while (numberIn>lastNumberIn+numberRows_/100) {
4865          lastNumberIn = numberIn;
4866          // we need to maximize chance of doing good
4867          int iRow;
4868          for (iRow=0;iRow<numberRows_;iRow++) {
4869            double lowerBound = rowLower_[iRow];
4870            double upperBound = rowUpper_[iRow];
4871            if (getRowStatus(iRow)==basic) {
4872              // see if we can find a column to pivot on
4873              int j;
4874              // down is amount pi can go down
4875              double maximumDown = COIN_DBL_MAX;
4876              double maximumUp = COIN_DBL_MAX;
4877              double minimumDown =0.0;
4878              double minimumUp =0.0;
4879              int iUp=-1;
4880              int iDown=-1;
4881              int iUpB=-1;
4882              int iDownB=-1;
4883              if (lowerBound<-1.0e20)
4884                maximumUp = -1.0;
4885              if (upperBound>1.0e20)
4886                maximumDown = -1.0;
4887              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4888                int iColumn = column[j];
4889                double value = elementByRow[j];
4890                double djValue = dj[iColumn];
4891                /* way -
4892                   -3 - okay at upper bound with negative dj
4893                   -2 - marginal at upper bound with zero dj - can only decrease
4894                   -1 - bad at upper bound
4895                   0 - we can never pivot on this row
4896                   1 - bad at lower bound
4897                   2 - marginal at lower bound with zero dj - can only increase
4898                   3 - okay at lower bound with positive dj
4899                   100 - fine we can just ignore
4900                */
4901                if (way[iColumn]!=100) {
4902                  switch(way[iColumn]) {
4903                   
4904                  case -3:
4905                    if (value>0.0) {
4906                      if (maximumDown*value>-djValue) {
4907                        maximumDown = -djValue/value;
4908                        iDown = iColumn;
4909                      }
4910                    } else {
4911                      if (-maximumUp*value>-djValue) {
4912                        maximumUp = djValue/value;
4913                        iUp = iColumn;
4914                      }
4915                    }
4916                    break;
4917                  case -2:
4918                    if (value>0.0) {
4919                      maximumDown = 0.0;
4920                    } else {
4921                      maximumUp = 0.0;
4922                    }
4923                    break;
4924                  case -1:
4925                    // see if could be satisfied
4926                    // dj value > 0
4927                    if (value>0.0) {
4928                      maximumDown=0.0;
4929                      if (maximumUp*value<djValue-dualTolerance) {
4930                        maximumUp = 0.0; // would improve but not enough
4931                      } else {
4932                        if (minimumUp*value<djValue) {
4933                          minimumUp = djValue/value;
4934                          iUpB = iColumn;
4935                        }
4936                      }
4937                    } else {
4938                      maximumUp=0.0;
4939                      if (-maximumDown*value<djValue-dualTolerance) {
4940                        maximumDown = 0.0; // would improve but not enough
4941                      } else {
4942                        if (-minimumDown*value<djValue) {
4943                          minimumDown = -djValue/value;
4944                          iDownB = iColumn;
4945                        }
4946                      }
4947                    }
4948                   
4949                    break;
4950                  case 0:
4951                    maximumDown = -1.0;
4952                    maximumUp=-1.0;
4953                    break;
4954                  case 1:
4955                    // see if could be satisfied
4956                    // dj value < 0
4957                    if (value>0.0) {
4958                      maximumUp=0.0;
4959                      if (maximumDown*value<-djValue-dualTolerance) {
4960                        maximumDown = 0.0; // would improve but not enough
4961                      } else {
4962                        if (minimumDown*value<-djValue) {
4963                          minimumDown = -djValue/value;
4964                          iDownB = iColumn;
4965                        }
4966                      }
4967                    } else {
4968                      maximumDown=0.0;
4969                      if (-maximumUp*value<-djValue-dualTolerance) {
4970                        maximumUp = 0.0; // would improve but not enough
4971                      } else {
4972                        if (-minimumUp*value<-djValue) {
4973                          minimumUp = djValue/value;
4974                          iUpB = iColumn;
4975                        }
4976                      }
4977                    }
4978                   
4979                    break;
4980                  case 2:
4981                    if (value>0.0) {
4982                      maximumUp = 0.0;
4983                    } else {
4984                      maximumDown = 0.0;
4985                    }
4986                   
4987                    break;
4988                  case 3:
4989                    if (value>0.0) {
4990                      if (maximumUp*value>djValue) {
4991                        maximumUp = djValue/value;
4992                        iUp = iColumn;
4993                      }
4994                    } else {
4995                      if (-maximumDown*value>djValue) {
4996                        maximumDown = -djValue/value;
4997                        iDown = iColumn;
4998                      }
4999                    }
5000                   
5001                    break;
5002                  default:
5003                    break;
5004                  }
5005                }
5006              }
5007              if (iUpB>=0)
5008                iUp=iUpB;
5009              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
5010                iUp=-1;
5011              if (iDownB>=0)
5012                iDown=iDownB;
5013              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
5014                iDown=-1;
5015              if (iUp>=0||iDown>=0) {
5016                // do something
5017                if (iUp>=0&&iDown>=0) {
5018                  if (maximumDown>maximumUp)
5019                    iUp=-1;
5020                }
5021                double change;
5022                int kColumn;
5023                if (iUp>=0) {
5024                  kColumn=iUp;
5025                  change=maximumUp;
5026                  // just do minimum if was dual infeasible
5027                  // ? only if maximum large?
5028                  if (minimumUp>0.0)
5029                    change=minimumUp;
5030                  setRowStatus(iRow,atUpperBound);
5031                } else {
5032                  kColumn=iDown;
5033                  change=-maximumDown;
5034                  // just do minimum if was dual infeasible
5035                  // ? only if maximum large?
5036                  if (minimumDown>0.0)
5037                    change=-minimumDown;
5038                  setRowStatus(iRow,atLowerBound);
5039                }
5040                assert (fabs(change)<1.0e20);
5041                setColumnStatus(kColumn,basic);
5042                numberIn++;
5043                pi[iRow]=change;
5044                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
5045                  int iColumn = column[j];
5046                  double value = elementByRow[j];
5047                  double djValue = dj[iColumn]-change*value;
5048                  dj[iColumn]=djValue;
5049                  if (abs(way[iColumn])==1) {
5050                    numberBad--;
5051                    /*if (!numberBad)
5052                      printf("Became dual feasible at row %d out of %d\n",
5053                      iRow, numberRows_);*/
5054                    lastNumberIn=-1000000;
5055                  }
5056                  int thisWay = 100;
5057                  double lowerBound = columnLower_[iColumn];
5058                  double upperBound = columnUpper_[iColumn];
5059                  if (upperBound>lowerBound) {
5060                    switch(getColumnStatus(iColumn)) {
5061                     
5062                    case basic:
5063                      thisWay=0;
5064                    case isFixed:
5065                      break;
5066                    case isFree:
5067                    case superBasic:
5068                      if (djValue<-dualTolerance) 
5069                        thisWay = 1;
5070                      else if (djValue>dualTolerance) 
5071                        thisWay = -1;
5072                      else
5073                        { thisWay =0; abort();}
5074                      break;
5075                    case atUpperBound:
5076                      if (djValue>dualTolerance) 
5077                        { thisWay =-1; abort();}
5078                      else if (djValue<-dualTolerance) 
5079                        thisWay = -3;
5080                      else
5081                        thisWay = -2;
5082                      break;
5083                    case atLowerBound:
5084                      if (djValue<-dualTolerance) 
5085                        { thisWay =1; abort();}
5086                      else if (djValue>dualTolerance) 
5087                        thisWay = 3;
5088                      else
5089                        thisWay = 2;
5090                      break;
5091                    }
5092                  }
5093                  way[iColumn] = thisWay;
5094                }
5095              }
5096            }
5097          }
5098          if (numberIn==lastNumberIn||numberBad||pivot<2)
5099            break;
5100          if (!(--numberPasses))
5101            break;
5102          //printf("%d put in so far\n",numberIn);
5103        }
5104        // last attempt to flip
5105        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5106          double lowerBound = columnLower_[iColumn];
5107          double upperBound = columnUpper_[iColumn];
5108          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
5109            double djValue=dj[iColumn];
5110            switch(getColumnStatus(iColumn)) {
5111             
5112            case basic:
5113            case ClpSimplex::isFixed:
5114              break;
5115            case isFree:
5116            case superBasic:
5117              break;
5118            case atUpperBound:
5119              if (djValue>dualTolerance) {
5120                setColumnStatus(iColumn,atUpperBound);
5121                solution[iColumn]=upperBound;
5122              } 
5123              break;
5124            case atLowerBound:
5125              if (djValue<-dualTolerance) {
5126                setColumnStatus(iColumn,atUpperBound);
5127                solution[iColumn]=upperBound;
5128              }
5129              break;
5130            }
5131          }
5132        }
5133        delete [] pi;
5134        delete [] dj;
5135        delete [] way;
5136        handler_->message(CLP_CRASH,messages_)
5137          <<numberIn
5138          <<numberBad
5139          <<CoinMessageEol;
5140        returnCode =  -1;
5141      }
5142    } else {
5143      delete [] dj;
5144      returnCode =  -1;
5145    }
5146    //cleanStatus();
5147  }
5148  return returnCode;
5149}
5150/* Pivot in a variable and out a variable.  Returns 0 if okay,
5151   1 if inaccuracy forced re-factorization, -1 if would be singular.
5152   Also updates primal/dual infeasibilities.
5153   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
5154*/
5155int ClpSimplex::pivot()
5156{
5157  // scaling not allowed
5158  assert (!scalingFlag_);
5159  // assume In_ and Out_ are correct and directionOut_ set
5160  // (or In_ if flip
5161  lowerIn_ = lower_[sequenceIn_];
5162  valueIn_ = solution_[sequenceIn_];
5163  upperIn_ = upper_[sequenceIn_];
5164  dualIn_ = dj_[sequenceIn_];
5165  lowerOut_ = lower_[sequenceOut_];
5166  valueOut_ = solution_[sequenceOut_];
5167  upperOut_ = upper_[sequenceOut_];
5168  // for now assume primal is feasible (or in dual)
5169  dualOut_ = dj_[sequenceOut_];
5170  assert(fabs(dualOut_)<1.0e-6);
5171  bool roundAgain = true;
5172  int returnCode=0;
5173  while (roundAgain) {
5174    roundAgain=false;
5175    unpack(rowArray_[1]);
5176    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
5177    alpha_=0.0;
5178    int i;
5179    int * index = rowArray_[1]->getIndices();
5180    int number = rowArray_[1]->getNumElements();
5181    double * element = rowArray_[1]->denseVector();
5182    for (i=0;i<number;i++) {
5183      int ii = index[i];
5184      if ( pivotVariable_[ii]==sequenceOut_) {
5185        pivotRow_=ii;
5186        alpha_=element[pivotRow_];
5187        break;
5188      }
5189    }
5190    assert (fabs(alpha_)>1.0e-12);
5191    // we are going to subtract movement from current basic
5192    double movement;
5193    // see where incoming will go to
5194    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
5195      // flip so go to bound
5196      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
5197    } else {
5198      // get where outgoing needs to get to
5199      double outValue = (directionOut_<0) ? upperOut_ : lowerOut_;
5200      // solutionOut_ - movement*alpha_ == outValue
5201      movement = (valueOut_-outValue)/alpha_;
5202      // set directionIn_ correctly
5203      directionIn_ = (movement>0) ? 1 :-1;
5204    }
5205    // update primal solution
5206    for (i=0;i<number;i++) {
5207      int ii = index[i];
5208      // get column
5209      int ij = pivotVariable_[ii];
5210      solution_[ij] -= movement*element[ii];
5211      element[ii]=0.0;
5212    }
5213    rowArray_[1]->setNumElements(0);
5214    // see where something went to
5215    if (sequenceOut_<0) {
5216      if (directionIn_<0) {
5217        assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
5218        solution_[sequenceIn_]=upperIn_;
5219      } else {
5220        assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
5221        solution_[sequenceIn_]=lowerIn_;
5222      }
5223    } else {
5224      if (directionOut_<0) {
5225        assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
5226        solution_[sequenceOut_]=upperOut_;
5227      } else {
5228        assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
5229        solution_[sequenceOut_]=lowerOut_;
5230      }
5231      solution_[sequenceIn_]=valueIn_+movement;
5232    }
5233    double objectiveChange = dualIn_*movement;
5234    // update duals
5235    if (pivotRow_>=0) {
5236      assert (fabs(alpha_)>1.0e-8);
5237      double multiplier = dualIn_/alpha_;
5238      rowArray_[0]->insert(pivotRow_,multiplier);
5239      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
5240      // put row of tableau in rowArray[0] and columnArray[0]
5241      matrix_->transposeTimes(this,-1.0,
5242                              rowArray_[0],columnArray_[1],columnArray_[0]);
5243      // update column djs
5244      int i;
5245      int * index = columnArray_[0]->getIndices();
5246      int number = columnArray_[0]->getNumElements();
5247      double * element = columnArray_[0]->denseVector();
5248      for (i=0;i<number;i++) {
5249        int ii = index[i];
5250        dj_[ii] += element[ii];
5251        element[ii]=0.0;
5252      }
5253      columnArray_[0]->setNumElements(0);
5254      // and row djs
5255      index = rowArray_[0]->getIndices();
5256      number = rowArray_[0]->getNumElements();
5257      element = rowArray_[0]->denseVector();
5258      for (i=0;i<number;i++) {
5259        int ii = index[i];
5260        dj_[ii+numberColumns_] += element[ii];
5261        dual_[ii] = dj_[ii+numberColumns_];
5262        element[ii]=0.0;
5263      }
5264      rowArray_[0]->setNumElements(0);
5265      // check incoming
5266      assert (fabs(dj_[sequenceIn_])<1.0e-6);
5267    }
5268   
5269    // if stable replace in basis
5270    int updateStatus = factorization_->replaceColumn(this,
5271                                                     rowArray_[2],
5272                                                     rowArray_[1],
5273                                                     pivotRow_,
5274                                                     alpha_);
5275    bool takePivot=true;
5276    // if no pivots, bad update but reasonable alpha - take and invert
5277    if (updateStatus==2&&
5278        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
5279      updateStatus=4;
5280    if (updateStatus==1||updateStatus==4) {
5281      // slight error
5282      if (factorization_->pivots()>5||updateStatus==4) {
5283        returnCode=-1;
5284      }
5285    } else if (updateStatus==2) {
5286      // major error
5287      rowArray_[1]->clear();
5288      takePivot=false;
5289      if (factorization_->pivots()) {
5290        // refactorize here
5291        statusOfProblem();
5292        roundAgain=true;
5293      } else {
5294        returnCode=1;
5295      }
5296    } else if (updateStatus==3) {
5297      // out of memory
5298      // increase space if not many iterations
5299      if (factorization_->pivots()<
5300          0.5*factorization_->maximumPivots()&&
5301          factorization_->pivots()<200)
5302        factorization_->areaFactor(
5303                                   factorization_->areaFactor() * 1.1);
5304      returnCode =-1; // factorize now
5305    }
5306    if (takePivot) {
5307      int save = algorithm_;
5308      // make simple so always primal
5309      algorithm_=1;
5310      housekeeping(objectiveChange);
5311      algorithm_=save;
5312    }
5313  }
5314  if (returnCode == -1) {
5315    // refactorize here
5316    statusOfProblem();
5317  } else {
5318    // just for now - recompute anyway
5319    gutsOfSolution(NULL,NULL);
5320  }
5321  return returnCode;
5322}
5323
5324/* Pivot in a variable and choose an outgoing one.  Assumes primal
5325   feasible - will not go through a bound.  Returns step length in theta
5326   Returns ray in ray_ (or NULL if no pivot)
5327   Return codes as before but -1 means no acceptable pivot
5328*/
5329int ClpSimplex::primalPivotResult()
5330{
5331  assert (sequenceIn_>=0);
5332  valueIn_=solution_[sequenceIn_];
5333  lowerIn_=lower_[sequenceIn_];
5334  upperIn_=upper_[sequenceIn_];
5335  dualIn_=dj_[sequenceIn_];
5336
5337  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
5338  if (returnCode<0&&returnCode>-4) {
5339    return 0;
5340  } else {
5341    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
5342           returnCode);
5343    return -1;
5344  }
5345}
5346 
5347/* Pivot out a variable and choose an incoing one.  Assumes dual
5348   feasible - will not go through a reduced cost. 
5349   Returns step length in theta
5350   Returns ray in ray_ (or NULL if no pivot)
5351   Return codes as before but -1 means no acceptable pivot
5352*/
5353int 
5354ClpSimplex::dualPivotResult()
5355{
5356  return ((ClpSimplexDual *) this)->pivotResult();
5357}
5358// Factorization frequency
5359int 
5360ClpSimplex::factorizationFrequency() const
5361{
5362  if (factorization_)
5363    return factorization_->maximumPivots();
5364  else 
5365    return -1;
5366}
5367void 
5368ClpSimplex::setFactorizationFrequency(int value)
5369{
5370  if (factorization_)
5371    factorization_->maximumPivots(value);
5372}
5373// Common bits of coding for dual and primal
5374int 
5375ClpSimplex::startup(int ifValuesPass)
5376{
5377  // sanity check
5378  // bad if empty (trap here to avoid using bad matrix_)
5379  if (!matrix_||(!matrix_->getNumElements()&&objective_->type()<2)) {
5380    int infeasNumber[2];
5381    double infeasSum[2];
5382    problemStatus_=emptyProblem(infeasNumber,infeasSum);
5383    numberDualInfeasibilities_=infeasNumber[0];
5384    sumDualInfeasibilities_=infeasSum[0];
5385    numberPrimalInfeasibilities_=infeasNumber[1];
5386    sumPrimalInfeasibilities_=infeasSum[1];
5387    return 2;
5388  }
5389  pivotRow_=-1;
5390  sequenceIn_=-1;
5391  sequenceOut_=-1;
5392  secondaryStatus_=0;
5393
5394  primalTolerance_=dblParam_[ClpPrimalTolerance];
5395  dualTolerance_=dblParam_[ClpDualTolerance];
5396  if (problemStatus_!=10)
5397    numberIterations_=0;
5398
5399  // put in standard form (and make row copy)
5400  // create modifiable copies of model rim and do optional scaling
5401  bool goodMatrix=createRim(7+8+16,true);
5402
5403  if (goodMatrix) {
5404    // Model looks okay
5405    // Do initial factorization
5406    // and set certain stuff
5407    // We can either set increasing rows so ...IsBasic gives pivot row
5408    // or we can just increment iBasic one by one
5409    // for now let ...iBasic give pivot row
5410    factorization_->increasingRows(2);
5411    // row activities have negative sign
5412    factorization_->slackValue(-1.0);
5413    factorization_->zeroTolerance(1.0e-13);
5414    // Switch off dense (unless special option set)
5415    int saveThreshold = factorization_->denseThreshold();
5416    factorization_->setDenseThreshold(0);
5417    // If values pass then perturb (otherwise may be optimal so leave a bit)
5418    if (ifValuesPass) {
5419      // do perturbation if asked for
5420     
5421      if (perturbation_<100) {
5422        if (algorithm_>0&&(objective_->type()<2||!objective_->activated())) {
5423          ((ClpSimplexPrimal *) this)->perturb(0);
5424        } else if (algorithm_<0) {
5425        ((ClpSimplexDual *) this)->perturb();
5426        }
5427      }
5428    }
5429    // for primal we will change bounds using infeasibilityCost_
5430    if (nonLinearCost_==NULL&&algorithm_>0) {
5431      // get a valid nonlinear cost function
5432      delete nonLinearCost_;
5433      nonLinearCost_= new ClpNonLinearCost(this);
5434    }
5435   
5436    // loop round to clean up solution if values pass
5437    int numberThrownOut = -1;
5438    int totalNumberThrownOut=0;
5439    while(numberThrownOut) {
5440      int status = internalFactorize(ifValuesPass ? 10 : 0);
5441      if (status<0)
5442        return 1; // some error
5443      else
5444        numberThrownOut = status;
5445     
5446      // for this we need clean basis so it is after factorize
5447      if (!numberThrownOut)
5448        numberThrownOut = gutsOfSolution(  NULL,NULL,
5449                                         ifValuesPass!=0);
5450      else
5451        matrix_->rhsOffset(this,true); // redo rhs offset
5452      totalNumberThrownOut+= numberThrownOut;
5453     
5454    }
5455   
5456    if (totalNumberThrownOut)
5457      handler_->message(CLP_SINGULARITIES,messages_)
5458        <<totalNumberThrownOut
5459        <<CoinMessageEol;
5460    // Switch back dense
5461    factorization_->setDenseThreshold(saveThreshold);
5462   
5463    problemStatus_ = -1;
5464   
5465    // number of times we have declared optimality
5466    numberTimesOptimal_=0;
5467
5468    return 0;
5469  } else {
5470    // bad matrix
5471    return 2;
5472  }
5473   
5474}
5475
5476
5477void 
5478ClpSimplex::finish()
5479{
5480  // Get rid of some arrays and empty factorization
5481  deleteRim();
5482  // Skip message if changing algorithms
5483  if (problemStatus_!=10) {
5484    if (problemStatus_==-1)
5485      problemStatus_=4;
5486    assert(problemStatus_>=0&&problemStatus_<6);
5487    handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
5488      <<objectiveValue()
5489      <<CoinMessageEol;
5490  }
5491  factorization_->relaxAccuracyCheck(1.0);
5492  // get rid of any network stuff - could do more
5493  factorization_->cleanUp();
5494}
5495// Save data
5496ClpDataSave
5497ClpSimplex::saveData() 
5498{
5499  ClpDataSave saved;
5500  saved.dualBound_ = dualBound_;
5501  saved.infeasibilityCost_ = infeasibilityCost_;
5502  saved.sparseThreshold_ = factorization_->sparseThreshold();
5503  saved.pivotTolerance_ = factorization_->pivotTolerance();
5504  saved.perturbation_ = perturbation_;
5505  // Progress indicator
5506  delete progress_;
5507  progress_ = new ClpSimplexProgress (this);
5508  return saved;
5509}
5510// Restore data
5511void 
5512ClpSimplex::restoreData(ClpDataSave saved)
5513{
5514  factorization_->sparseThreshold(saved.sparseThreshold_);
5515  factorization_->pivotTolerance(saved.pivotTolerance_);
5516  perturbation_ = saved.perturbation_;
5517  infeasibilityCost_ = saved.infeasibilityCost_;
5518  dualBound_ = saved.dualBound_;
5519  delete progress_;
5520  progress_=NULL;
5521}
5522// To flag a variable (not inline to allow for column generation)
5523void 
5524ClpSimplex::setFlagged( int sequence)
5525{
5526  status_[sequence] |= 64;
5527  matrix_->generalExpanded(this,7,sequence);
5528  lastFlaggedIteration_=numberIterations_;
5529}
5530/* Factorizes and returns true if optimal.  Used by user */
5531bool
5532ClpSimplex::statusOfProblem(bool initial)
5533{
5534  // is factorization okay?
5535  if (initial) {
5536    // First time - allow singularities
5537    int numberThrownOut = -1;
5538    int totalNumberThrownOut=0;
5539    while(numberThrownOut) {
5540      int status = internalFactorize(0);
5541      if (status<0)
5542        return false; // some error
5543      else
5544        numberThrownOut = status;
5545     
5546      // for this we need clean basis so it is after factorize
5547      //if (!numberThrownOut)
5548      //numberThrownOut = gutsOfSolution(  NULL,NULL,
5549      //                                   false);
5550      //else
5551      //matrix_->rhsOffset(this,true); // redo rhs offset
5552      totalNumberThrownOut+= numberThrownOut;
5553     
5554    }
5555   
5556    if (totalNumberThrownOut)
5557      handler_->message(CLP_SINGULARITIES,messages_)
5558        <<totalNumberThrownOut
5559        <<CoinMessageEol;
5560  } else {
5561#ifndef NDEBUG   
5562    int returnCode=internalFactorize(1);
5563    assert (!returnCode);
5564#else
5565    internalFactorize(1);
5566#endif
5567  }
5568  // put back original costs and then check
5569  // also move to work arrays
5570  createRim(4+32);
5571  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
5572  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
5573  gutsOfSolution(NULL,NULL);
5574  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
5575  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
5576  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
5577  deleteRim(-1);
5578  return (primalFeasible()&&dualFeasible());
5579}
5580/* Return model - updates any scalars */
5581void 
5582ClpSimplex::returnModel(ClpSimplex & otherModel)
5583{
5584  ClpModel::returnModel(otherModel);
5585  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
5586  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
5587  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
5588  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
5589  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
5590  otherModel.columnDualSequence_ = columnDualSequence_;
5591  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
5592  otherModel.rowDualSequence_ = rowDualSequence_;
5593  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5594  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
5595  otherModel.largestPrimalError_ = largestPrimalError_;
5596  otherModel.largestDualError_ = largestDualError_;
5597  otherModel.largestSolutionError_ = largestSolutionError_;
5598  otherModel.alpha_ = alpha_;
5599  otherModel.theta_ = theta_;
5600  otherModel.lowerIn_ = lowerIn_;
5601  otherModel.valueIn_ = valueIn_;
5602  otherModel.upperIn_ = upperIn_;
5603  otherModel.dualIn_ = dualIn_;
5604  otherModel.sequenceIn_ = sequenceIn_;
5605  otherModel.directionIn_ = directionIn_;
5606  otherModel.lowerOut_ = lowerOut_;
5607  otherModel.valueOut_ = valueOut_;
5608  otherModel.upperOut_ = upperOut_;
5609  otherModel.dualOut_ = dualOut_;
5610  otherModel.sequenceOut_ = sequenceOut_;
5611  otherModel.directionOut_ = directionOut_;
5612  otherModel.pivotRow_ = pivotRow_;
5613  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
5614  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
5615  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
5616    numberDualInfeasibilitiesWithoutFree_;
5617  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
5618  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
5619  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
5620  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
5621  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
5622}
5623/* Constructs a non linear cost from list of non-linearities (columns only)
5624   First lower of each column is taken as real lower
5625   Last lower is taken as real upper and cost ignored
5626   
5627   Returns nonzero if bad data e.g. lowers not monotonic
5628*/
5629int 
5630ClpSimplex::createPiecewiseLinearCosts(const int * starts,
5631                                       const double * lower, const double * gradient)
5632{
5633  delete nonLinearCost_;
5634  // Set up feasible bounds and check monotonicity
5635  int iColumn;
5636  int returnCode=0;
5637
5638  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5639    int iIndex = starts[iColumn];
5640    int end = starts[iColumn+1]-1;
5641    columnLower_[iColumn] = lower[iIndex];
5642    columnUpper_[iColumn] = lower[end];
5643    double value = columnLower_[iColumn];
5644    iIndex++;
5645    for (;iIndex<end;iIndex++) {
5646      if (lower[iIndex]<value)
5647        returnCode++; // not monotonic
5648      value=lower[iIndex];
5649    }
5650  }
5651  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
5652  specialOptions_ |= 2; // say keep
5653  return returnCode;
5654}
5655/* For advanced use.  When doing iterative solves things can get
5656   nasty so on values pass if incoming solution has largest
5657   infeasibility < incomingInfeasibility throw out variables
5658   from basis until largest infeasibility < allowedInfeasibility
5659   or incoming largest infeasibility.
5660   If allowedInfeasibility>= incomingInfeasibility this is
5661   always possible altough you may end up with an all slack basis.
5662   
5663   Defaults are 1.0,10.0
5664*/
5665void 
5666ClpSimplex::setValuesPassAction(float incomingInfeasibility,
5667                                float allowedInfeasibility)
5668{
5669  incomingInfeasibility_=incomingInfeasibility;
5670  allowedInfeasibility_=allowedInfeasibility;
5671  assert(incomingInfeasibility_>=0.0);
5672  assert(allowedInfeasibility_>=incomingInfeasibility_);
5673}