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

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

out CoinCopy? and CoinFill?

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