source: branches/devel-1/ClpSimplex.cpp @ 2351

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

Some tidying up

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 89.2 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include <math.h>
12
13#include "CoinHelperFunctions.hpp"
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16#include "ClpPackedMatrix.hpp"
17#include "CoinIndexedVector.hpp"
18#include "ClpDualRowDantzig.hpp"
19#include "ClpDualRowSteepest.hpp"
20#include "ClpPrimalColumnDantzig.hpp"
21#include "ClpPrimalColumnSteepest.hpp"
22#include "ClpNonLinearCost.hpp"
23#include "ClpMessage.hpp"
24#include <cfloat>
25
26#include <string>
27#include <stdio.h>
28#include <iostream>
29//#############################################################################
30
31ClpSimplex::ClpSimplex () :
32
33  ClpModel(),
34  columnPrimalInfeasibility_(0.0),
35  columnPrimalSequence_(-2),
36  rowPrimalInfeasibility_(0.0),
37  rowPrimalSequence_(-2), 
38  columnDualInfeasibility_(0.0),
39  columnDualSequence_(-2),
40  rowDualInfeasibility_(0.0),
41  rowDualSequence_(-2),
42  primalToleranceToGetOptimal_(-1.0),
43  remainingDualInfeasibility_(0.0),
44  largeValue_(1.0e15),
45  largestPrimalError_(0.0),
46  largestDualError_(0.0),
47  largestSolutionError_(0.0),
48  dualBound_(1.0e6),
49  lower_(NULL),
50  rowLowerWork_(NULL),
51  columnLowerWork_(NULL),
52  upper_(NULL),
53  rowUpperWork_(NULL),
54  columnUpperWork_(NULL),
55  cost_(NULL),
56  rowObjectiveWork_(NULL),
57  objectiveWork_(NULL),
58  alpha_(0.0),
59  theta_(0.0),
60  lowerIn_(0.0),
61  valueIn_(0.0),
62  upperIn_(0.0),
63  dualIn_(0.0),
64  sequenceIn_(-1),
65  directionIn_(-1),
66  lowerOut_(-1),
67  valueOut_(-1),
68  upperOut_(-1),
69  dualOut_(-1),
70  sequenceOut_(-1),
71  directionOut_(-1),
72  pivotRow_(-1),
73  dj_(NULL),
74  rowReducedCost_(NULL),
75  reducedCostWork_(NULL),
76  solution_(NULL),
77  rowActivityWork_(NULL),
78  columnActivityWork_(NULL),
79  dualTolerance_(0.0),
80  primalTolerance_(0.0),
81  sumDualInfeasibilities_(0.0),
82  numberDualInfeasibilities_(0),
83  numberDualInfeasibilitiesWithoutFree_(0),
84  sumPrimalInfeasibilities_(0.0),
85  numberPrimalInfeasibilities_(0),
86  pivotVariable_(NULL),
87  factorization_(NULL),
88  numberRefinements_(0),
89  rowScale_(NULL),
90  savedSolution_(NULL),
91  columnScale_(NULL),
92  scalingFlag_(1),
93  numberTimesOptimal_(0),
94  changeMade_(1),
95  algorithm_(0),
96  forceFactorization_(-1),
97  perturbation_(100),
98  infeasibilityCost_(1.0e7),
99  nonLinearCost_(NULL),
100  specialOptions_(0),
101  lastBadIteration_(-999999),
102  numberFake_(0),
103  progressFlag_(0),
104  sumOfRelaxedDualInfeasibilities_(0.0),
105  sumOfRelaxedPrimalInfeasibilities_(0.0)
106
107{
108  int i;
109  for (i=0;i<6;i++) {
110    rowArray_[i]=NULL;
111    columnArray_[i]=NULL;
112  }
113  saveStatus_=NULL;
114  // get an empty factorization so we can set tolerances etc
115  factorization_ = new ClpFactorization();
116  // Say sparse
117  factorization_->sparseThreshold(1);
118  // say Steepest pricing
119  dualRowPivot_ = new ClpDualRowSteepest();
120  // say Steepest pricing
121  primalColumnPivot_ = new ClpPrimalColumnSteepest();
122 
123}
124
125
126//-----------------------------------------------------------------------------
127
128ClpSimplex::~ClpSimplex ()
129{
130  gutsOfDelete(0);
131}
132//#############################################################################
133void ClpSimplex::setLargeValue( double value) 
134{
135  if (value>0.0&&value<DBL_MAX)
136    largeValue_=value;
137}
138int
139ClpSimplex::gutsOfSolution ( const double * rowActivities,
140                             const double * columnActivities,
141                             bool valuesPass)
142{
143
144
145  // if values pass, save values of basic variables
146  double * save = NULL;
147  double oldValue=0.0;
148  if (valuesPass) {
149    assert(algorithm_>0); // only primal at present
150    assert(nonLinearCost_);
151    int iRow;
152    checkPrimalSolution( rowActivities, columnActivities);
153    // get correct bounds on all variables
154    nonLinearCost_->checkInfeasibilities();
155    oldValue = nonLinearCost_->largestInfeasibility();
156    save = new double[numberRows_];
157    for (iRow=0;iRow<numberRows_;iRow++) {
158      int iPivot=pivotVariable_[iRow];
159      save[iRow] = solution_[iPivot];
160    }
161  }
162  // do work
163  computePrimals(rowActivities, columnActivities);
164  double objectiveModification = 0.0;
165  if (algorithm_>0&&nonLinearCost_!=NULL) {
166    // primal algorithm
167    // get correct bounds on all variables
168    nonLinearCost_->checkInfeasibilities();
169    objectiveModification += nonLinearCost_->changeInCost();
170    if (nonLinearCost_->numberInfeasibilities())
171      handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
172        <<nonLinearCost_->changeInCost()
173        <<nonLinearCost_->numberInfeasibilities()
174        <<CoinMessageEol;
175  }
176  if (valuesPass) {
177#ifdef CLP_DEBUG
178    std::cout<<"Largest given infeasibility "<<oldValue
179             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
180#endif
181    int numberOut=0;
182    if (oldValue<1.0
183        &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
184        largestPrimalError_>1.0e-3) {
185      // throw out up to 1000 structurals
186      int iRow;
187      int * sort = new int[numberRows_];
188      // first put back solution and store difference
189      for (iRow=0;iRow<numberRows_;iRow++) {
190        int iPivot=pivotVariable_[iRow];
191        double difference = fabs(solution_[iPivot]=save[iRow]);
192        solution_[iPivot]=save[iRow];
193        save[iRow]=difference;
194      }
195      for (iRow=0;iRow<numberRows_;iRow++) {
196        int iPivot=pivotVariable_[iRow];
197
198        if (iPivot<numberColumns_) {
199          // column
200          double difference= save[iRow];
201          if (difference>1.0e-4) {
202            sort[numberOut]=iPivot;
203            save[numberOut++]=difference;
204          }
205        }
206      }
207      CoinSort_2(save, save + numberOut, sort,
208                 CoinFirstGreater_2<double, int>());
209      numberOut = min(1000,numberOut);
210      for (iRow=0;iRow<numberOut;iRow++) {
211        int iColumn=sort[iRow];
212        setColumnStatus(iColumn,superBasic);
213
214      }
215      delete [] sort;
216    }
217    delete [] save;
218    if (numberOut)
219      return numberOut;
220  }
221
222  computeDuals();
223
224  // now check solutions
225  checkPrimalSolution( rowActivities, columnActivities);
226  objectiveValue_ += objectiveModification;
227  checkDualSolution();
228  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
229                               largestDualError_>1.0e-2)) 
230    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
231      <<largestPrimalError_
232      <<largestDualError_
233      <<CoinMessageEol;
234  return 0;
235}
236void
237ClpSimplex::computePrimals ( const double * rowActivities,
238                                     const double * columnActivities)
239{
240
241  //work space
242  CoinIndexedVector  * workSpace = rowArray_[0];
243
244  double * array = new double [numberRows_];
245  double * save = new double [numberRows_];
246  double * previous = new double [numberRows_];
247
248  // accumulate non basic stuff
249  ClpFillN(array,numberRows_,0.0);
250
251  int iRow;
252  // order is this way for scaling
253  // Use whole matrix every time to make it easier for ClpMatrixBase
254  // So zero out basic
255  if (columnActivities!=columnActivityWork_)
256    ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
257  if (rowActivities!=rowActivityWork_)
258    ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
259  for (iRow=0;iRow<numberRows_;iRow++) {
260    int iPivot=pivotVariable_[iRow];
261    solution_[iPivot] = 0.0;
262  }
263  times(-1.0,columnActivityWork_,array);
264
265  for (iRow=0;iRow<numberRows_;iRow++) {
266    array[iRow] += rowActivityWork_[iRow];
267  }
268
269  // Ftran adjusted RHS and iterate to improve accuracy
270  double lastError=DBL_MAX;
271  int iRefine;
272  double * work = workSpace->denseVector();
273  factorization_->updateColumn(workSpace,array);
274
275  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
276    // save in case we want to restore
277    ClpDisjointCopyN ( array, numberRows_ , save);
278
279    // put solution in correct place
280    for (iRow=0;iRow<numberRows_;iRow++) {
281      int iPivot=pivotVariable_[iRow];
282      solution_[iPivot] = array[iRow];
283    }
284    // check Ax == b  (for all)
285    times(-1.0,columnActivityWork_,work);
286    for (iRow=0;iRow<numberRows_;iRow++) {
287      work[iRow] += rowActivityWork_[iRow];
288    }
289
290    largestPrimalError_=0.0;
291    for (iRow=0;iRow<numberRows_;iRow++) {
292      if (fabs(work[iRow])>largestPrimalError_) {
293        largestPrimalError_=fabs(work[iRow]);
294      }
295      //work[iRow] -= save[iRow];
296    }
297    if (largestPrimalError_>=lastError) {
298      // restore
299      double * temp = array;
300      array = previous;
301      previous=temp;
302      break;
303    }
304    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
305      // try and make better
306      // save this
307      double * temp = array;
308      array = previous;
309      previous=temp;
310      double multiplier = 131072.0;
311      for (iRow=0;iRow<numberRows_;iRow++) {
312        array[iRow] = multiplier*work[iRow];
313        work[iRow]=0.0;
314      }
315      lastError=largestPrimalError_;
316      factorization_->updateColumn(workSpace,array);
317      multiplier = 1.0/multiplier;
318      for (iRow=0;iRow<numberRows_;iRow++) {
319        array[iRow] = previous[iRow] + multiplier*array[iRow];
320        work[iRow]=0.0;
321      }
322    }
323  }
324
325  // solution as accurate as we are going to get
326  ClpFillN(work,numberRows_,0.0);
327  // put solution in correct place
328  for (iRow=0;iRow<numberRows_;iRow++) {
329    int iPivot=pivotVariable_[iRow];
330    solution_[iPivot] = array[iRow];
331  }
332  delete [] previous;
333  delete [] array;
334  delete [] save;
335}
336// now dual side
337void
338ClpSimplex::computeDuals()
339{
340  double slackValue = factorization_->slackValue();
341  //work space
342  CoinIndexedVector  * workSpace = rowArray_[0];
343
344  double * array = new double [numberRows_];
345  double * save = new double [numberRows_];
346  double * previous = new double [numberRows_];
347
348  int iRow;
349#ifdef CLP_DEBUG
350  workSpace->checkClear();
351#endif
352  for (iRow=0;iRow<numberRows_;iRow++) {
353    int iPivot=pivotVariable_[iRow];
354    array[iRow]=cost_[iPivot];
355  }
356  ClpDisjointCopyN ( array, numberRows_ , save);
357
358  // Btran basic costs and get as accurate as possible
359  double lastError=DBL_MAX;
360  int iRefine;
361  double * work = workSpace->denseVector();
362  factorization_->updateColumnTranspose(workSpace,array);
363  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
364    // check basic reduced costs zero
365    largestDualError_=0.0;
366    // would be faster to do just for basic but this reduces code
367    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
368    transposeTimes(-1.0,array,reducedCostWork_);
369    for (iRow=0;iRow<numberRows_;iRow++) {
370      int iPivot=pivotVariable_[iRow];
371      if (iPivot>=numberColumns_) {
372        // slack
373        //work[iRow] += slackValue*array[iPivot-numberColumns_];
374        //work[iRow] += rowObjectiveWork_[iPivot-numberColumns_];
375        work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
376        - slackValue*array[iPivot-numberColumns_];
377      } else {
378        // column
379        work[iRow] = reducedCostWork_[iPivot];
380      }
381      if (fabs(work[iRow])>largestDualError_) {
382        largestDualError_=fabs(work[iRow]);
383      }
384    }
385    if (largestDualError_>=lastError) {
386      // restore
387      double * temp = array;
388      array = previous;
389      previous=temp;
390      break;
391    }
392    if (iRefine<numberRefinements_&&largestDualError_>1.0e-20) {
393      // try and make better
394      // save this
395      double * temp = array;
396      array = previous;
397      previous=temp;
398      double multiplier = 131072.0;
399      for (iRow=0;iRow<numberRows_;iRow++) {
400        array[iRow] = multiplier*work[iRow];
401        work[iRow]=0.0;
402      }
403      lastError=largestDualError_;
404      factorization_->updateColumnTranspose(workSpace,array);
405      multiplier = 1.0/multiplier;
406      for (iRow=0;iRow<numberRows_;iRow++) {
407        array[iRow] = previous[iRow] + multiplier*array[iRow];
408        work[iRow]=0.0;
409      }
410    }
411  }
412  ClpFillN(work,numberRows_,0.0);
413  // now look at dual solution
414  for (iRow=0;iRow<numberRows_;iRow++) {
415    // slack
416    double value = array[iRow];
417    dual_[iRow]=value;
418    value = rowObjectiveWork_[iRow] - value*slackValue;
419    rowReducedCost_[iRow]=value;
420  }
421  ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
422  transposeTimes(-1.0,dual_,reducedCostWork_);
423  delete [] previous;
424  delete [] array;
425  delete [] save;
426}
427/* Given an existing factorization computes and checks
428   primal and dual solutions.  Uses input arrays for variables at
429   bounds.  Returns feasibility states */
430int ClpSimplex::getSolution ( const double * rowActivities,
431                               const double * columnActivities)
432{
433  if (!factorization_->status()) {
434    // put in standard form
435    createRim(7+8+16);
436    // do work
437    gutsOfSolution ( rowActivities, columnActivities);
438    // release extra memory
439    deleteRim(false);
440  }
441  return factorization_->status();
442}
443/* Given an existing factorization computes and checks
444   primal and dual solutions.  Uses current problem arrays for
445   bounds.  Returns feasibility states */
446int ClpSimplex::getSolution ( )
447{
448  double * rowActivities = new double[numberRows_];
449  double * columnActivities = new double[numberColumns_];
450  ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
451  ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
452  int status = getSolution( rowActivities, columnActivities);
453  delete [] rowActivities;
454  delete [] columnActivities;
455  return status;
456}
457// Factorizes using current basis.  This is for external use
458// Return codes are as from ClpFactorization
459int ClpSimplex::factorize ()
460{
461  // put in standard form
462  createRim(7+8+16);
463  // do work
464  int status = internalFactorize(-1);
465  // release extra memory
466  deleteRim(false);
467
468  return status;
469}
470
471/* Factorizes using current basis. 
472      solveType - 1 iterating, 0 initial, -1 external
473      If 10 added then in primal values pass
474*/
475/* Return codes are as from ClpFactorization unless initial factorization
476   when total number of singularities is returned */
477int ClpSimplex::internalFactorize ( int solveType)
478{
479
480  int iRow,iColumn;
481  int totalSlacks=numberRows_;
482
483  bool valuesPass=false;
484  if (solveType>=10) {
485    valuesPass=true;
486    solveType -= 10;
487  }
488  if (!solveType) {
489    if (!valuesPass) {
490      // not values pass so set to bounds
491      bool allSlack=true;
492      if (status_) {
493        for (iRow=0;iRow<numberRows_;iRow++) {
494          if (getRowStatus(iRow)!=basic) {
495            allSlack=false;
496            break;
497          }
498        }
499      }
500      if (!allSlack) {
501        // set values from warm start (if sensible)
502        int numberBasic=0;
503        for (iRow=0;iRow<numberRows_;iRow++) {
504          switch(getRowStatus(iRow)) {
505           
506          case basic:
507            numberBasic++;
508            break;
509          case isFree:
510            assert(rowLowerWork_[iRow]<-largeValue_);
511            assert(rowUpperWork_[iRow]>largeValue_);
512            rowActivityWork_[iRow]=0.0;
513            break;
514          case atUpperBound:
515            rowActivityWork_[iRow]=rowUpperWork_[iRow];
516            if (rowActivityWork_[iRow]>largeValue_) {
517              if (rowLowerWork_[iRow]>-largeValue_) {
518                rowActivityWork_[iRow]=rowLowerWork_[iRow];
519                setRowStatus(iRow,atLowerBound);
520              } else {
521                // say free
522                setRowStatus(iRow,isFree);
523                rowActivityWork_[iRow]=0.0;
524              }
525            }
526            break;
527          case atLowerBound:
528            rowActivityWork_[iRow]=rowLowerWork_[iRow];
529            if (rowActivityWork_[iRow]<-largeValue_) {
530              if (rowUpperWork_[iRow]<largeValue_) {
531                rowActivityWork_[iRow]=rowUpperWork_[iRow];
532                setRowStatus(iRow,atUpperBound);
533              } else {
534                // say free
535                setRowStatus(iRow,isFree);
536                rowActivityWork_[iRow]=0.0;
537              }
538            }
539            break;
540          case superBasic:
541            if (rowUpperWork_[iRow]>largeValue_) {
542              if (rowLowerWork_[iRow]>-largeValue_) {
543                rowActivityWork_[iRow]=rowLowerWork_[iRow];
544                setRowStatus(iRow,atLowerBound);
545              } else {
546                // say free
547                setRowStatus(iRow,isFree);
548                rowActivityWork_[iRow]=0.0;
549              }
550            } else {
551              if (rowLowerWork_[iRow]>-largeValue_) {
552                // set to nearest
553                if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
554                    <fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])) {
555                  rowActivityWork_[iRow]=rowLowerWork_[iRow];
556                  setRowStatus(iRow,atLowerBound);
557                } else {
558                  rowActivityWork_[iRow]=rowUpperWork_[iRow];
559                  setRowStatus(iRow,atUpperBound);
560                }
561              } else {
562                rowActivityWork_[iRow]=rowUpperWork_[iRow];
563                setRowStatus(iRow,atUpperBound);
564              }
565            }
566            break;
567          }
568        }
569        totalSlacks=numberBasic;
570        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
571          switch(getColumnStatus(iColumn)) {
572           
573          case basic:
574            if (numberBasic==numberRows_) {
575              // take out of basis
576              if (columnLowerWork_[iColumn]>-largeValue_) {
577                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
578                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
579                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
580                  setColumnStatus(iColumn,atLowerBound);
581                } else {
582                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
583                  setColumnStatus(iColumn,atUpperBound);
584                }
585              } else if (columnUpperWork_[iColumn]<largeValue_) {
586                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
587                setColumnStatus(iColumn,atUpperBound);
588              } else {
589                columnActivityWork_[iColumn]=0.0;
590                setColumnStatus(iColumn,isFree);
591              }
592            } else {
593              numberBasic++;
594            }
595            break;
596          case isFree:
597            columnActivityWork_[iColumn]=0.0;
598            break;
599          case atUpperBound:
600            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
601            if (columnActivityWork_[iColumn]>largeValue_) {
602              assert(columnLowerWork_[iColumn]>-largeValue_);
603              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
604              setColumnStatus(iColumn,atLowerBound);
605            }
606            break;
607          case atLowerBound:
608            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
609            if (columnActivityWork_[iColumn]<-largeValue_) {
610              assert(columnUpperWork_[iColumn]<largeValue_);
611              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
612              setColumnStatus(iColumn,atUpperBound);
613            }
614            break;
615          case superBasic:
616            if (columnUpperWork_[iColumn]>largeValue_) {
617              if (columnLowerWork_[iColumn]>-largeValue_) {
618                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
619                setColumnStatus(iColumn,atLowerBound);
620              } else {
621                // say free
622                setColumnStatus(iColumn,isFree);
623                columnActivityWork_[iColumn]=0.0;
624              }
625            } else {
626              if (columnLowerWork_[iColumn]>-largeValue_) {
627                // set to nearest
628                if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
629                    <fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])) {
630                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
631                  setColumnStatus(iColumn,atLowerBound);
632                } else {
633                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
634                  setColumnStatus(iColumn,atUpperBound);
635                }
636              } else {
637                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
638                setColumnStatus(iColumn,atUpperBound);
639              }
640            }
641            break;
642          }
643        }
644      } else {
645        //#define TESTFREE
646        // all slack basis
647        int numberBasic=0;
648        // changed to put free variables in basis
649        if (!status_) {
650          createStatus();
651        }
652        for (iRow=0;iRow<numberRows_;iRow++) {
653          double lower=rowLowerWork_[iRow];
654          double upper=rowUpperWork_[iRow];
655          if (lower>-largeValue_||upper<largeValue_) {
656            if (fabs(lower)<=fabs(upper)) {
657              setRowStatus(iRow,atLowerBound);
658              rowActivityWork_[iRow]=lower;
659            } else {
660              setRowStatus(iRow,atUpperBound);
661              rowActivityWork_[iRow]=upper;
662            }
663          } else {
664            setRowStatus(iRow,isFree);
665            rowActivityWork_[iRow]=0.0;
666          }
667#ifdef TESTFREE
668          setRowStatus(iRow,basic);
669          numberBasic++;
670#endif
671        }
672        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
673          double lower=columnLowerWork_[iColumn];
674          double upper=columnUpperWork_[iColumn];
675          double big_bound = largeValue_;
676          if (lower>-big_bound||upper<big_bound) {
677            if (fabs(lower)<=fabs(upper)) {
678              setColumnStatus(iColumn,atLowerBound);
679              columnActivityWork_[iColumn]=lower;
680            } else {
681              setColumnStatus(iColumn,atUpperBound);
682              columnActivityWork_[iColumn]=upper;
683            }
684          } else {
685#ifndef TESTFREE
686            numberBasic++;
687            setColumnStatus(iColumn,basic);
688#else
689            setColumnStatus(iColumn,isFree);
690#endif
691            columnActivityWork_[iColumn]=0.0;
692          }
693        }
694        assert(numberBasic<=numberRows_); // problems if too many free
695        if (!numberBasic) {
696          // might as well do all slack basis
697          for (iRow=0;iRow<numberRows_;iRow++) {
698            setRowStatus(iRow,basic);
699          }
700        }
701      }
702    } else {
703      // values pass has less coding
704      // make row activities correct
705      memset(rowActivityWork_,0,numberRows_*sizeof(double));
706      times(1.0,columnActivityWork_,rowActivityWork_);
707      if (status_) {
708        // set values from warm start (if sensible)
709        int numberBasic=0;
710        for (iRow=0;iRow<numberRows_;iRow++) {
711          if (getRowStatus(iRow)==basic) 
712            numberBasic++;
713          else {
714            setRowStatus(iRow,superBasic);
715            // but put to bound if close
716            if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
717                <=primalTolerance_)
718              rowActivityWork_[iRow]=rowLowerWork_[iRow];
719            else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
720                <=primalTolerance_)
721              rowActivityWork_[iRow]=rowUpperWork_[iRow];
722          }
723         
724        }
725        totalSlacks=numberBasic;
726        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
727          if (getColumnStatus(iColumn)==basic) {
728            if (numberBasic==numberRows_) {
729              // take out of basis
730              setColumnStatus(iColumn,superBasic);
731              // but put to bound if close
732              if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
733                  <=primalTolerance_)
734                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
735              else if (fabs(columnActivityWork_[iColumn]
736                            -columnUpperWork_[iColumn])
737                       <=primalTolerance_)
738                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
739            } else 
740              numberBasic++;
741          } else {
742            setColumnStatus(iColumn,superBasic);
743            // but put to bound if close
744            if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
745                <=primalTolerance_)
746              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
747            else if (fabs(columnActivityWork_[iColumn]
748                          -columnUpperWork_[iColumn])
749                     <=primalTolerance_)
750              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
751          }
752        }
753      } else {
754        // all slack basis
755        int numberBasic=0;
756        if (!status_) {
757          createStatus();
758        }
759        for (iRow=0;iRow<numberRows_;iRow++) {
760          setRowStatus(iRow,basic);
761          numberBasic++;
762        }
763        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
764          setColumnStatus(iColumn,superBasic);
765          // but put to bound if close
766          if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
767              <=primalTolerance_)
768            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
769          else if (fabs(columnActivityWork_[iColumn]
770                        -columnUpperWork_[iColumn])
771                   <=primalTolerance_)
772            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
773        }
774      }
775    }
776    numberRefinements_=1;
777  }
778  int status=-99;
779  int * rowIsBasic = new int[numberRows_];
780  int * columnIsBasic = new int[numberColumns_];
781  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
782  while (status<-98) {
783 
784    int i;
785    int numberBasic=0;
786    for (i=0;i<numberRows_;i++) {
787      if (getRowStatus(i) == basic) {
788        rowIsBasic[i]=1;
789        numberBasic++;
790      } else {
791        rowIsBasic[i]=-1;
792      }
793    }
794    for (i=0;i<numberColumns_;i++) {
795      if (getColumnStatus(i) == basic) {
796        columnIsBasic[i]=1;
797        numberBasic++;
798      } else {
799        columnIsBasic[i]=-1;
800      }
801    }
802    assert (numberBasic<=numberRows_);
803    while (status==-99) {
804      status =  factorization_->factorize(this,matrix_,
805                                          numberRows_,numberColumns_,
806                                          rowIsBasic, columnIsBasic,
807                                          0.0);
808      if (status==-99) {
809        // get more memory
810        factorization_->areaFactor(2.0*factorization_->areaFactor());
811      }
812    }
813    if (!status) {
814      // do pivot information
815      for (i=0;i<numberRows_;i++) {
816        if (getRowStatus(i) == basic) {
817          pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
818        }
819      }
820      for (i=0;i<numberColumns_;i++) {
821        if (getColumnStatus(i) == basic) {
822          pivotVariable_[columnIsBasic[i]]=i;
823        }
824      }
825    } else {
826      // leave pivotVariable_ in useful form for cleaning basis
827      for (i=0;i<numberRows_;i++) {
828        pivotVariable_[i]=-1;
829      }
830      for (i=0;i<numberRows_;i++) {
831        if (getRowStatus(i) == basic) {
832          int iPivot = rowIsBasic[i];
833          if (iPivot>=0) 
834            pivotVariable_[iPivot]=i+numberColumns_;
835        }
836      }
837      for (i=0;i<numberColumns_;i++) {
838        if (getColumnStatus(i) == basic) {
839          int iPivot = columnIsBasic[i];
840          if (iPivot>=0) 
841            pivotVariable_[iPivot]=i;
842        }
843      }
844    }
845    if (status==-1) {
846      if (!solveType) {
847        //redo basis - first take ALL columns out
848        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
849          if (getColumnStatus(iColumn)==
850              basic) {
851            // take out
852            if (!valuesPass) {
853              double lower=columnLowerWork_[iColumn];
854              double upper=columnUpperWork_[iColumn];
855              double value=columnActivityWork_[iColumn];
856              if (lower>-largeValue_||upper<largeValue_) {
857                if (fabs(value-lower)<fabs(value-upper)) {
858                  setColumnStatus(iColumn,atLowerBound);
859                  columnActivityWork_[iColumn]=lower;
860                } else {
861                  setColumnStatus(iColumn,atUpperBound);
862                  columnActivityWork_[iColumn]=upper;
863                }
864              } else {
865                setColumnStatus(iColumn,isFree);
866              }
867            } else {
868              setColumnStatus(iColumn,superBasic);
869            }
870          }
871        }
872        for (iRow=0;iRow<numberRows_;iRow++) {
873          int iSequence=pivotVariable_[iRow];
874          if (iSequence>=0) {
875            // basic
876            if (iSequence>=numberColumns_) {
877              // slack in - leave
878              assert (iSequence-numberColumns_==iRow);
879            } else {
880              // put back structural
881              setColumnStatus(iSequence,basic);
882            }
883          } else {
884            // put in slack
885            setRowStatus(iRow,basic);
886          }
887        }
888        // signal repeat
889        status=-99;
890      }
891    }
892  } 
893  delete [] rowIsBasic;
894  delete [] columnIsBasic;
895  if (status) {
896    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
897      <<status
898      <<CoinMessageEol;
899    return -1;
900  } else if (!solveType) {
901    // Initial basis - return number of singularities
902    int numberSlacks=0;
903    for (iRow=0;iRow<numberRows_;iRow++) {
904      if (getRowStatus(iRow) == basic)
905        numberSlacks++;
906    }
907    status= max(numberSlacks-totalSlacks,0);
908  }
909
910  // sparse methods
911  if (factorization_->sparseThreshold()) {
912    // get default value
913    factorization_->sparseThreshold(0);
914    factorization_->goSparse();
915  }
916 
917  return status;
918}
919/*
920   This does basis housekeeping and does values for in/out variables.
921   Can also decide to re-factorize
922*/
923int 
924ClpSimplex::housekeeping(double objectiveChange)
925{
926  numberIterations_++;
927  changeMade_++; // something has happened
928  // incoming variable
929
930  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
931    <<directionOut_
932    <<directionIn_<<theta_
933    <<dualOut_<<dualIn_<<alpha_
934    <<CoinMessageEol;
935  if (getStatus(sequenceIn_)==isFree) {
936    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
937      <<sequenceIn_
938      <<CoinMessageEol;
939  }
940  // change of incoming
941  if (algorithm_<0) {
942    dualOut_ /= alpha_;
943    dualOut_ *= -directionOut_;
944  }
945  char rowcol[]={'R','C'};
946  double cost = cost_[sequenceIn_];
947  double value=valueIn_;
948  if (pivotRow_>=0)
949    pivotVariable_[pivotRow_]=sequenceIn();
950  if (directionIn_==-1) {
951    // as if from upper bound
952    if (sequenceIn_!=sequenceOut_) {
953      // variable becoming basic
954      setStatus(sequenceIn_,basic);
955      if (algorithm_<0) {
956        value = upperIn_+dualOut_;
957        dj_[sequenceIn_]=0.0;
958      } else {
959        value = valueIn_-fabs(theta_);
960      }
961    } else {
962      value=lowerIn_;
963      setStatus(sequenceIn_, atLowerBound);
964    }
965  } else {
966    // as if from lower bound
967    if (sequenceIn_!=sequenceOut_) {
968      // variable becoming basic
969      setStatus(sequenceIn_,basic);
970      if (algorithm_<0) {
971        value = lowerIn_+dualOut_;
972        dj_[sequenceIn_]=0.0;
973      } else {
974        value = valueIn_+fabs(theta_);
975      }
976    } else {
977      value=upperIn_;
978      setStatus(sequenceIn_, atUpperBound);
979    }
980  }
981  if (algorithm_<0)
982    objectiveChange += cost*(value-valueIn_);
983  else
984    objectiveChange += dualIn_*(value-valueIn_);
985  solution_[sequenceIn_]=value;
986
987  // outgoing
988  if (sequenceIn_!=sequenceOut_) {
989    assert( getStatus(sequenceOut_)== basic);
990    if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
991      progressFlag_ |= 1; // making real progress
992    if (algorithm_<0) {
993      if (directionOut_>0) {
994        value = lowerOut_;
995        setStatus(sequenceOut_,atLowerBound);
996        dj_[sequenceOut_]=theta_;
997      } else {
998        value = upperOut_;
999        setStatus(sequenceOut_,atUpperBound);
1000        dj_[sequenceOut_]=-theta_;
1001      }
1002      solution_[sequenceOut_]=value;
1003    } else {
1004      if (directionOut_>0) {
1005        value = lowerOut_;
1006      } else {
1007        value = upperOut_;
1008      }
1009      double lowerValue = lower_[sequenceOut_];
1010      double upperValue = upper_[sequenceOut_];
1011      assert(value>=lowerValue-primalTolerance_&&
1012             value<=upperValue+primalTolerance_);
1013      // may not be exactly at bound and bounds may have changed
1014      if (value<=lowerValue+primalTolerance_) {
1015        setStatus(sequenceOut_,atLowerBound);
1016      } else if (value>=upperValue-primalTolerance_) {
1017        setStatus(sequenceOut_,atUpperBound);
1018      } else {
1019        printf("*** variable wandered off bound %g %g %g!\n",
1020               lowerValue,value,upperValue);
1021        if (value-lowerValue<=upperValue-value) {
1022          setStatus(sequenceOut_,atLowerBound);
1023          value=lowerValue;
1024        } else {
1025          setStatus(sequenceOut_,atUpperBound);
1026          value=upperValue;
1027        }
1028      }
1029      solution_[sequenceOut_]=value;
1030      nonLinearCost_->setOne(sequenceOut_,value);
1031    }
1032  }
1033  // change cost and bounds on incoming if primal
1034  if (algorithm_>0)
1035    nonLinearCost_->setOne(sequenceIn_,solution_[sequenceIn_]); 
1036  objectiveValue_ += objectiveChange;
1037  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
1038    <<numberIterations_<<objectiveValue_
1039    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
1040    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
1041  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
1042  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
1043  handler_->message()<<CoinMessageEol;
1044  if (numberIterations_>=maximumIterations_)
1045    return 2;
1046  // only time to re-factorize if one before real time
1047  // this is so user won't be surprised that maximumPivots has exact meaning
1048  if (factorization_->pivots()==factorization_->maximumPivots()) {
1049    return 1;
1050  } else {
1051    if (forceFactorization_>0&&
1052        factorization_->pivots()==forceFactorization_) {
1053      // relax
1054      forceFactorization_ *= 2;
1055      if (forceFactorization_>factorization_->maximumPivots())
1056        forceFactorization_ = -1; //off
1057      return 1;
1058    } else {
1059      // carry on iterating
1060      return 0;
1061    }
1062  }
1063}
1064// Copy constructor.
1065ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
1066  ClpModel(rhs),
1067  columnPrimalInfeasibility_(0.0),
1068  columnPrimalSequence_(-2),
1069  rowPrimalInfeasibility_(0.0),
1070  rowPrimalSequence_(-2), 
1071  columnDualInfeasibility_(0.0),
1072  columnDualSequence_(-2),
1073  rowDualInfeasibility_(0.0),
1074  rowDualSequence_(-2),
1075  primalToleranceToGetOptimal_(-1.0),
1076  remainingDualInfeasibility_(0.0),
1077  largeValue_(1.0e15),
1078  largestPrimalError_(0.0),
1079  largestDualError_(0.0),
1080  largestSolutionError_(0.0),
1081  dualBound_(1.0e7),
1082  lower_(NULL),
1083  rowLowerWork_(NULL),
1084  columnLowerWork_(NULL),
1085  upper_(NULL),
1086  rowUpperWork_(NULL),
1087  columnUpperWork_(NULL),
1088  cost_(NULL),
1089  rowObjectiveWork_(NULL),
1090  objectiveWork_(NULL),
1091  alpha_(0.0),
1092  theta_(0.0),
1093  lowerIn_(0.0),
1094  valueIn_(0.0),
1095  upperIn_(0.0),
1096  dualIn_(0.0),
1097  sequenceIn_(-1),
1098  directionIn_(-1),
1099  lowerOut_(-1),
1100  valueOut_(-1),
1101  upperOut_(-1),
1102  dualOut_(-1),
1103  sequenceOut_(-1),
1104  directionOut_(-1),
1105  pivotRow_(-1),
1106  dj_(NULL),
1107  rowReducedCost_(NULL),
1108  reducedCostWork_(NULL),
1109  solution_(NULL),
1110  rowActivityWork_(NULL),
1111  columnActivityWork_(NULL),
1112  dualTolerance_(0.0),
1113  primalTolerance_(0.0),
1114  sumDualInfeasibilities_(0.0),
1115  numberDualInfeasibilities_(0),
1116  numberDualInfeasibilitiesWithoutFree_(0),
1117  sumPrimalInfeasibilities_(0.0),
1118  numberPrimalInfeasibilities_(0),
1119  pivotVariable_(NULL),
1120  factorization_(NULL),
1121  numberRefinements_(0),
1122  rowScale_(NULL),
1123  savedSolution_(NULL),
1124  columnScale_(NULL),
1125  scalingFlag_(0),
1126  numberTimesOptimal_(0),
1127  changeMade_(1),
1128  algorithm_(0),
1129  forceFactorization_(-1),
1130  perturbation_(100),
1131  infeasibilityCost_(1.0e7),
1132  nonLinearCost_(NULL),
1133  specialOptions_(0),
1134  lastBadIteration_(-999999),
1135  numberFake_(0),
1136  progressFlag_(0),
1137  sumOfRelaxedDualInfeasibilities_(0.0),
1138  sumOfRelaxedPrimalInfeasibilities_(0.0)
1139{
1140  int i;
1141  for (i=0;i<6;i++) {
1142    rowArray_[i]=NULL;
1143    columnArray_[i]=NULL;
1144  }
1145  saveStatus_=NULL;
1146  factorization_ = NULL;
1147  dualRowPivot_ = NULL;
1148  primalColumnPivot_ = NULL;
1149  gutsOfDelete(0);
1150  gutsOfCopy(rhs);
1151}
1152// Copy constructor from model
1153ClpSimplex::ClpSimplex(const ClpModel &rhs) :
1154  ClpModel(rhs),
1155  columnPrimalInfeasibility_(0.0),
1156  columnPrimalSequence_(-2),
1157  rowPrimalInfeasibility_(0.0),
1158  rowPrimalSequence_(-2), 
1159  columnDualInfeasibility_(0.0),
1160  columnDualSequence_(-2),
1161  rowDualInfeasibility_(0.0),
1162  rowDualSequence_(-2),
1163  primalToleranceToGetOptimal_(-1.0),
1164  remainingDualInfeasibility_(0.0),
1165  largeValue_(1.0e15),
1166  largestPrimalError_(0.0),
1167  largestDualError_(0.0),
1168  largestSolutionError_(0.0),
1169  dualBound_(1.0e7),
1170  lower_(NULL),
1171  rowLowerWork_(NULL),
1172  columnLowerWork_(NULL),
1173  upper_(NULL),
1174  rowUpperWork_(NULL),
1175  columnUpperWork_(NULL),
1176  cost_(NULL),
1177  rowObjectiveWork_(NULL),
1178  objectiveWork_(NULL),
1179  alpha_(0.0),
1180  theta_(0.0),
1181  lowerIn_(0.0),
1182  valueIn_(0.0),
1183  upperIn_(0.0),
1184  dualIn_(0.0),
1185  sequenceIn_(-1),
1186  directionIn_(-1),
1187  lowerOut_(-1),
1188  valueOut_(-1),
1189  upperOut_(-1),
1190  dualOut_(-1),
1191  sequenceOut_(-1),
1192  directionOut_(-1),
1193  pivotRow_(-1),
1194  dj_(NULL),
1195  rowReducedCost_(NULL),
1196  reducedCostWork_(NULL),
1197  solution_(NULL),
1198  rowActivityWork_(NULL),
1199  columnActivityWork_(NULL),
1200  dualTolerance_(0.0),
1201  primalTolerance_(0.0),
1202  sumDualInfeasibilities_(0.0),
1203  numberDualInfeasibilities_(0),
1204  numberDualInfeasibilitiesWithoutFree_(0),
1205  sumPrimalInfeasibilities_(0.0),
1206  numberPrimalInfeasibilities_(0),
1207  pivotVariable_(NULL),
1208  factorization_(NULL),
1209  numberRefinements_(0),
1210  rowScale_(NULL),
1211  savedSolution_(NULL),
1212  columnScale_(NULL),
1213  scalingFlag_(0),
1214  numberTimesOptimal_(0),
1215  changeMade_(1),
1216  algorithm_(0),
1217  forceFactorization_(-1),
1218  perturbation_(100),
1219  infeasibilityCost_(1.0e7),
1220  nonLinearCost_(NULL),
1221  specialOptions_(0),
1222  lastBadIteration_(-999999),
1223  numberFake_(0),
1224  progressFlag_(0),
1225  sumOfRelaxedDualInfeasibilities_(0.0),
1226  sumOfRelaxedPrimalInfeasibilities_(0.0)
1227{
1228  int i;
1229  for (i=0;i<6;i++) {
1230    rowArray_[i]=NULL;
1231    columnArray_[i]=NULL;
1232  }
1233  saveStatus_=NULL;
1234  // get an empty factorization so we can set tolerances etc
1235  factorization_ = new ClpFactorization();
1236  // say Steepest pricing
1237  dualRowPivot_ = new ClpDualRowSteepest();
1238  // say Dantzig pricing
1239  primalColumnPivot_ = new ClpPrimalColumnDantzig();
1240 
1241}
1242// Assignment operator. This copies the data
1243ClpSimplex & 
1244ClpSimplex::operator=(const ClpSimplex & rhs)
1245{
1246  if (this != &rhs) {
1247    gutsOfDelete(0);
1248    ClpModel::operator=(rhs);
1249    gutsOfCopy(rhs);
1250  }
1251  return *this;
1252}
1253void 
1254ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1255{
1256  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
1257  rowLowerWork_ = lower_+numberColumns_;
1258  columnLowerWork_ = lower_;
1259  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
1260  rowUpperWork_ = upper_+numberColumns_;
1261  columnUpperWork_ = upper_;
1262  cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
1263  objectiveWork_ = cost_;
1264  rowObjectiveWork_ = cost_+numberColumns_;
1265  dj_ = ClpCopyOfArray(rhs.dj_,numberRows_+numberColumns_);
1266  if (dj_) {
1267    reducedCostWork_ = dj_;
1268    rowReducedCost_ = dj_+numberColumns_;
1269  }
1270  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
1271  if (solution_) {
1272    columnActivityWork_ = solution_;
1273    rowActivityWork_ = solution_+numberColumns_;
1274  }
1275  if (rhs.pivotVariable_) {
1276    pivotVariable_ = new int[numberRows_];
1277    ClpDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
1278  } else {
1279    pivotVariable_=NULL;
1280  }
1281  if (rhs.factorization_) {
1282    factorization_ = new ClpFactorization(*rhs.factorization_);
1283  } else {
1284    factorization_=NULL;
1285  }
1286  rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_);
1287  savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1288  columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_);
1289  int i;
1290  for (i=0;i<6;i++) {
1291    rowArray_[i]=NULL;
1292    if (rhs.rowArray_[i]) 
1293      rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
1294    columnArray_[i]=NULL;
1295    if (rhs.columnArray_[i]) 
1296      columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
1297  }
1298  if (rhs.saveStatus_) {
1299    saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1300  }
1301  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1302  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1303  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1304  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1305  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1306  columnDualSequence_ = rhs.columnDualSequence_;
1307  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1308  rowDualSequence_ = rhs.rowDualSequence_;
1309  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1310  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1311  largeValue_ = rhs.largeValue_;
1312  largestPrimalError_ = rhs.largestPrimalError_;
1313  largestDualError_ = rhs.largestDualError_;
1314  largestSolutionError_ = rhs.largestSolutionError_;
1315  dualBound_ = rhs.dualBound_;
1316  alpha_ = rhs.alpha_;
1317  theta_ = rhs.theta_;
1318  lowerIn_ = rhs.lowerIn_;
1319  valueIn_ = rhs.valueIn_;
1320  upperIn_ = rhs.upperIn_;
1321  dualIn_ = rhs.dualIn_;
1322  sequenceIn_ = rhs.sequenceIn_;
1323  directionIn_ = rhs.directionIn_;
1324  lowerOut_ = rhs.lowerOut_;
1325  valueOut_ = rhs.valueOut_;
1326  upperOut_ = rhs.upperOut_;
1327  dualOut_ = rhs.dualOut_;
1328  sequenceOut_ = rhs.sequenceOut_;
1329  directionOut_ = rhs.directionOut_;
1330  pivotRow_ = rhs.pivotRow_;
1331  numberRefinements_ = rhs.numberRefinements_;
1332  dualTolerance_ = rhs.dualTolerance_;
1333  primalTolerance_ = rhs.primalTolerance_;
1334  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1335  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1336  numberDualInfeasibilitiesWithoutFree_ = 
1337    rhs.numberDualInfeasibilitiesWithoutFree_;
1338  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1339  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1340  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1341  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1342  scalingFlag_ = rhs.scalingFlag_;
1343  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1344  changeMade_ = rhs.changeMade_;
1345  algorithm_ = rhs.algorithm_;
1346  forceFactorization_ = rhs.forceFactorization_;
1347  perturbation_ = rhs.perturbation_;
1348  infeasibilityCost_ = rhs.infeasibilityCost_;
1349  specialOptions_ = rhs.specialOptions_;
1350  lastBadIteration_ = rhs.lastBadIteration_;
1351  numberFake_ = rhs.numberFake_;
1352  progressFlag_ = rhs.progressFlag_;
1353  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
1354  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
1355  if (rhs.nonLinearCost_!=NULL)
1356    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1357}
1358// type == 0 do everything, most + pivot data, 2 factorization data as well
1359void 
1360ClpSimplex::gutsOfDelete(int type)
1361{
1362  delete [] lower_;
1363  lower_=NULL;
1364  rowLowerWork_=NULL;
1365  columnLowerWork_=NULL;
1366  delete [] upper_;
1367  upper_=NULL;
1368  rowUpperWork_=NULL;
1369  columnUpperWork_=NULL;
1370  delete [] cost_;
1371  cost_=NULL;
1372  objectiveWork_=NULL;
1373  rowObjectiveWork_=NULL;
1374  delete [] dj_;
1375  dj_=NULL;
1376  reducedCostWork_=NULL;
1377  rowReducedCost_=NULL;
1378  delete [] solution_;
1379  solution_=NULL;
1380  rowActivityWork_=NULL;
1381  columnActivityWork_=NULL;
1382  delete [] rowScale_;
1383  rowScale_ = NULL;
1384  delete [] savedSolution_;
1385  savedSolution_ = NULL;
1386  delete [] columnScale_;
1387  columnScale_ = NULL;
1388  delete nonLinearCost_;
1389  nonLinearCost_ = NULL;
1390  int i;
1391  for (i=0;i<6;i++) {
1392    delete rowArray_[i];
1393    rowArray_[i]=NULL;
1394    delete columnArray_[i];
1395    columnArray_[i]=NULL;
1396  }
1397  delete rowCopy_;
1398  rowCopy_=NULL;
1399  delete [] saveStatus_;
1400  saveStatus_=NULL;
1401  if (!type) {
1402    // delete everything
1403    delete factorization_;
1404    factorization_ = NULL;
1405    delete [] pivotVariable_;
1406    pivotVariable_=NULL;
1407    delete dualRowPivot_;
1408    dualRowPivot_ = NULL;
1409    delete primalColumnPivot_;
1410    primalColumnPivot_ = NULL;
1411  } else {
1412    // delete any size information in methods
1413    if (type>1) {
1414      factorization_->clearArrays();
1415      delete [] pivotVariable_;
1416      pivotVariable_=NULL;
1417    }
1418    dualRowPivot_->clearArrays();
1419    primalColumnPivot_->clearArrays();
1420  }
1421}
1422// This sets largest infeasibility and most infeasible
1423void 
1424ClpSimplex::checkPrimalSolution(const double * rowActivities,
1425                                        const double * columnActivities)
1426{
1427  double * solution;
1428  int iRow,iColumn;
1429
1430  objectiveValue_ = 0.0;
1431  // now look at primal solution
1432  columnPrimalInfeasibility_=0.0;
1433  columnPrimalSequence_=-1;
1434  rowPrimalInfeasibility_=0.0;
1435  rowPrimalSequence_=-1;
1436  largestSolutionError_=0.0;
1437  solution = rowActivityWork_;
1438  sumPrimalInfeasibilities_=0.0;
1439  numberPrimalInfeasibilities_=0;
1440  double primalTolerance = primalTolerance_;
1441  double relaxedTolerance=dualTolerance_;
1442  // we can't really trust infeasibilities if there is primal error
1443  double error = min(1.0e-3,largestPrimalError_);
1444  relaxedTolerance = max(relaxedTolerance, error);
1445  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
1446
1447  for (iRow=0;iRow<numberRows_;iRow++) {
1448    double infeasibility=0.0;
1449    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1450    if (solution[iRow]>rowUpperWork_[iRow]) {
1451      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1452    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1453      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1454    }
1455    if (infeasibility>primalTolerance) {
1456      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1457      if (infeasibility>relaxedTolerance) 
1458        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1459      numberPrimalInfeasibilities_ ++;
1460    }
1461    if (infeasibility>rowPrimalInfeasibility_) {
1462      rowPrimalInfeasibility_=infeasibility;
1463      rowPrimalSequence_=iRow;
1464    }
1465    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1466    if (infeasibility>largestSolutionError_)
1467      largestSolutionError_=infeasibility;
1468    if (rowLowerWork_[iRow]!=rowUpperWork_[iRow])
1469      clearFixed(iRow+numberColumns_);
1470    else
1471      setFixed(iRow+numberColumns_);
1472  }
1473  solution = columnActivityWork_;
1474  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1475    double infeasibility=0.0;
1476    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1477    if (solution[iColumn]>columnUpperWork_[iColumn]) {
1478      infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1479    } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1480      infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1481    }
1482    if (infeasibility>columnPrimalInfeasibility_) {
1483      columnPrimalInfeasibility_=infeasibility;
1484      columnPrimalSequence_=iColumn;
1485    }
1486    if (infeasibility>primalTolerance) {
1487      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1488      if (infeasibility>relaxedTolerance) 
1489        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1490      numberPrimalInfeasibilities_ ++;
1491    }
1492    infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1493    if (infeasibility>largestSolutionError_)
1494      largestSolutionError_=infeasibility;
1495    if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
1496      clearFixed(iColumn);
1497    else
1498      setFixed(iColumn);
1499  }
1500}
1501void 
1502ClpSimplex::checkDualSolution()
1503{
1504
1505  double * solution;
1506  int iRow,iColumn;
1507  sumDualInfeasibilities_=0.0;
1508  numberDualInfeasibilities_=0;
1509  numberDualInfeasibilitiesWithoutFree_=0;
1510  columnDualInfeasibility_=0.0;
1511  columnDualSequence_=-1;
1512  rowDualInfeasibility_=0.0;
1513  rowDualSequence_=-1;
1514  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
1515                                   columnPrimalInfeasibility_);
1516  remainingDualInfeasibility_=0.0;
1517  solution = rowActivityWork_;
1518  double relaxedTolerance=dualTolerance_;
1519  // we can't really trust infeasibilities if there is dual error
1520  double error = min(1.0e-3,largestDualError_);
1521  relaxedTolerance = max(relaxedTolerance, error);
1522  sumOfRelaxedDualInfeasibilities_ = 0.0;
1523
1524  for (iRow=0;iRow<numberRows_;iRow++) {
1525    if (getRowStatus(iRow) != basic) {
1526      // not basic
1527      double value = rowReducedCost_[iRow];
1528      double distance;
1529      distance = rowUpperWork_[iRow]-solution[iRow];
1530      if (distance>primalTolerance_) {
1531        // should not be negative
1532        if (value<0.0) {
1533          value = - value;
1534          if (value>rowDualInfeasibility_) {
1535            rowDualInfeasibility_=value;
1536            rowDualSequence_=iRow;
1537          }
1538          if (value>dualTolerance_) {
1539            sumDualInfeasibilities_ += value-dualTolerance_;
1540            if (value>relaxedTolerance) 
1541              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1542            numberDualInfeasibilities_ ++;
1543            if (getRowStatus(iRow) != isFree) 
1544              numberDualInfeasibilitiesWithoutFree_ ++;
1545            // maybe we can make feasible by increasing tolerance
1546            if (distance<largeValue_) {
1547              if (distance>primalToleranceToGetOptimal_)
1548                primalToleranceToGetOptimal_=distance;
1549            } else {
1550              //gap too big for any tolerance
1551              remainingDualInfeasibility_=
1552                max(remainingDualInfeasibility_,value);
1553            }
1554          }
1555        }
1556      }
1557      distance = solution[iRow] -rowLowerWork_[iRow];
1558      if (distance>primalTolerance_) {
1559        // should not be positive
1560        if (value>0.0) {
1561          if (value>rowDualInfeasibility_) {
1562            rowDualInfeasibility_=value;
1563            rowDualSequence_=iRow;
1564          }
1565          if (value>dualTolerance_) {
1566            sumDualInfeasibilities_ += value-dualTolerance_;
1567            if (value>relaxedTolerance) 
1568              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1569            numberDualInfeasibilities_ ++;
1570            if (getRowStatus(iRow) != isFree) 
1571              numberDualInfeasibilitiesWithoutFree_ ++;
1572            // maybe we can make feasible by increasing tolerance
1573            if (distance<largeValue_&&
1574                distance>primalToleranceToGetOptimal_)
1575              primalToleranceToGetOptimal_=distance;
1576          }
1577        }
1578      }
1579    }
1580  }
1581  solution = columnActivityWork_;
1582  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1583    if (getColumnStatus(iColumn) != basic) {
1584      // not basic
1585      double value = reducedCostWork_[iColumn];
1586      double distance;
1587      distance = columnUpperWork_[iColumn]-solution[iColumn];
1588      if (distance>primalTolerance_) {
1589        // should not be negative
1590        if (value<0.0) {
1591          value = - value;
1592          if (value>columnDualInfeasibility_) {
1593            columnDualInfeasibility_=value;
1594            columnDualSequence_=iColumn;
1595          }
1596          if (value>dualTolerance_) {
1597            sumDualInfeasibilities_ += value-dualTolerance_;
1598            if (value>relaxedTolerance) 
1599              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1600            numberDualInfeasibilities_ ++;
1601            if (getColumnStatus(iColumn) != isFree) 
1602              numberDualInfeasibilitiesWithoutFree_ ++;
1603            // maybe we can make feasible by increasing tolerance
1604            if (distance<largeValue_) {
1605              if (distance>primalToleranceToGetOptimal_)
1606                primalToleranceToGetOptimal_=distance;
1607            } else {
1608              //gap too big for any tolerance
1609              remainingDualInfeasibility_=
1610                max(remainingDualInfeasibility_,value);
1611            }
1612          }
1613        }
1614      }
1615      distance = solution[iColumn] -columnLowerWork_[iColumn];
1616      if (distance>primalTolerance_) {
1617        // should not be positive
1618        if (value>0.0) {
1619          if (value>columnDualInfeasibility_) {
1620            columnDualInfeasibility_=value;
1621            columnDualSequence_=iColumn;
1622          }
1623          if (value>dualTolerance_) {
1624            sumDualInfeasibilities_ += value-dualTolerance_;
1625            if (value>relaxedTolerance) 
1626              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1627            numberDualInfeasibilities_ ++;
1628            if (getColumnStatus(iColumn) != isFree) 
1629              numberDualInfeasibilitiesWithoutFree_ ++;
1630            // maybe we can make feasible by increasing tolerance
1631            if (distance<largeValue_&&
1632                distance>primalToleranceToGetOptimal_)
1633              primalToleranceToGetOptimal_=distance;
1634          }
1635        }
1636      }
1637    }
1638  }
1639}
1640/*
1641  Unpacks one column of the matrix into indexed array
1642*/
1643void 
1644ClpSimplex::unpack(CoinIndexedVector * rowArray)
1645{
1646  rowArray->clear();
1647  if (sequenceIn_>=numberColumns_) {
1648    //slack
1649    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
1650  } else {
1651    // column
1652    matrix_->unpack(this,rowArray,sequenceIn_);
1653  }
1654}
1655void 
1656ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence)
1657{
1658  rowArray->clear();
1659  if (sequence>=numberColumns_) {
1660    //slack
1661    rowArray->insert(sequence-numberColumns_,-1.0);
1662  } else {
1663    // column
1664    matrix_->unpack(this,rowArray,sequence);
1665  }
1666}
1667bool
1668ClpSimplex::createRim(int what,bool makeRowCopy)
1669{
1670  bool goodMatrix=true;
1671  int i;
1672  if ((what&16)!=0) {
1673    // move information to work arrays
1674    if (optimizationDirection_<0.0) {
1675      // reverse all dual signs
1676      for (i=0;i<numberColumns_;i++) 
1677        reducedCost_[i] = -reducedCost_[i];
1678      for (i=0;i<numberRows_;i++) 
1679        dual_[i] = -dual_[i];
1680    }
1681    // row reduced costs
1682    if (!dj_) {
1683      dj_ = new double[numberRows_+numberColumns_];
1684      reducedCostWork_ = dj_;
1685      rowReducedCost_ = dj_+numberColumns_;
1686      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1687      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
1688    }
1689    if (!solution_) {
1690      solution_ = new double[numberRows_+numberColumns_];
1691      columnActivityWork_ = solution_;
1692      rowActivityWork_ = solution_+numberColumns_;
1693      memcpy(columnActivityWork_,columnActivity_,
1694             numberColumns_*sizeof(double));
1695      memcpy(rowActivityWork_,rowActivity_,
1696             numberRows_*sizeof(double));
1697    }
1698    //check matrix
1699    if (!matrix_->allElementsInRange(this,1.0e-20,1.0e20)) {
1700      problemStatus_=4;
1701      goodMatrix= false;
1702    }
1703    if (makeRowCopy) {
1704      delete rowCopy_;
1705      // may return NULL if can't give row copy
1706      rowCopy_ = matrix_->reverseOrderedCopy();
1707    }
1708  }
1709  if ((what&4)!=0) {
1710    delete [] cost_;
1711    cost_ = new double[numberColumns_+numberRows_];
1712    objectiveWork_ = cost_;
1713    rowObjectiveWork_ = cost_+numberColumns_;
1714    memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
1715    if (rowObjective_)
1716      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1717    else
1718      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1719  }
1720  if ((what&1)!=0) {
1721    delete [] lower_;
1722    delete [] upper_;
1723    lower_ = new double[numberColumns_+numberRows_];
1724    upper_ = new double[numberColumns_+numberRows_];
1725    rowLowerWork_ = lower_+numberColumns_;
1726    columnLowerWork_ = lower_;
1727    rowUpperWork_ = upper_+numberColumns_;
1728    columnUpperWork_ = upper_;
1729    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1730    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1731    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1732    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1733    // clean up any mismatches on infinity
1734    for (i=0;i<numberColumns_;i++) {
1735      if (columnLowerWork_[i]<-1.0e30)
1736        columnLowerWork_[i] = -DBL_MAX;
1737      if (columnUpperWork_[i]>1.0e30)
1738        columnUpperWork_[i] = DBL_MAX;
1739    }
1740    // clean up any mismatches on infinity
1741    for (i=0;i<numberRows_;i++) {
1742      if (rowLowerWork_[i]<-1.0e30)
1743        rowLowerWork_[i] = -DBL_MAX;
1744      if (rowUpperWork_[i]>1.0e30)
1745        rowUpperWork_[i] = DBL_MAX;
1746    }
1747  }
1748  // do scaling if needed
1749  if (scalingFlag_>0&&!rowScale_) {
1750    if (matrix_->scale(this))
1751      scalingFlag_=-scalingFlag_; // not scaled after all
1752  }
1753  if ((what&4)!=0) {
1754    double direction = optimizationDirection_;
1755    // direction is actually scale out not scale in
1756    if (direction)
1757      direction = 1.0/direction;
1758    // but also scale by scale factors
1759    // not really sure about row scaling
1760    if (!rowScale_) {
1761      if (direction!=1.0) {
1762        for (i=0;i<numberRows_;i++)
1763          rowObjectiveWork_[i] *= direction;
1764        for (i=0;i<numberColumns_;i++)
1765          objectiveWork_[i] *= direction;
1766      }
1767    } else {
1768      for (i=0;i<numberRows_;i++)
1769        rowObjectiveWork_[i] *= direction/rowScale_[i];
1770      for (i=0;i<numberColumns_;i++)
1771        objectiveWork_[i] *= direction*columnScale_[i];
1772    }
1773  }
1774  if ((what&1)!=0&&rowScale_) {
1775    for (i=0;i<numberColumns_;i++) {
1776      double multiplier = 1.0/columnScale_[i];
1777      if (columnLowerWork_[i]>-1.0e50)
1778        columnLowerWork_[i] *= multiplier;
1779      if (columnUpperWork_[i]<1.0e50)
1780        columnUpperWork_[i] *= multiplier;
1781    }
1782    for (i=0;i<numberRows_;i++) {
1783      if (rowLowerWork_[i]>-1.0e50)
1784        rowLowerWork_[i] *= rowScale_[i];
1785      if (rowUpperWork_[i]<1.0e50)
1786        rowUpperWork_[i] *= rowScale_[i];
1787    }
1788  }
1789  if ((what&8)!=0&&rowScale_) {
1790    // on entry
1791    for (i=0;i<numberColumns_;i++) {
1792      columnActivityWork_[i] /= columnScale_[i];
1793      reducedCostWork_[i] *= columnScale_[i];
1794    }
1795    for (i=0;i<numberRows_;i++) {
1796      rowActivityWork_[i] *= rowScale_[i];
1797      dual_[i] /= rowScale_[i];
1798      rowReducedCost_[i] = dual_[i];
1799    }
1800  }
1801 
1802  if ((what&16)!=0) {
1803    // check rim of problem okay
1804    if (!sanityCheck())
1805      goodMatrix=false;
1806  } 
1807  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1808  // maybe we need to move scales to SimplexModel for factorization?
1809  if ((what&8)!=0&&!pivotVariable_) {
1810    pivotVariable_=new int[numberRows_];
1811  }
1812
1813  if ((what&16)!=0&&!rowArray_[2]) {
1814    // get some arrays
1815    int iRow,iColumn;
1816    // these are "indexed" arrays so we always know where nonzeros are
1817    /**********************************************************
1818      rowArray_[3] is long enough for columns as well
1819    *********************************************************/
1820    for (iRow=0;iRow<4;iRow++) {
1821      if (!rowArray_[iRow]) {
1822        rowArray_[iRow]=new CoinIndexedVector();
1823        int length =numberRows_+factorization_->maximumPivots();
1824        if (iRow==3)
1825          length = max(length,numberColumns_);
1826        rowArray_[iRow]->reserve(length);
1827      }
1828    }
1829   
1830    for (iColumn=0;iColumn<2;iColumn++) {
1831      if (!columnArray_[iColumn]) {
1832        columnArray_[iColumn]=new CoinIndexedVector();
1833        columnArray_[iColumn]->reserve(numberColumns_);
1834      }
1835    }   
1836  }
1837  return goodMatrix;
1838}
1839void
1840ClpSimplex::deleteRim(bool getRidOfFactorizationData)
1841{
1842  int i;
1843  if (problemStatus_!=1&&problemStatus_!=2) {
1844    delete [] ray_;
1845    ray_=NULL;
1846  }
1847  if (rowScale_) {
1848    for (i=0;i<numberColumns_;i++) {
1849      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1850      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1851    }
1852    for (i=0;i<numberRows_;i++) {
1853      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1854      dual_[i] *= rowScale_[i];
1855    }
1856    if (problemStatus_==2) {
1857      for (i=0;i<numberColumns_;i++) {
1858        ray_[i] *= columnScale_[i];
1859      }
1860    } else if (problemStatus_==1) {
1861      for (i=0;i<numberRows_;i++) {
1862        ray_[i] *= rowScale_[i];
1863      }
1864    }
1865  } else {
1866    for (i=0;i<numberColumns_;i++) {
1867      columnActivity_[i] = columnActivityWork_[i];
1868      reducedCost_[i] = reducedCostWork_[i];
1869    }
1870    for (i=0;i<numberRows_;i++) {
1871      rowActivity_[i] = rowActivityWork_[i];
1872    }
1873  }
1874  // direction may have been modified by scaling - clean up
1875  if (optimizationDirection_>0.0) {
1876    optimizationDirection_ = 1.0;
1877  }  else if (optimizationDirection_<0.0) {
1878    optimizationDirection_ = -1.0;
1879    // and reverse all dual signs
1880    for (i=0;i<numberColumns_;i++) 
1881      reducedCost_[i] = -reducedCost_[i];
1882    for (i=0;i<numberRows_;i++) 
1883      dual_[i] = -dual_[i];
1884  }
1885  // scaling may have been turned off
1886  scalingFlag_ = abs(scalingFlag_);
1887  if(getRidOfFactorizationData)
1888    gutsOfDelete(2);
1889  else
1890    gutsOfDelete(1);
1891}
1892void 
1893ClpSimplex::setDualBound(double value)
1894{
1895  if (value>0.0)
1896    dualBound_=value;
1897}
1898void 
1899ClpSimplex::setInfeasibilityCost(double value)
1900{
1901  if (value>0.0)
1902    infeasibilityCost_=value;
1903}
1904void ClpSimplex::setnumberRefinements( int value) 
1905{
1906  if (value>=0&&value<10)
1907    numberRefinements_=value;
1908}
1909// Sets row pivot choice algorithm in dual
1910void 
1911ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
1912{
1913  delete dualRowPivot_;
1914  dualRowPivot_ = choice.clone(true);
1915}
1916// Sets row pivot choice algorithm in dual
1917void 
1918ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
1919{
1920  delete primalColumnPivot_;
1921  primalColumnPivot_ = choice.clone(true);
1922}
1923// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
1924void 
1925ClpSimplex::scaling(int mode)
1926{
1927  if (mode>0&&mode<3) {
1928    scalingFlag_=mode;
1929  } else if (!mode) {
1930    scalingFlag_=0;
1931    delete [] rowScale_;
1932    rowScale_ = NULL;
1933    delete [] columnScale_;
1934    columnScale_ = NULL;
1935  }
1936}
1937// Passes in factorization
1938void 
1939ClpSimplex::setFactorization( ClpFactorization & factorization)
1940{
1941  delete factorization_;
1942  factorization_= new ClpFactorization(factorization);
1943}
1944void 
1945ClpSimplex::times(double scalar,
1946                  const double * x, double * y) const
1947{
1948  if (rowScale_)
1949    matrix_->times(scalar,x,y,rowScale_,columnScale_);
1950  else
1951    matrix_->times(scalar,x,y);
1952}
1953void 
1954ClpSimplex::transposeTimes(double scalar,
1955                           const double * x, double * y) const 
1956{
1957  if (rowScale_)
1958    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
1959  else
1960    matrix_->transposeTimes(scalar,x,y);
1961}
1962/* Perturbation:
1963   -50 to +50 - perturb by this power of ten (-6 sounds good)
1964   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1965   101 - we are perturbed
1966   102 - don't try perturbing again
1967   default is 100
1968*/
1969void 
1970ClpSimplex::setPerturbation(int value)
1971{
1972  if(value<=100&&value >=-50) {
1973    perturbation_=value;
1974  } 
1975}
1976// Sparsity on or off
1977bool 
1978ClpSimplex::sparseFactorization() const
1979{
1980  return factorization_->sparseThreshold();
1981}
1982void 
1983ClpSimplex::setSparseFactorization(bool value)
1984{
1985  if (value) {
1986    if (!factorization_->sparseThreshold())
1987      factorization_->goSparse();
1988  } else {
1989    factorization_->sparseThreshold(0);
1990  }
1991}
1992/* Tightens primal bounds to make dual faster.  Unless
1993   fixed, bounds are slightly looser than they could be.
1994   This is to make dual go faster and is probably not needed
1995   with a presolve.  Returns non-zero if problem infeasible
1996*/
1997int 
1998ClpSimplex::tightenPrimalBounds()
1999{
2000 
2001  // Get a row copy in standard format
2002  CoinPackedMatrix copy;
2003  copy.reverseOrderedCopyOf(*matrix());
2004  // get matrix data pointers
2005  const int * column = copy.getIndices();
2006  const CoinBigIndex * rowStart = copy.getVectorStarts();
2007  const int * rowLength = copy.getVectorLengths(); 
2008  const double * element = copy.getElements();
2009  int numberChanged=1,iPass=0;
2010  double large = largeValue(); // treat bounds > this as infinite
2011  int numberInfeasible=0;
2012  int totalTightened = 0;
2013
2014  double tolerance = primalTolerance();
2015
2016
2017  // Save column bounds
2018  double * saveLower = new double [numberColumns_];
2019  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
2020  double * saveUpper = new double [numberColumns_];
2021  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
2022#define MAXPASS 10 
2023
2024  // Loop round seeing if we can tighten bounds
2025  // Would be faster to have a stack of possible rows
2026  // and we put altered rows back on stack
2027
2028  int iRow, iColumn;
2029
2030  while(numberChanged) {
2031
2032    numberChanged = 0; // Bounds tightened this pass
2033   
2034    if (iPass==MAXPASS) break;
2035    iPass++;
2036   
2037    for (iRow = 0; iRow < numberRows_; iRow++) {
2038
2039      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
2040
2041        // possible row
2042        int infiniteUpper = 0;
2043        int infiniteLower = 0;
2044        double maximumUp = 0.0;
2045        double maximumDown = 0.0;
2046        double newBound;
2047        CoinBigIndex rStart = rowStart[iRow];
2048        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2049        CoinBigIndex j;
2050
2051        // Compute possible lower and upper ranges
2052
2053        for (j = rStart; j < rEnd; ++j) {
2054          double value=element[j];
2055          iColumn = column[j];
2056          if (value > 0.0) {
2057            if (columnUpper_[iColumn] >= large) {
2058              maximumUp = 1e31;
2059              ++infiniteUpper;
2060            } else {
2061              maximumUp += columnUpper_[iColumn] * value;
2062            }
2063            if (columnLower_[iColumn] <= -large) {
2064              maximumDown = -1e31;
2065              ++infiniteLower;
2066            } else {
2067              maximumDown += columnLower_[iColumn] * value;
2068            }
2069          } else if (value<0.0) {
2070            if (columnUpper_[iColumn] >= large) {
2071              maximumDown = -1e31;
2072              ++infiniteLower;
2073            } else {
2074              maximumDown += columnUpper_[iColumn] * value;
2075            }
2076            if (columnLower_[iColumn] <= -large) {
2077              maximumUp = 1e31;
2078              ++infiniteUpper;
2079            } else {
2080              maximumUp += columnLower_[iColumn] * value;
2081            }
2082          }
2083        }
2084        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2085            maximumDown >= rowLower_[iRow] - tolerance) {
2086
2087          // Row is redundant - make totally free
2088          rowLower_[iRow]=-DBL_MAX;
2089          rowUpper_[iRow]=DBL_MAX;
2090
2091        } else {
2092          if (maximumUp < rowLower_[iRow] -tolerance ||
2093              maximumDown > rowUpper_[iRow]+tolerance) {
2094            // problem is infeasible - exit at once
2095            numberInfeasible++;
2096            break;
2097          }
2098
2099          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2100            for (j = rStart; j < rEnd; ++j) {
2101              double value=element[j];
2102              iColumn = column[j];
2103              if (value > 0.0) {
2104                if (columnUpper_[iColumn] < large) {
2105                  newBound = columnUpper_[iColumn] + 
2106                    (rowLower_[iRow] - maximumUp) / value;
2107                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2108                    // Tighten the lower bound
2109
2110                    columnLower_[iColumn] = newBound;
2111                    ++numberChanged;
2112
2113                    // check infeasible (relaxed)
2114                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2115                        -100.0*tolerance) 
2116                      numberInfeasible++;
2117                    infiniteLower=1; // skip looking at other bound
2118                  }
2119                }
2120              } else {
2121                if (columnLower_[iColumn] > -large) {
2122                  newBound = columnLower_[iColumn] + 
2123                    (rowLower_[iRow] - maximumUp) / value;
2124                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2125                    // Tighten the upper bound
2126
2127                    columnUpper_[iColumn] = newBound;
2128                    ++numberChanged;
2129
2130                    // check infeasible (relaxed)
2131                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2132                        -100.0*tolerance) 
2133                      numberInfeasible++;
2134                    infiniteLower=1; // skip looking at other bound
2135                  }
2136                }
2137              }
2138            }
2139          }
2140         
2141          // Try other way
2142          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2143            for (j = rStart; j < rEnd; ++j) {
2144              double value=element[j];
2145              iColumn = column[j];
2146              if (value < 0.0) {
2147                if (columnUpper_[iColumn] < large) {
2148                  newBound = columnUpper_[iColumn] + 
2149                    (rowUpper_[iRow] - maximumDown) / value;
2150                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2151                    // Tighten the lower bound
2152
2153                    columnLower_[iColumn] = newBound;
2154                    ++numberChanged;
2155
2156                    // check infeasible (relaxed)
2157                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2158                        -100.0*tolerance) 
2159                      numberInfeasible++;
2160                  }
2161                } 
2162              } else {
2163                if (columnLower_[iColumn] > -large) {
2164                  newBound = columnLower_[iColumn] + 
2165                    (rowUpper_[iRow] - maximumDown) / value;
2166                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2167                    // Tighten the upper bound
2168
2169                    columnUpper_[iColumn] = newBound;
2170                    ++numberChanged;
2171
2172                    // check infeasible (relaxed)
2173                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2174                        -100.0*tolerance) 
2175                      numberInfeasible++;
2176                  }
2177                }
2178              }
2179            }
2180          }
2181        }
2182      }
2183    }
2184    totalTightened += numberChanged;
2185    if (numberInfeasible) break;
2186  }
2187  if (!numberInfeasible) {
2188    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2189      <<totalTightened
2190      <<CoinMessageEol;
2191    // Set bounds slightly loose
2192    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2193      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2194        // fix
2195        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2196          columnLower_[iColumn]=columnUpper_[iColumn];
2197        else
2198          columnUpper_[iColumn]=columnLower_[iColumn];
2199      } else {
2200        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2201          // relax a bit
2202          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2203                                      saveUpper[iColumn]);
2204        }
2205        if (columnLower_[iColumn]>saveLower[iColumn]) {
2206          // relax a bit
2207          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2208                                      saveLower[iColumn]);
2209        }
2210      }
2211    }
2212  } else {
2213    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2214      <<numberInfeasible
2215      <<CoinMessageEol;
2216    // restore column bounds
2217    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2218    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2219  }
2220  delete [] saveLower;
2221  delete [] saveUpper;
2222  return (numberInfeasible);
2223}
2224// dual
2225#include "ClpSimplexDual.hpp"
2226int ClpSimplex::dual ( )
2227{
2228  return ((ClpSimplexDual *) this)->dual();
2229}
2230// primal
2231#include "ClpSimplexPrimal.hpp"
2232int ClpSimplex::primal (int ifValuesPass )
2233{
2234  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2235}
2236/* For strong branching.  On input lower and upper are new bounds
2237   while on output they are objective function values (>1.0e50 infeasible).
2238   Return code is 0 if nothing interesting, -1 if infeasible both
2239   ways and +1 if infeasible one way (check values to see which one(s))
2240*/
2241int ClpSimplex::strongBranching(int numberVariables,const int * variables,
2242                                double * newLower, double * newUpper,
2243                                bool stopOnFirstInfeasible,
2244                                bool alwaysFinish)
2245{
2246  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
2247                                                    newLower,  newUpper,
2248                                                    stopOnFirstInfeasible,
2249                                                    alwaysFinish);
2250}
2251/* Borrow model.  This is so we dont have to copy large amounts
2252   of data around.  It assumes a derived class wants to overwrite
2253   an empty model with a real one - while it does an algorithm.
2254   This is same as ClpModel one, but sets scaling on etc. */
2255void 
2256ClpSimplex::borrowModel(ClpModel & otherModel) 
2257{
2258  ClpModel::borrowModel(otherModel);
2259  createStatus();
2260  scaling();
2261  ClpDualRowSteepest steep1;
2262  setDualRowPivotAlgorithm(steep1);
2263  ClpPrimalColumnSteepest steep2;
2264  setPrimalColumnPivotAlgorithm(steep2);
2265}
2266typedef struct {
2267  double optimizationDirection;
2268  double dblParam[ClpLastDblParam];
2269  double objectiveValue;
2270  double dualBound;
2271  double dualTolerance;
2272  double primalTolerance;
2273  double sumDualInfeasibilities;
2274  double sumPrimalInfeasibilities;
2275  double infeasibilityCost;
2276  int numberRows;
2277  int numberColumns;
2278  int intParam[ClpLastIntParam];
2279  int numberIterations;
2280  int problemStatus;
2281  int maximumIterations;
2282  int lengthNames;
2283  int numberDualInfeasibilities;
2284  int numberDualInfeasibilitiesWithoutFree;
2285  int numberPrimalInfeasibilities;
2286  int numberRefinements;
2287  int scalingFlag;
2288  int algorithm;
2289  unsigned int specialOptions;
2290  int dualPivotChoice;
2291  int primalPivotChoice;
2292  int matrixStorageChoice;
2293} Clp_scalars;
2294
2295int outDoubleArray(double * array, int length, FILE * fp)
2296{
2297  int numberWritten;
2298  if (array&&length) {
2299    numberWritten = fwrite(&length,sizeof(int),1,fp);
2300    if (numberWritten!=1)
2301      return 1;
2302    numberWritten = fwrite(array,sizeof(double),length,fp);
2303    if (numberWritten!=length)
2304      return 1;
2305  } else {
2306    length = 0;
2307    numberWritten = fwrite(&length,sizeof(int),1,fp);
2308    if (numberWritten!=1)
2309      return 1;
2310  }
2311  return 0;
2312}
2313// Save model to file, returns 0 if success
2314int
2315ClpSimplex::saveModel(const char * fileName)
2316{
2317  FILE * fp = fopen(fileName,"wb");
2318  if (fp) {
2319    Clp_scalars scalars;
2320    int i;
2321    CoinBigIndex numberWritten;
2322    // Fill in scalars
2323    scalars.optimizationDirection = optimizationDirection_;
2324    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2325    scalars.objectiveValue = objectiveValue_;
2326    scalars.dualBound = dualBound_;
2327    scalars.dualTolerance = dualTolerance_;
2328    scalars.primalTolerance = primalTolerance_;
2329    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2330    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2331    scalars.infeasibilityCost = infeasibilityCost_;
2332    scalars.numberRows = numberRows_;
2333    scalars.numberColumns = numberColumns_;
2334    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2335    scalars.numberIterations = numberIterations_;
2336    scalars.problemStatus = problemStatus_;
2337    scalars.maximumIterations = maximumIterations_;
2338    scalars.lengthNames = lengthNames_;
2339    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2340    scalars.numberDualInfeasibilitiesWithoutFree
2341      = numberDualInfeasibilitiesWithoutFree_;
2342    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2343    scalars.numberRefinements = numberRefinements_;
2344    scalars.scalingFlag = scalingFlag_;
2345    scalars.algorithm = algorithm_;
2346    scalars.specialOptions = specialOptions_;
2347    scalars.dualPivotChoice = dualRowPivot_->type();
2348    scalars.primalPivotChoice = primalColumnPivot_->type();
2349    scalars.matrixStorageChoice = matrix_->type();
2350
2351    // put out scalars
2352    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2353    if (numberWritten!=1)
2354      return 1;
2355    // strings
2356    CoinBigIndex length;
2357    for (i=0;i<ClpLastStrParam;i++) {
2358      length = strParam_[i].size();
2359      numberWritten = fwrite(&length,sizeof(int),1,fp);
2360      if (numberWritten!=1)
2361        return 1;
2362      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2363      if (numberWritten!=1)
2364        return 1;
2365    }
2366    // arrays - in no particular order
2367    if (outDoubleArray(rowActivity_,numberRows_,fp))
2368        return 1;
2369    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2370        return 1;
2371    if (outDoubleArray(dual_,numberRows_,fp))
2372        return 1;
2373    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2374        return 1;
2375    if (outDoubleArray(rowLower_,numberRows_,fp))
2376        return 1;
2377    if (outDoubleArray(rowUpper_,numberRows_,fp))
2378        return 1;
2379    if (outDoubleArray(objective_,numberColumns_,fp))
2380        return 1;
2381    if (outDoubleArray(rowObjective_,numberRows_,fp))
2382        return 1;
2383    if (outDoubleArray(columnLower_,numberColumns_,fp))
2384        return 1;
2385    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2386        return 1;
2387    if (ray_) {
2388      if (problemStatus_==1)
2389        if (outDoubleArray(ray_,numberRows_,fp))
2390          return 1;
2391      else if (problemStatus_==2)
2392        if (outDoubleArray(ray_,numberColumns_,fp))
2393          return 1;
2394      else
2395        if (outDoubleArray(NULL,0,fp))
2396          return 1;
2397    } else {
2398      if (outDoubleArray(NULL,0,fp))
2399        return 1;
2400    }
2401    if (status_&&(numberRows_+numberColumns_)>0) {
2402      length = numberRows_+numberColumns_;
2403      numberWritten = fwrite(&length,sizeof(int),1,fp);
2404      if (numberWritten!=1)
2405        return 1;
2406      numberWritten = fwrite(status_,sizeof(char),length, fp);
2407      if (numberWritten!=length)
2408        return 1;
2409    } else {
2410      length = 0;
2411      numberWritten = fwrite(&length,sizeof(int),1,fp);
2412      if (numberWritten!=1)
2413        return 1;
2414    }
2415    if (lengthNames_) {
2416      assert (numberRows_ == (int) rowNames_.size());
2417      for (i=0;i<numberRows_;i++) {
2418        length = rowNames_[i].size();
2419        numberWritten = fwrite(&length,sizeof(int),1,fp);
2420        if (numberWritten!=1)
2421          return 1;
2422        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
2423        if (numberWritten!=1)
2424          return 1;
2425      }
2426      assert (numberColumns_ == (int) columnNames_.size());
2427      for (i=0;i<numberColumns_;i++) {
2428        length = columnNames_[i].size();
2429        numberWritten = fwrite(&length,sizeof(int),1,fp);
2430        if (numberWritten!=1)
2431          return 1;
2432        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
2433        if (numberWritten!=1)
2434          return 1;
2435      }
2436    }
2437    // just standard type at present
2438    assert (matrix_->type()==1);
2439    assert (matrix_->getNumCols() == numberColumns_);
2440    assert (matrix_->getNumRows() == numberRows_);
2441    // we are going to save with gaps
2442    length = matrix_->getVectorStarts()[numberColumns_-1]
2443      + matrix_->getVectorLengths()[numberColumns_-1];
2444    numberWritten = fwrite(&length,sizeof(int),1,fp);
2445    if (numberWritten!=1)
2446      return 1;
2447    numberWritten = fwrite(matrix_->getElements(),
2448                           sizeof(double),length,fp);
2449    if (numberWritten!=length)
2450      return 1;
2451    numberWritten = fwrite(matrix_->getIndices(),
2452                           sizeof(int),length,fp);
2453    if (numberWritten!=length)
2454      return 1;
2455    numberWritten = fwrite(matrix_->getVectorStarts(),
2456                           sizeof(int),numberColumns_,fp);
2457    if (numberWritten!=numberColumns_)
2458      return 1;
2459    numberWritten = fwrite(matrix_->getVectorLengths(),
2460                           sizeof(int),numberColumns_,fp);
2461    if (numberWritten!=numberColumns_)
2462      return 1;
2463    // finished
2464    fclose(fp);
2465    return 0;
2466  } else {
2467    return -1;
2468  }
2469}
2470
2471int inDoubleArray(double * &array, int length, FILE * fp)
2472{
2473  int numberRead;
2474  int length2;
2475  numberRead = fread(&length2,sizeof(int),1,fp);
2476  if (numberRead!=1)
2477    return 1;
2478  if (length2) {
2479    // lengths must match
2480    if (length!=length2)
2481      return 2;
2482    array = new double[length];
2483    numberRead = fread(array,sizeof(double),length,fp);
2484    if (numberRead!=length)
2485      return 1;
2486  } 
2487  return 0;
2488}
2489/* Restore model from file, returns 0 if success,
2490   deletes current model */
2491int 
2492ClpSimplex::restoreModel(const char * fileName)
2493{
2494  FILE * fp = fopen(fileName,"rb");
2495  if (fp) {
2496    // Get rid of current model
2497    ClpModel::gutsOfDelete();
2498    gutsOfDelete(0);
2499    int i;
2500    for (i=0;i<6;i++) {
2501      rowArray_[i]=NULL;
2502      columnArray_[i]=NULL;
2503    }
2504    // get an empty factorization so we can set tolerances etc
2505    factorization_ = new ClpFactorization();
2506    // Say sparse
2507    factorization_->sparseThreshold(1);
2508    Clp_scalars scalars;
2509    CoinBigIndex numberRead;
2510
2511    // get scalars
2512    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
2513    if (numberRead!=1)
2514      return 1;
2515    // Fill in scalars
2516    optimizationDirection_ = scalars.optimizationDirection;
2517    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
2518    objectiveValue_ = scalars.objectiveValue;
2519    dualBound_ = scalars.dualBound;
2520    dualTolerance_ = scalars.dualTolerance;
2521    primalTolerance_ = scalars.primalTolerance;
2522    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
2523    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
2524    infeasibilityCost_ = scalars.infeasibilityCost;
2525    numberRows_ = scalars.numberRows;
2526    numberColumns_ = scalars.numberColumns;
2527    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
2528    numberIterations_ = scalars.numberIterations;
2529    problemStatus_ = scalars.problemStatus;
2530    maximumIterations_ = scalars.maximumIterations;
2531    lengthNames_ = scalars.lengthNames;
2532    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
2533    numberDualInfeasibilitiesWithoutFree_
2534      = scalars.numberDualInfeasibilitiesWithoutFree;
2535    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
2536    numberRefinements_ = scalars.numberRefinements;
2537    scalingFlag_ = scalars.scalingFlag;
2538    algorithm_ = scalars.algorithm;
2539    specialOptions_ = scalars.specialOptions;
2540    // strings
2541    CoinBigIndex length;
2542    for (i=0;i<ClpLastStrParam;i++) {
2543      numberRead = fread(&length,sizeof(int),1,fp);
2544      if (numberRead!=1)
2545        return 1;
2546      char * array = new char[length+1];
2547      numberRead = fread(array,length,1,fp);
2548      if (numberRead!=1)
2549        return 1;
2550      array[length]='\0';
2551      strParam_[i]=array;
2552      delete [] array;
2553    }
2554    // arrays - in no particular order
2555    if (inDoubleArray(rowActivity_,numberRows_,fp))
2556        return 1;
2557    if (inDoubleArray(columnActivity_,numberColumns_,fp))
2558        return 1;
2559    if (inDoubleArray(dual_,numberRows_,fp))
2560        return 1;
2561    if (inDoubleArray(reducedCost_,numberColumns_,fp))
2562        return 1;
2563    if (inDoubleArray(rowLower_,numberRows_,fp))
2564        return 1;
2565    if (inDoubleArray(rowUpper_,numberRows_,fp))
2566        return 1;
2567    if (inDoubleArray(objective_,numberColumns_,fp))
2568        return 1;
2569    if (inDoubleArray(rowObjective_,numberRows_,fp))
2570        return 1;
2571    if (inDoubleArray(columnLower_,numberColumns_,fp))
2572        return 1;
2573    if (inDoubleArray(columnUpper_,numberColumns_,fp))
2574        return 1;
2575    if (problemStatus_==1) {
2576      if (inDoubleArray(ray_,numberRows_,fp))
2577        return 1;
2578    } else if (problemStatus_==2) {
2579      if (inDoubleArray(ray_,numberColumns_,fp))
2580        return 1;
2581    } else {
2582      // ray should be null
2583      numberRead = fread(&length,sizeof(int),1,fp);
2584      if (numberRead!=1)
2585        return 1;
2586      if (length)
2587        return 2;
2588    }
2589    delete [] status_;
2590    status_=NULL;
2591    // status region
2592    numberRead = fread(&length,sizeof(int),1,fp);
2593    if (numberRead!=1)
2594        return 1;
2595    if (length) {
2596      if (length!=numberRows_+numberColumns_)
2597        return 1;
2598      status_ = new char unsigned[length];
2599      numberRead = fread(status_,sizeof(char),length, fp);
2600      if (numberRead!=length)
2601        return 1;
2602    }
2603    if (lengthNames_) {
2604      char * array = new char[lengthNames_+1];
2605      rowNames_.resize(0);
2606      for (i=0;i<numberRows_;i++) {
2607        numberRead = fread(&length,sizeof(int),1,fp);
2608        if (numberRead!=1)
2609          return 1;
2610        numberRead = fread(array,length,1,fp);
2611        if (numberRead!=1)
2612          return 1;
2613        rowNames_[i]=array;
2614      }
2615      columnNames_.resize(0);
2616      for (i=0;i<numberColumns_;i++) {
2617        numberRead = fread(&length,sizeof(int),1,fp);
2618        if (numberRead!=1)
2619          return 1;
2620        numberRead = fread(array,length,1,fp);
2621        if (numberRead!=1)
2622          return 1;
2623        columnNames_[i]=array;
2624      }
2625      delete [] array;
2626    }
2627    // Pivot choices
2628    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
2629    delete dualRowPivot_;
2630    switch ((scalars.dualPivotChoice&63)) {
2631    case 1:
2632      // Dantzig
2633      dualRowPivot_ = new ClpDualRowDantzig();
2634      break;
2635    case 2:
2636      // Steepest - use mode
2637      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
2638      break;
2639    default:
2640      abort();
2641    }
2642    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
2643    delete primalColumnPivot_;
2644    switch ((scalars.primalPivotChoice&63)) {
2645    case 1:
2646      // Dantzig
2647      primalColumnPivot_ = new ClpPrimalColumnDantzig();
2648      break;
2649    case 2:
2650      // Steepest - use mode
2651      primalColumnPivot_
2652        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
2653      break;
2654    default:
2655      abort();
2656    }
2657    assert(scalars.matrixStorageChoice==1);
2658    delete matrix_;
2659    // get arrays
2660    numberRead = fread(&length,sizeof(int),1,fp);
2661    if (numberRead!=1)
2662      return 1;
2663    double * elements = new double[length];
2664    int * indices = new int[length];
2665    CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
2666    int * lengths = new int[numberColumns_];
2667    numberRead = fread(elements, sizeof(double),length,fp);
2668    if (numberRead!=length)
2669      return 1;
2670    numberRead = fread(indices, sizeof(int),length,fp);
2671    if (numberRead!=length)
2672      return 1;
2673    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
2674    if (numberRead!=numberColumns_)
2675      return 1;
2676    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
2677    if (numberRead!=numberColumns_)
2678      return 1;
2679    // assign matrix
2680    CoinPackedMatrix * matrix = new CoinPackedMatrix();
2681    matrix->assignMatrix(true, numberRows_, numberColumns_,
2682                         length, elements, indices, starts, lengths);
2683    // and transfer to Clp
2684    matrix_ = new ClpPackedMatrix(matrix);
2685    // finished
2686    fclose(fp);
2687    return 0;
2688  } else {
2689    return -1;
2690  }
2691  return 0;
2692}
2693// value of incoming variable (in Dual)
2694double 
2695ClpSimplex::valueIncomingDual() const
2696{
2697  // Need value of incoming for list of infeasibilities as may be infeasible
2698  double valueIncoming = (dualOut_/alpha_)*directionOut_;
2699  if (directionIn_==-1)
2700    valueIncoming = upperIn_-valueIncoming;
2701  else
2702    valueIncoming = lowerIn_-valueIncoming;
2703  return valueIncoming;
2704}
2705//#############################################################################
2706
2707ClpSimplexProgress::ClpSimplexProgress () 
2708{
2709  int i;
2710  for (i=0;i<CLP_PROGRESS;i++) {
2711    objective_[i] = 0.0;
2712    infeasibility_[i] = -1.0; // set to an impossible value
2713    numberInfeasibilities_[i]=-1; 
2714  }
2715  numberTimes_ = 0;
2716  numberBadTimes_ = 0;
2717  model_ = NULL;
2718}
2719
2720
2721//-----------------------------------------------------------------------------
2722
2723ClpSimplexProgress::~ClpSimplexProgress ()
2724{
2725}
2726// Copy constructor.
2727ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2728{
2729  int i;
2730  for (i=0;i<CLP_PROGRESS;i++) {
2731    objective_[i] = rhs.objective_[i];
2732    infeasibility_[i] = rhs.infeasibility_[i];
2733    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2734  }
2735  numberTimes_ = rhs.numberTimes_;
2736  numberBadTimes_ = rhs.numberBadTimes_;
2737  model_ = rhs.model_;
2738}
2739// Copy constructor.from model
2740ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2741{
2742  int i;
2743  for (i=0;i<CLP_PROGRESS;i++) {
2744    objective_[i] = 0.0;
2745    infeasibility_[i] = -1.0; // set to an impossible value
2746    numberInfeasibilities_[i]=-1; 
2747  }
2748  numberTimes_ = 0;
2749  numberBadTimes_ = 0;
2750  model_ = model;
2751}
2752// Assignment operator. This copies the data
2753ClpSimplexProgress & 
2754ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2755{
2756  if (this != &rhs) {
2757    int i;
2758    for (i=0;i<CLP_PROGRESS;i++) {
2759      objective_[i] = rhs.objective_[i];
2760      infeasibility_[i] = rhs.infeasibility_[i];
2761      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2762    }
2763    numberTimes_ = rhs.numberTimes_;
2764    numberBadTimes_ = rhs.numberBadTimes_;
2765    model_ = rhs.model_;
2766  }
2767  return *this;
2768}
2769// Seems to be something odd about exact comparison of doubles on linux
2770static bool equalDouble(double value1, double value2)
2771{
2772  assert(sizeof(unsigned int)*2==sizeof(double));
2773  unsigned int *i1 = (unsigned int *) &value1;
2774  unsigned int *i2 = (unsigned int *) &value2;
2775  return (i1[0]==i2[0]&&i1[1]==i2[1]);
2776}
2777int
2778ClpSimplexProgress::looping()
2779{
2780  assert(model_);
2781  double objective = model_->objectiveValue();
2782  double infeasibility;
2783  int numberInfeasibilities;
2784  if (model_->algorithm()<0) {
2785    // dual
2786    infeasibility = model_->sumPrimalInfeasibilities();
2787    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2788  } else {
2789    //primal
2790    infeasibility = model_->sumDualInfeasibilities();
2791    numberInfeasibilities = model_->numberDualInfeasibilities();
2792  }
2793  int i;
2794  int numberMatched=0;
2795  int matched=0;
2796
2797  for (i=0;i<CLP_PROGRESS;i++) {
2798    bool matchedOnObjective = equalDouble(objective,objective_[i]);
2799    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
2800    bool matchedOnInfeasibilities = 
2801      (numberInfeasibilities==numberInfeasibilities_[i]);
2802
2803    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
2804      matched |= (1<<i);
2805      numberMatched++;
2806      // here mainly to get over compiler bug?
2807      if (model_->messageHandler()->logLevel()>10)
2808        printf("%d %d %d %d %d loop check\n",i,numberMatched,
2809               matchedOnObjective, matchedOnInfeasibility, 
2810               matchedOnInfeasibilities);
2811    }
2812    if (i) {
2813      objective_[i-1] = objective_[i];
2814      infeasibility_[i-1] = infeasibility_[i];
2815      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2816    }
2817  }
2818  objective_[CLP_PROGRESS-1] = objective;
2819  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2820  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2821  if (model_->progressFlag())
2822    numberMatched=0;
2823  numberTimes_++;
2824  if (numberTimes_<10)
2825    numberMatched=0;
2826  // skip if just last time as may be checking something
2827  if (matched==(1<<(CLP_PROGRESS-1)))
2828    numberMatched=0;
2829  if (numberMatched) {
2830    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2831      <<numberMatched
2832      <<matched
2833      <<numberTimes_
2834      <<CoinMessageEol;
2835    numberBadTimes_++;
2836    if (numberBadTimes_<10) {
2837      if (model_->algorithm()<0) {
2838        // dual - change tolerance
2839        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2840        // if infeasible increase dual bound
2841        if (model_->dualBound()<1.0e19) {
2842          model_->setDualBound(model_->dualBound()*1.1);
2843        }
2844      } else {
2845        // primal - change tolerance
2846        model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2847        // if infeasible increase infeasibility cost
2848        if (model_->infeasibilityCost()<1.0e19) {
2849          model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
2850        }
2851      }
2852    } else {
2853      model_->messageHandler()->message(CLP_LOOP,model_->messages())
2854        <<CoinMessageEol;
2855      // look at solution and maybe declare victory
2856      if (infeasibility<1.0e-4)
2857        return 0;
2858      else
2859        return 3;
2860    }
2861  }
2862  return -1;
2863}
2864// Sanity check on input data - returns true if okay
2865bool 
2866ClpSimplex::sanityCheck()
2867{
2868  // bad if empty
2869  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
2870    handler_->message(CLP_EMPTY_PROBLEM,messages_)
2871      <<numberRows_
2872      <<numberColumns_
2873      <<matrix_->getNumElements()
2874      <<CoinMessageEol;
2875    problemStatus_=4;
2876    return false;
2877  }
2878  int numberBad , numberBadBounds;
2879  double largestBound, smallestBound, minimumGap;
2880  double smallestObj, largestObj;
2881  int firstBad;
2882  int modifiedBounds=0;
2883  int i;
2884  numberBad=0;
2885  numberBadBounds=0;
2886  firstBad=-1;
2887  minimumGap=1.0e100;
2888  smallestBound=1.0e100;
2889  largestBound=0.0;
2890  smallestObj=1.0e100;
2891  largestObj=0.0;
2892  // If bounds are too close - fix
2893  double fixTolerance = 1.1*primalTolerance_;
2894  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
2895    double value;
2896    value = fabs(cost_[i]);
2897    if (value>1.0e50) {
2898      numberBad++;
2899      if (firstBad<0)
2900        firstBad=i;
2901    } else if (value) {
2902      if (value>largestObj)
2903        largestObj=value;
2904      if (value<smallestObj)
2905        smallestObj=value;
2906    }
2907    value=upper_[i]-lower_[i];
2908    if (value<-primalTolerance_) {
2909      numberBadBounds++;
2910      if (firstBad<0)
2911        firstBad=i;
2912    } else if (value<=fixTolerance) {
2913      if (value) {
2914        // modify
2915        upper_[i] = lower_[i];
2916        modifiedBounds++;
2917      }
2918    } else {
2919      if (value<minimumGap)
2920        minimumGap=value;
2921    }
2922    if (lower_[i]>-1.0e100&&lower_[i]) {
2923      value = fabs(lower_[i]);
2924      if (value>largestBound)
2925        largestBound=value;
2926      if (value<smallestBound)
2927        smallestBound=value;
2928    }
2929    if (upper_[i]<1.0e100&&upper_[i]) {
2930      value = fabs(upper_[i]);
2931      if (value>largestBound)
2932        largestBound=value;
2933      if (value<smallestBound)
2934        smallestBound=value;
2935    }
2936  }
2937  if (largestBound)
2938    handler_->message(CLP_RIMSTATISTICS3,messages_)
2939      <<smallestBound
2940      <<largestBound
2941      <<minimumGap
2942      <<CoinMessageEol;
2943  minimumGap=1.0e100;
2944  smallestBound=1.0e100;
2945  largestBound=0.0;
2946  for (i=0;i<numberColumns_;i++) {
2947    double value;
2948    value = fabs(cost_[i]);
2949    if (value>1.0e50) {
2950      numberBad++;
2951      if (firstBad<0)
2952        firstBad=i;
2953    } else if (value) {
2954      if (value>largestObj)
2955        largestObj=value;
2956      if (value<smallestObj)
2957        smallestObj=value;
2958    }
2959    value=upper_[i]-lower_[i];
2960    if (value<-primalTolerance_) {
2961      numberBadBounds++;
2962      if (firstBad<0)
2963        firstBad=i;
2964    } else if (value<=fixTolerance) {
2965      if (value) {
2966        // modify
2967        upper_[i] = lower_[i];
2968        modifiedBounds++;
2969      }
2970    } else {
2971      if (value<minimumGap)
2972        minimumGap=value;
2973    }
2974    if (lower_[i]>-1.0e100&&lower_[i]) {
2975      value = fabs(lower_[i]);
2976      if (value>largestBound)
2977        largestBound=value;
2978      if (value<smallestBound)
2979        smallestBound=value;
2980    }
2981    if (upper_[i]<1.0e100&&upper_[i]) {
2982      value = fabs(upper_[i]);
2983      if (value>largestBound)
2984        largestBound=value;
2985      if (value<smallestBound)
2986        smallestBound=value;
2987    }
2988  }
2989  char rowcol[]={'R','C'};
2990  if (numberBad) {
2991    handler_->message(CLP_BAD_BOUNDS,messages_)
2992      <<numberBad
2993      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
2994      <<CoinMessageEol;
2995    problemStatus_=4;
2996    return false;
2997  }
2998  if (modifiedBounds)
2999    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
3000      <<modifiedBounds
3001      <<CoinMessageEol;
3002  handler_->message(CLP_RIMSTATISTICS1,messages_)
3003    <<smallestObj
3004    <<largestObj
3005    <<CoinMessageEol;
3006  if (largestBound)
3007    handler_->message(CLP_RIMSTATISTICS2,messages_)
3008      <<smallestBound
3009      <<largestBound
3010      <<minimumGap
3011      <<CoinMessageEol;
3012  return true;
3013}
3014// Set up status array (for OsiClp)
3015void 
3016ClpSimplex::createStatus() 
3017{
3018  if(!status_)
3019    status_ = new unsigned char [numberColumns_+numberRows_];
3020  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3021  int i;
3022  // set column status to one nearest zero
3023  for (i=0;i<numberColumns_;i++) {
3024#if 0
3025    if (columnLower_[i]>=0.0) {
3026      setColumnStatus(i,atLowerBound);
3027    } else if (columnUpper_[i]<=0.0) {
3028      setColumnStatus(i,atUpperBound);
3029    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
3030      // free
3031      setColumnStatus(i,isFree);
3032    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
3033      setColumnStatus(i,atLowerBound);
3034    } else {
3035      setColumnStatus(i,atUpperBound);
3036    }
3037#else
3038    setColumnStatus(i,atLowerBound);
3039#endif
3040  }
3041  for (i=0;i<numberRows_;i++) {
3042    setRowStatus(i,basic);
3043  }
3044}
3045/* Loads a problem (the constraints on the
3046   rows are given by lower and upper bounds). If a pointer is 0 then the
3047   following values are the default:
3048   <ul>
3049   <li> <code>colub</code>: all columns have upper bound infinity
3050   <li> <code>collb</code>: all columns have lower bound 0
3051   <li> <code>rowub</code>: all rows have upper bound infinity
3052   <li> <code>rowlb</code>: all rows have lower bound -infinity
3053   <li> <code>obj</code>: all variables have 0 objective coefficient
3054   </ul>
3055*/
3056void 
3057ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
3058                    const double* collb, const double* colub,   
3059                    const double* obj,
3060                    const double* rowlb, const double* rowub,
3061                    const double * rowObjective)
3062{
3063  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3064                        rowObjective);
3065  createStatus();
3066}
3067void 
3068ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
3069                    const double* collb, const double* colub,   
3070                    const double* obj,
3071                    const double* rowlb, const double* rowub,
3072                    const double * rowObjective)
3073{
3074  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3075                        rowObjective);
3076  createStatus();
3077}
3078
3079/* Just like the other loadProblem() method except that the matrix is
3080   given in a standard column major ordered format (without gaps). */
3081void 
3082ClpSimplex::loadProblem (  const int numcols, const int numrows,
3083                    const CoinBigIndex* start, const int* index,
3084                    const double* value,
3085                    const double* collb, const double* colub,   
3086                    const double* obj,
3087                    const double* rowlb, const double* rowub,
3088                    const double * rowObjective)
3089{
3090  ClpModel::loadProblem(numcols, numrows, start, index, value,
3091                          collb, colub, obj, rowlb, rowub,
3092                          rowObjective);
3093  createStatus();
3094}
3095void 
3096ClpSimplex::loadProblem (  const int numcols, const int numrows,
3097                           const CoinBigIndex* start, const int* index,
3098                           const double* value,const int * length,
3099                           const double* collb, const double* colub,   
3100                           const double* obj,
3101                           const double* rowlb, const double* rowub,
3102                           const double * rowObjective)
3103{
3104  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3105                          collb, colub, obj, rowlb, rowub,
3106                          rowObjective);
3107  createStatus();
3108}
3109// Read an mps file from the given filename
3110int 
3111ClpSimplex::readMps(const char *filename,
3112            bool keepNames,
3113            bool ignoreErrors)
3114{
3115  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3116  createStatus();
3117  return status;
3118}
3119// Just check solution (for external use)
3120void 
3121ClpSimplex::checkSolution()
3122{
3123  // put in standard form
3124  createRim(7+8+16);
3125  dualTolerance_=dblParam_[ClpDualTolerance];
3126  primalTolerance_=dblParam_[ClpPrimalTolerance];
3127  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
3128  checkDualSolution();
3129#ifdef CLP_DEBUG
3130  int i;
3131  double value=0.0;
3132  for (i=0;i<numberRows_+numberColumns_;i++)
3133    value += dj_[i]*solution_[i];
3134  printf("dual value %g, primal %g\n",value,objectiveValue_);
3135#endif
3136  // release extra memory
3137  deleteRim(false);
3138}
Note: See TracBrowser for help on using the repository browser.