source: trunk/ClpSimplex.cpp @ 171

Last change on this file since 171 was 171, checked in by forrest, 18 years ago

Dual values pass fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 121.3 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpPackedMatrix.hpp"
15#include "CoinIndexedVector.hpp"
16#include "ClpDualRowDantzig.hpp"
17#include "ClpDualRowSteepest.hpp"
18#include "ClpPrimalColumnDantzig.hpp"
19#include "ClpPrimalColumnSteepest.hpp"
20#include "ClpNonLinearCost.hpp"
21#include "ClpMessage.hpp"
22#include "ClpLinearObjective.hpp"
23#include <cfloat>
24
25#include <string>
26#include <stdio.h>
27#include <iostream>
28//#############################################################################
29
30ClpSimplex::ClpSimplex () :
31
32  ClpModel(),
33  columnPrimalInfeasibility_(0.0),
34  rowPrimalInfeasibility_(0.0),
35  columnPrimalSequence_(-2),
36  rowPrimalSequence_(-2), 
37  columnDualInfeasibility_(0.0),
38  rowDualInfeasibility_(0.0),
39  columnDualSequence_(-2),
40  rowDualSequence_(-2),
41  primalToleranceToGetOptimal_(-1.0),
42  remainingDualInfeasibility_(0.0),
43  largeValue_(1.0e15),
44  largestPrimalError_(0.0),
45  largestDualError_(0.0),
46  largestSolutionError_(0.0),
47  dualBound_(1.0e10),
48  alpha_(0.0),
49  theta_(0.0),
50  lowerIn_(0.0),
51  valueIn_(0.0),
52  upperIn_(0.0),
53  dualIn_(0.0),
54  lowerOut_(-1),
55  valueOut_(-1),
56  upperOut_(-1),
57  dualOut_(-1),
58  dualTolerance_(0.0),
59  primalTolerance_(0.0),
60  sumDualInfeasibilities_(0.0),
61  sumPrimalInfeasibilities_(0.0),
62  infeasibilityCost_(1.0e10),
63  sumOfRelaxedDualInfeasibilities_(0.0),
64  sumOfRelaxedPrimalInfeasibilities_(0.0),
65  lower_(NULL),
66  rowLowerWork_(NULL),
67  columnLowerWork_(NULL),
68  upper_(NULL),
69  rowUpperWork_(NULL),
70  columnUpperWork_(NULL),
71  cost_(NULL),
72  rowObjectiveWork_(NULL),
73  objectiveWork_(NULL),
74  sequenceIn_(-1),
75  directionIn_(-1),
76  sequenceOut_(-1),
77  directionOut_(-1),
78  pivotRow_(-1),
79  lastGoodIteration_(-100),
80  dj_(NULL),
81  rowReducedCost_(NULL),
82  reducedCostWork_(NULL),
83  solution_(NULL),
84  rowActivityWork_(NULL),
85  columnActivityWork_(NULL),
86  numberDualInfeasibilities_(0),
87  numberDualInfeasibilitiesWithoutFree_(0),
88  numberPrimalInfeasibilities_(100),
89  numberRefinements_(0),
90  pivotVariable_(NULL),
91  factorization_(NULL),
92  rowScale_(NULL),
93  savedSolution_(NULL),
94  columnScale_(NULL),
95  scalingFlag_(1),
96  numberTimesOptimal_(0),
97  changeMade_(1),
98  algorithm_(0),
99  forceFactorization_(-1),
100  perturbation_(100),
101  nonLinearCost_(NULL),
102  specialOptions_(0),
103  lastBadIteration_(-999999),
104  numberFake_(0),
105  progressFlag_(0),
106  firstFree_(-1)
107
108{
109  int i;
110  for (i=0;i<6;i++) {
111    rowArray_[i]=NULL;
112    columnArray_[i]=NULL;
113  }
114  saveStatus_=NULL;
115  // get an empty factorization so we can set tolerances etc
116  factorization_ = new ClpFactorization();
117  // Say sparse
118  factorization_->sparseThreshold(1);
119  // say Steepest pricing
120  dualRowPivot_ = new ClpDualRowSteepest();
121  // say Steepest pricing
122  primalColumnPivot_ = new ClpPrimalColumnSteepest();
123  solveType_=1; // say simplex based life form
124 
125}
126
127
128//-----------------------------------------------------------------------------
129
130ClpSimplex::~ClpSimplex ()
131{
132  gutsOfDelete(0);
133  delete nonLinearCost_;
134}
135//#############################################################################
136void ClpSimplex::setLargeValue( double value) 
137{
138  if (value>0.0&&value<COIN_DBL_MAX)
139    largeValue_=value;
140}
141int
142ClpSimplex::gutsOfSolution ( const double * rowActivities,
143                             const double * columnActivities,
144                             double * givenDuals,
145                             const double * givenPrimals,
146                             bool valuesPass)
147{
148
149
150  // if values pass, save values of basic variables
151  double * save = NULL;
152  double oldValue=0.0;
153  if (valuesPass) {
154    assert(algorithm_>0); // only primal at present
155    assert(nonLinearCost_);
156    int iRow;
157    checkPrimalSolution( rowActivities, columnActivities);
158    // get correct bounds on all variables
159    nonLinearCost_->checkInfeasibilities();
160    oldValue = nonLinearCost_->largestInfeasibility();
161    save = new double[numberRows_];
162    for (iRow=0;iRow<numberRows_;iRow++) {
163      int iPivot=pivotVariable_[iRow];
164      save[iRow] = solution_[iPivot];
165    }
166  }
167  // do work
168  computePrimals(rowActivities, columnActivities);
169  // If necessary - override results
170  if (givenPrimals) {
171    memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double));
172    memset(rowActivityWork_,0,numberRows_*sizeof(double));
173    times(-1.0,columnActivityWork_,rowActivityWork_);
174  }
175  double objectiveModification = 0.0;
176  if (algorithm_>0&&nonLinearCost_!=NULL) {
177    // primal algorithm
178    // get correct bounds on all variables
179    nonLinearCost_->checkInfeasibilities();
180    objectiveModification += nonLinearCost_->changeInCost();
181    if (nonLinearCost_->numberInfeasibilities())
182      handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
183        <<nonLinearCost_->changeInCost()
184        <<nonLinearCost_->numberInfeasibilities()
185        <<CoinMessageEol;
186  }
187  if (valuesPass) {
188#ifdef CLP_DEBUG
189    std::cout<<"Largest given infeasibility "<<oldValue
190             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
191#endif
192    int numberOut=0;
193    if (oldValue<1.0
194        &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
195        largestPrimalError_>1.0e-3) {
196      // throw out up to 1000 structurals
197      int iRow;
198      int * sort = new int[numberRows_];
199      // first put back solution and store difference
200      for (iRow=0;iRow<numberRows_;iRow++) {
201        int iPivot=pivotVariable_[iRow];
202        double difference = fabs(solution_[iPivot]=save[iRow]);
203        solution_[iPivot]=save[iRow];
204        save[iRow]=difference;
205      }
206      for (iRow=0;iRow<numberRows_;iRow++) {
207        int iPivot=pivotVariable_[iRow];
208
209        if (iPivot<numberColumns_) {
210          // column
211          double difference= save[iRow];
212          if (difference>1.0e-4) {
213            sort[numberOut]=iPivot;
214            save[numberOut++]=difference;
215          }
216        }
217      }
218      CoinSort_2(save, save + numberOut, sort,
219                 CoinFirstGreater_2<double, int>());
220      numberOut = min(1000,numberOut);
221      for (iRow=0;iRow<numberOut;iRow++) {
222        int iColumn=sort[iRow];
223        setColumnStatus(iColumn,superBasic);
224
225      }
226      delete [] sort;
227    }
228    delete [] save;
229    if (numberOut)
230      return numberOut;
231  }
232
233  computeDuals(givenDuals);
234
235  // now check solutions
236  checkPrimalSolution( rowActivities, columnActivities);
237  objectiveValue_ += objectiveModification;
238  checkDualSolution();
239  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
240                               largestDualError_>1.0e-2)) 
241    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
242      <<largestPrimalError_
243      <<largestDualError_
244      <<CoinMessageEol;
245  return 0;
246}
247void
248ClpSimplex::computePrimals ( const double * rowActivities,
249                                     const double * columnActivities)
250{
251
252  //work space
253  CoinIndexedVector  * workSpace = rowArray_[0];
254
255  double * array = new double [numberRows_+1]; // +1 for network
256  double * save = new double [numberRows_+1];
257  double * previous = new double [numberRows_+1];
258
259  // accumulate non basic stuff
260  ClpFillN(array,numberRows_,0.0);
261
262  int iRow;
263  // order is this way for scaling
264  // Use whole matrix every time to make it easier for ClpMatrixBase
265  // So zero out basic
266  if (columnActivities!=columnActivityWork_)
267    ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
268  if (rowActivities!=rowActivityWork_)
269    ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
270  for (iRow=0;iRow<numberRows_;iRow++) {
271    int iPivot=pivotVariable_[iRow];
272    solution_[iPivot] = 0.0;
273  }
274  times(-1.0,columnActivityWork_,array);
275
276  for (iRow=0;iRow<numberRows_;iRow++) {
277    array[iRow] += rowActivityWork_[iRow];
278  }
279
280  // Ftran adjusted RHS and iterate to improve accuracy
281  double lastError=COIN_DBL_MAX;
282  int iRefine;
283  double * work = workSpace->denseVector();
284  factorization_->updateColumn(workSpace,array);
285
286  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
287    // save in case we want to restore
288    ClpDisjointCopyN ( array, numberRows_ , save);
289
290    // put solution in correct place
291    for (iRow=0;iRow<numberRows_;iRow++) {
292      int iPivot=pivotVariable_[iRow];
293      solution_[iPivot] = array[iRow];
294    }
295    // check Ax == b  (for all)
296    times(-1.0,columnActivityWork_,work);
297    for (iRow=0;iRow<numberRows_;iRow++) {
298      work[iRow] += rowActivityWork_[iRow];
299    }
300
301    largestPrimalError_=0.0;
302    for (iRow=0;iRow<numberRows_;iRow++) {
303      if (fabs(work[iRow])>largestPrimalError_) {
304        largestPrimalError_=fabs(work[iRow]);
305      }
306      //work[iRow] -= save[iRow];
307    }
308    if (largestPrimalError_>=lastError) {
309      // restore
310      double * temp = array;
311      array = previous;
312      previous=temp;
313      break;
314    }
315    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
316      // try and make better
317      // save this
318      double * temp = array;
319      array = previous;
320      previous=temp;
321      double multiplier = 131072.0;
322      for (iRow=0;iRow<numberRows_;iRow++) {
323        array[iRow] = multiplier*work[iRow];
324        work[iRow]=0.0;
325      }
326      lastError=largestPrimalError_;
327      factorization_->updateColumn(workSpace,array);
328      multiplier = 1.0/multiplier;
329      for (iRow=0;iRow<numberRows_;iRow++) {
330        array[iRow] = previous[iRow] + multiplier*array[iRow];
331        work[iRow]=0.0;
332      }
333    }
334  }
335
336  // solution as accurate as we are going to get
337  ClpFillN(work,numberRows_,0.0);
338  // put solution in correct place
339  for (iRow=0;iRow<numberRows_;iRow++) {
340    int iPivot=pivotVariable_[iRow];
341    solution_[iPivot] = array[iRow];
342  }
343  delete [] previous;
344  delete [] array;
345  delete [] save;
346}
347// now dual side
348void
349ClpSimplex::computeDuals(double * givenDjs)
350{
351  //work space
352  CoinIndexedVector  * workSpace = rowArray_[0];
353
354  double * array = new double [numberRows_+1]; // +1 for network
355  double * save = new double [numberRows_+1];
356  double * previous = new double [numberRows_+1];
357
358  int iRow;
359#ifdef CLP_DEBUG
360  workSpace->checkClear();
361#endif
362  if (!givenDjs) {
363    for (iRow=0;iRow<numberRows_;iRow++) {
364      int iPivot=pivotVariable_[iRow];
365      array[iRow]=cost_[iPivot];
366    }
367  } else {
368    // dual values pass - djs may not be zero
369    for (iRow=0;iRow<numberRows_;iRow++) {
370      int iPivot=pivotVariable_[iRow];
371      // make sure zero if done
372      if (!pivoted(iPivot))
373        givenDjs[iPivot]=0.0;
374      array[iRow]=cost_[iPivot]-givenDjs[iPivot];
375    }
376  }
377  ClpDisjointCopyN ( array, numberRows_ , save);
378
379  // Btran basic costs and get as accurate as possible
380  double lastError=COIN_DBL_MAX;
381  int iRefine;
382  double * work = workSpace->denseVector();
383  factorization_->updateColumnTranspose(workSpace,array);
384  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
385    // check basic reduced costs zero
386    largestDualError_=0.0;
387    // would be faster to do just for basic but this reduces code
388    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
389    transposeTimes(-1.0,array,reducedCostWork_);
390    if (!givenDjs) {
391      for (iRow=0;iRow<numberRows_;iRow++) {
392        int iPivot=pivotVariable_[iRow];
393        if (iPivot>=numberColumns_) {
394          // slack
395          work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
396            + array[iPivot-numberColumns_];
397        } else {
398          // column
399        work[iRow] = reducedCostWork_[iPivot];
400        }
401        if (fabs(work[iRow])>largestDualError_) {
402          largestDualError_=fabs(work[iRow]);
403        }
404      }
405    } else {
406      for (iRow=0;iRow<numberRows_;iRow++) {
407        int iPivot=pivotVariable_[iRow];
408        if (iPivot>=numberColumns_) {
409          // slack
410          work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
411            + array[iPivot-numberColumns_]-givenDjs[iPivot];
412        } else {
413          // column
414          work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
415        }
416        if (fabs(work[iRow])>largestDualError_) {
417          largestDualError_=fabs(work[iRow]);
418          //assert (largestDualError_<1.0e-7);
419          //if (largestDualError_>1.0e-7)
420          //printf("large dual error %g\n",largestDualError_);
421        }
422      }
423    }
424    if (largestDualError_>=lastError) {
425      // restore
426      double * temp = array;
427      array = previous;
428      previous=temp;
429      break;
430    }
431    if (iRefine<numberRefinements_&&largestDualError_>1.0e-20
432        &&!givenDjs) {
433      // try and make better
434      // save this
435      double * temp = array;
436      array = previous;
437      previous=temp;
438      double multiplier = 131072.0;
439      for (iRow=0;iRow<numberRows_;iRow++) {
440        array[iRow] = multiplier*work[iRow];
441        work[iRow]=0.0;
442      }
443      lastError=largestDualError_;
444      factorization_->updateColumnTranspose(workSpace,array);
445      multiplier = 1.0/multiplier;
446      for (iRow=0;iRow<numberRows_;iRow++) {
447        array[iRow] = previous[iRow] + multiplier*array[iRow];
448        work[iRow]=0.0;
449      }
450    } else {
451      break;
452    }
453  }
454  ClpFillN(work,numberRows_,0.0);
455  // now look at dual solution
456  for (iRow=0;iRow<numberRows_;iRow++) {
457    // slack
458    double value = array[iRow];
459    dual_[iRow]=value;
460    value += rowObjectiveWork_[iRow];
461    rowReducedCost_[iRow]=value;
462  }
463  ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
464  transposeTimes(-1.0,dual_,reducedCostWork_);
465  // If necessary - override results
466  if (givenDjs) {
467    // restore accurate duals
468    memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
469  }
470
471  delete [] previous;
472  delete [] array;
473  delete [] save;
474}
475/* Given an existing factorization computes and checks
476   primal and dual solutions.  Uses input arrays for variables at
477   bounds.  Returns feasibility states */
478int ClpSimplex::getSolution ( const double * rowActivities,
479                               const double * columnActivities)
480{
481  if (!factorization_->status()) {
482    // put in standard form
483    createRim(7+8+16+32);
484    // do work
485    gutsOfSolution ( rowActivities, columnActivities,NULL,NULL);
486    // release extra memory
487    deleteRim(0);
488  }
489  return factorization_->status();
490}
491/* Given an existing factorization computes and checks
492   primal and dual solutions.  Uses current problem arrays for
493   bounds.  Returns feasibility states */
494int ClpSimplex::getSolution ( )
495{
496  double * rowActivities = new double[numberRows_];
497  double * columnActivities = new double[numberColumns_];
498  ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
499  ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
500  int status = getSolution( rowActivities, columnActivities);
501  delete [] rowActivities;
502  delete [] columnActivities;
503  return status;
504}
505// Factorizes using current basis.  This is for external use
506// Return codes are as from ClpFactorization
507int ClpSimplex::factorize ()
508{
509  // put in standard form
510  createRim(7+8+16+32,false);
511  // do work
512  int status = internalFactorize(-1);
513  // release extra memory
514  deleteRim(0);
515
516  return status;
517}
518
519/* Factorizes using current basis. 
520      solveType - 1 iterating, 0 initial, -1 external
521      If 10 added then in primal values pass
522*/
523/* Return codes are as from ClpFactorization unless initial factorization
524   when total number of singularities is returned */
525int ClpSimplex::internalFactorize ( int solveType)
526{
527
528  int iRow,iColumn;
529  int totalSlacks=numberRows_;
530
531  bool valuesPass=false;
532  if (solveType>=10) {
533    valuesPass=true;
534    solveType -= 10;
535  }
536#ifdef CLP_DEBUG
537  if (solveType==1) {
538    int numberFreeIn=0,numberFreeOut=0;
539    double biggestDj=0.0;
540    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
541      switch(getColumnStatus(iColumn)) {
542       
543      case basic:
544        if (columnLower_[iColumn]<-largeValue_
545            &&columnUpper_[iColumn]>largeValue_) 
546          numberFreeIn++;
547        break;
548      default:
549        if (columnLower_[iColumn]<-largeValue_
550            &&columnUpper_[iColumn]>largeValue_) {
551          numberFreeOut++;
552          biggestDj = max(fabs(dj_[iColumn]),biggestDj);
553        }
554        break;
555      }
556    }
557    if (numberFreeIn+numberFreeOut)
558      printf("%d in basis, %d out - largest dj %g\n",
559             numberFreeIn,numberFreeOut,biggestDj);
560  }
561#endif
562  if (!solveType) {
563    if (!valuesPass) {
564      // not values pass so set to bounds
565      bool allSlack=true;
566      if (status_) {
567        for (iRow=0;iRow<numberRows_;iRow++) {
568          if (getRowStatus(iRow)!=basic) {
569            allSlack=false;
570            break;
571          }
572        }
573      }
574      if (!allSlack) {
575        // set values from warm start (if sensible)
576        int numberBasic=0;
577        for (iRow=0;iRow<numberRows_;iRow++) {
578          switch(getRowStatus(iRow)) {
579           
580          case basic:
581            numberBasic++;
582            break;
583          case atUpperBound:
584            rowActivityWork_[iRow]=rowUpperWork_[iRow];
585            if (rowActivityWork_[iRow]>largeValue_) {
586              if (rowLowerWork_[iRow]>-largeValue_) {
587                rowActivityWork_[iRow]=rowLowerWork_[iRow];
588                setRowStatus(iRow,atLowerBound);
589              } else {
590                // say free
591                setRowStatus(iRow,isFree);
592                rowActivityWork_[iRow]=0.0;
593              }
594            }
595            break;
596          case ClpSimplex::isFixed:
597          case atLowerBound:
598            rowActivityWork_[iRow]=rowLowerWork_[iRow];
599            if (rowActivityWork_[iRow]<-largeValue_) {
600              if (rowUpperWork_[iRow]<largeValue_) {
601                rowActivityWork_[iRow]=rowUpperWork_[iRow];
602                setRowStatus(iRow,atUpperBound);
603              } else {
604                // say free
605                setRowStatus(iRow,isFree);
606                rowActivityWork_[iRow]=0.0;
607              }
608            }
609            break;
610          case isFree:
611            // may not be free after all
612            if (rowLowerWork_[iRow]<-largeValue_&&rowUpperWork_[iRow]>largeValue_)
613              break;
614            // not really free - fall through to superbasic
615          case superBasic:
616            if (rowUpperWork_[iRow]>largeValue_) {
617              if (rowLowerWork_[iRow]>-largeValue_) {
618                rowActivityWork_[iRow]=rowLowerWork_[iRow];
619                setRowStatus(iRow,atLowerBound);
620              } else {
621                // say free
622                setRowStatus(iRow,isFree);
623                rowActivityWork_[iRow]=0.0;
624              }
625            } else {
626              if (rowLowerWork_[iRow]>-largeValue_) {
627                // set to nearest
628                if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
629                    <fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])) {
630                  rowActivityWork_[iRow]=rowLowerWork_[iRow];
631                  setRowStatus(iRow,atLowerBound);
632                } else {
633                  rowActivityWork_[iRow]=rowUpperWork_[iRow];
634                  setRowStatus(iRow,atUpperBound);
635                }
636              } else {
637                rowActivityWork_[iRow]=rowUpperWork_[iRow];
638                setRowStatus(iRow,atUpperBound);
639              }
640            }
641            break;
642          }
643        }
644        totalSlacks=numberBasic;
645
646        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
647          switch(getColumnStatus(iColumn)) {
648           
649          case basic:
650            if (numberBasic==numberRows_) {
651              // take out of basis
652              if (columnLowerWork_[iColumn]>-largeValue_) {
653                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
654                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
655                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
656                  setColumnStatus(iColumn,atLowerBound);
657                } else {
658                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
659                  setColumnStatus(iColumn,atUpperBound);
660                }
661              } else if (columnUpperWork_[iColumn]<largeValue_) {
662                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
663                setColumnStatus(iColumn,atUpperBound);
664              } else {
665                columnActivityWork_[iColumn]=0.0;
666                setColumnStatus(iColumn,isFree);
667              }
668            } else {
669              numberBasic++;
670            }
671            break;
672          case atUpperBound:
673            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
674            if (columnActivityWork_[iColumn]>largeValue_) {
675              if (columnLowerWork_[iColumn]<-largeValue_) {
676                columnActivityWork_[iColumn]=0.0;
677                setColumnStatus(iColumn,isFree);
678              } else {
679                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
680                setColumnStatus(iColumn,atLowerBound);
681              }
682            }
683            break;
684          case atLowerBound:
685          case ClpSimplex::isFixed:
686            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
687            if (columnActivityWork_[iColumn]<-largeValue_) {
688              if (columnUpperWork_[iColumn]>largeValue_) {
689                columnActivityWork_[iColumn]=0.0;
690                setColumnStatus(iColumn,isFree);
691              } else {
692                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
693                setColumnStatus(iColumn,atUpperBound);
694              }
695            }
696            break;
697          case isFree:
698            // may not be free after all
699            if (columnLowerWork_[iColumn]<-largeValue_&&columnUpperWork_[iColumn]>largeValue_)
700              break;
701            // not really free - fall through to superbasic
702          case superBasic:
703            if (columnUpperWork_[iColumn]>largeValue_) {
704              if (columnLowerWork_[iColumn]>-largeValue_) {
705                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
706                setColumnStatus(iColumn,atLowerBound);
707              } else {
708                // say free
709                setColumnStatus(iColumn,isFree);
710                columnActivityWork_[iColumn]=0.0;
711              }
712            } else {
713              if (columnLowerWork_[iColumn]>-largeValue_) {
714                // set to nearest
715                if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
716                    <fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])) {
717                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
718                  setColumnStatus(iColumn,atLowerBound);
719                } else {
720                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
721                  setColumnStatus(iColumn,atUpperBound);
722                }
723              } else {
724                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
725                setColumnStatus(iColumn,atUpperBound);
726              }
727            }
728            break;
729          }
730        }
731      } else {
732        // all slack basis
733        int numberBasic=0;
734        if (!status_) {
735          createStatus();
736        }
737        for (iRow=0;iRow<numberRows_;iRow++) {
738          double lower=rowLowerWork_[iRow];
739          double upper=rowUpperWork_[iRow];
740          if (lower>-largeValue_||upper<largeValue_) {
741            if (fabs(lower)<=fabs(upper)) {
742              rowActivityWork_[iRow]=lower;
743            } else {
744              rowActivityWork_[iRow]=upper;
745            }
746          } else {
747            rowActivityWork_[iRow]=0.0;
748          }
749          setRowStatus(iRow,basic);
750          numberBasic++;
751        }
752        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
753          double lower=columnLowerWork_[iColumn];
754          double upper=columnUpperWork_[iColumn];
755          double big_bound = largeValue_;
756          if (lower>-big_bound||upper<big_bound) {
757            if ((getColumnStatus(iColumn)==atLowerBound&&
758                 columnActivityWork_[iColumn]==lower)||
759                (getColumnStatus(iColumn)==atUpperBound&&
760                 columnActivityWork_[iColumn]==upper)) {
761              // status looks plausible
762            } else {
763              // set to sensible
764              if (fabs(lower)<=fabs(upper)) {
765                setColumnStatus(iColumn,atLowerBound);
766                columnActivityWork_[iColumn]=lower;
767              } else {
768                setColumnStatus(iColumn,atUpperBound);
769                columnActivityWork_[iColumn]=upper;
770              }
771            }
772          } else {
773            setColumnStatus(iColumn,isFree);
774            columnActivityWork_[iColumn]=0.0;
775          }
776        }
777      }
778    } else {
779      // values pass has less coding
780      // make row activities correct
781      memset(rowActivityWork_,0,numberRows_*sizeof(double));
782      times(1.0,columnActivityWork_,rowActivityWork_);
783      if (status_) {
784        // set values from warm start (if sensible)
785        int numberBasic=0;
786        for (iRow=0;iRow<numberRows_;iRow++) {
787          if (getRowStatus(iRow)==basic) 
788            numberBasic++;
789          else {
790            setRowStatus(iRow,superBasic);
791            // but put to bound if close
792            if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
793                <=primalTolerance_) {
794              rowActivityWork_[iRow]=rowLowerWork_[iRow];
795              setRowStatus(iRow,atLowerBound);
796            } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
797                       <=primalTolerance_) {
798              rowActivityWork_[iRow]=rowUpperWork_[iRow];
799              setRowStatus(iRow,atUpperBound);
800            }
801          }
802         
803        }
804        totalSlacks=numberBasic;
805        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
806          if (getColumnStatus(iColumn)==basic) {
807            if (numberBasic==numberRows_) {
808              // take out of basis
809              setColumnStatus(iColumn,superBasic);
810              // but put to bound if close
811              if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
812                  <=primalTolerance_) {
813                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
814                setColumnStatus(iColumn,atLowerBound);
815              } else if (fabs(columnActivityWork_[iColumn]
816                            -columnUpperWork_[iColumn])
817                         <=primalTolerance_) {
818                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
819                setColumnStatus(iColumn,atUpperBound);
820              }
821            } else 
822              numberBasic++;
823          } else {
824            setColumnStatus(iColumn,superBasic);
825            // but put to bound if close
826            if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
827                <=primalTolerance_) {
828              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
829              setColumnStatus(iColumn,atLowerBound);
830            } else if (fabs(columnActivityWork_[iColumn]
831                          -columnUpperWork_[iColumn])
832                       <=primalTolerance_) {
833              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
834              setColumnStatus(iColumn,atUpperBound);
835            }
836          }
837        }
838      } else {
839        // all slack basis
840        int numberBasic=0;
841        if (!status_) {
842          createStatus();
843        }
844        for (iRow=0;iRow<numberRows_;iRow++) {
845          setRowStatus(iRow,basic);
846          numberBasic++;
847        }
848        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
849          setColumnStatus(iColumn,superBasic);
850          // but put to bound if close
851          if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
852              <=primalTolerance_) {
853            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
854            setColumnStatus(iColumn,atLowerBound);
855          } else if (fabs(columnActivityWork_[iColumn]
856                        -columnUpperWork_[iColumn])
857                     <=primalTolerance_) {
858            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
859            setColumnStatus(iColumn,atUpperBound);
860          }
861        }
862      }
863    }
864    numberRefinements_=1;
865    // set fixed if they are
866    for (iRow=0;iRow<numberRows_;iRow++) {
867      if (getRowStatus(iRow)!=basic ) {
868        if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
869          rowActivityWork_[iRow]=rowLowerWork_[iRow];
870          setRowStatus(iRow,isFixed);
871        }
872      }
873    }
874    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
875      if (getColumnStatus(iColumn)!=basic ) {
876        if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
877          columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
878          setColumnStatus(iColumn,isFixed);
879        }
880      }
881    }
882  }
883  int status=-99;
884  int * rowIsBasic = new int[numberRows_];
885  int * columnIsBasic = new int[numberColumns_];
886  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
887  while (status<-98) {
888 
889    int i;
890    int numberBasic=0;
891    for (i=0;i<numberRows_;i++) {
892      if (getRowStatus(i) == basic) {
893        rowIsBasic[i]=1;
894        numberBasic++;
895      } else {
896        rowIsBasic[i]=-1;
897      }
898    }
899    for (i=0;i<numberColumns_;i++) {
900      if (getColumnStatus(i) == basic) {
901        columnIsBasic[i]=1;
902        numberBasic++;
903      } else {
904        columnIsBasic[i]=-1;
905      }
906    }
907    assert (numberBasic<=numberRows_);
908    while (status==-99) {
909      status =  factorization_->factorize(this,matrix_,
910                                          numberRows_,numberColumns_,
911                                          rowIsBasic, columnIsBasic,
912                                          0.0);
913      if (status==-99) {
914        // get more memory
915        factorization_->areaFactor(2.0*factorization_->areaFactor());
916      }
917    }
918    if (!status) {
919  /// See whether to redo pivot order
920      if (factorization_->needToReorder()||!numberIterations_) {
921        // do pivot information
922        for (i=0;i<numberRows_;i++) {
923          if (getRowStatus(i) == basic) {
924            pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
925          }
926        }
927        for (i=0;i<numberColumns_;i++) {
928          if (getColumnStatus(i) == basic) {
929            pivotVariable_[columnIsBasic[i]]=i;
930          }
931        }
932      }
933    } else {
934      // leave pivotVariable_ in useful form for cleaning basis
935      for (i=0;i<numberRows_;i++) {
936        pivotVariable_[i]=-1;
937      }
938      for (i=0;i<numberRows_;i++) {
939        if (getRowStatus(i) == basic) {
940          int iPivot = rowIsBasic[i];
941          if (iPivot>=0) 
942            pivotVariable_[iPivot]=i+numberColumns_;
943        }
944      }
945      for (i=0;i<numberColumns_;i++) {
946        if (getColumnStatus(i) == basic) {
947          int iPivot = columnIsBasic[i];
948          if (iPivot>=0) 
949            pivotVariable_[iPivot]=i;
950        }
951      }
952    }
953    if (status==-1) {
954      if (!solveType) {
955        //redo basis - first take ALL columns out
956        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
957          if (getColumnStatus(iColumn)==
958              basic) {
959            // take out
960            if (!valuesPass) {
961              double lower=columnLowerWork_[iColumn];
962              double upper=columnUpperWork_[iColumn];
963              double value=columnActivityWork_[iColumn];
964              if (lower>-largeValue_||upper<largeValue_) {
965                if (fabs(value-lower)<fabs(value-upper)) {
966                  setColumnStatus(iColumn,atLowerBound);
967                  columnActivityWork_[iColumn]=lower;
968                } else {
969                  setColumnStatus(iColumn,atUpperBound);
970                  columnActivityWork_[iColumn]=upper;
971                }
972              } else {
973                setColumnStatus(iColumn,isFree);
974              }
975            } else {
976              setColumnStatus(iColumn,superBasic);
977            }
978          }
979        }
980        for (iRow=0;iRow<numberRows_;iRow++) {
981          int iSequence=pivotVariable_[iRow];
982          if (iSequence>=0) {
983            // basic
984            if (iSequence>=numberColumns_) {
985              // slack in - leave
986              //assert (iSequence-numberColumns_==iRow);
987            } else {
988              // put back structural
989              setColumnStatus(iSequence,basic);
990            }
991          } else {
992            // put in slack
993            setRowStatus(iRow,basic);
994          }
995        }
996        // signal repeat
997        status=-99;
998        // set fixed if they are
999        for (iRow=0;iRow<numberRows_;iRow++) {
1000          if (getRowStatus(iRow)!=basic ) {
1001            if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
1002              rowActivityWork_[iRow]=rowLowerWork_[iRow];
1003              setRowStatus(iRow,isFixed);
1004            }
1005          }
1006        }
1007        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1008          if (getColumnStatus(iColumn)!=basic ) {
1009            if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
1010              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
1011              setColumnStatus(iColumn,isFixed);
1012            }
1013          }
1014        }
1015      }
1016    }
1017  } 
1018  delete [] rowIsBasic;
1019  delete [] columnIsBasic;
1020  if (status) {
1021    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
1022      <<status
1023      <<CoinMessageEol;
1024    return -1;
1025  } else if (!solveType) {
1026    // Initial basis - return number of singularities
1027    int numberSlacks=0;
1028    for (iRow=0;iRow<numberRows_;iRow++) {
1029      if (getRowStatus(iRow) == basic)
1030        numberSlacks++;
1031    }
1032    status= max(numberSlacks-totalSlacks,0);
1033  }
1034
1035  // sparse methods
1036  //if (factorization_->sparseThreshold()) {
1037    // get default value
1038    factorization_->sparseThreshold(0);
1039    factorization_->goSparse();
1040    //}
1041 
1042  return status;
1043}
1044/*
1045   This does basis housekeeping and does values for in/out variables.
1046   Can also decide to re-factorize
1047*/
1048int 
1049ClpSimplex::housekeeping(double objectiveChange)
1050{
1051  numberIterations_++;
1052  changeMade_++; // something has happened
1053  // incoming variable
1054
1055  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
1056    <<directionOut_
1057    <<directionIn_<<theta_
1058    <<dualOut_<<dualIn_<<alpha_
1059    <<CoinMessageEol;
1060  if (getStatus(sequenceIn_)==isFree) {
1061    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
1062      <<sequenceIn_
1063      <<CoinMessageEol;
1064  }
1065  // change of incoming
1066  char rowcol[]={'R','C'};
1067  if (pivotRow_>=0)
1068    pivotVariable_[pivotRow_]=sequenceIn();
1069  if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20)
1070    progressFlag_ |= 2; // making real progress
1071  solution_[sequenceIn_]=valueIn_;
1072  if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
1073    progressFlag_ |= 1; // making real progress
1074  if (sequenceIn_!=sequenceOut_) {
1075    //assert( getStatus(sequenceOut_)== basic);
1076    setStatus(sequenceIn_,basic);
1077    if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
1078      if (directionOut_>0) {
1079        // going to lower
1080        setStatus(sequenceOut_,atLowerBound);
1081      } else {
1082        // going to upper
1083        setStatus(sequenceOut_,atUpperBound);
1084      }
1085    } else {
1086      // fixed
1087      setStatus(sequenceOut_,isFixed);
1088    }
1089    solution_[sequenceOut_]=valueOut_;
1090  } else {
1091    // flip from bound to bound
1092    if (directionIn_==-1) {
1093      // as if from upper bound
1094      setStatus(sequenceIn_, atLowerBound);
1095    } else {
1096      // as if from lower bound
1097      setStatus(sequenceIn_, atUpperBound);
1098    }
1099  }
1100  objectiveValue_ += objectiveChange;
1101  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
1102    <<numberIterations_<<objectiveValue()
1103    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
1104    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
1105  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
1106  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
1107  handler_->message()<<CoinMessageEol;
1108  if (hitMaximumIterations())
1109    return 2;
1110  // only time to re-factorize if one before real time
1111  // this is so user won't be surprised that maximumPivots has exact meaning
1112  if (factorization_->pivots()==factorization_->maximumPivots()) {
1113    return 1;
1114  } else {
1115    if (forceFactorization_>0&&
1116        factorization_->pivots()==forceFactorization_) {
1117      // relax
1118      forceFactorization_ = (1+3*forceFactorization_)/2;
1119      if (forceFactorization_>factorization_->maximumPivots())
1120        forceFactorization_ = -1; //off
1121      return 1;
1122    } else {
1123      // carry on iterating
1124      return 0;
1125    }
1126  }
1127}
1128// Copy constructor.
1129ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
1130  ClpModel(rhs),
1131  columnPrimalInfeasibility_(0.0),
1132  rowPrimalInfeasibility_(0.0),
1133  columnPrimalSequence_(-2),
1134  rowPrimalSequence_(-2), 
1135  columnDualInfeasibility_(0.0),
1136  rowDualInfeasibility_(0.0),
1137  columnDualSequence_(-2),
1138  rowDualSequence_(-2),
1139  primalToleranceToGetOptimal_(-1.0),
1140  remainingDualInfeasibility_(0.0),
1141  largeValue_(1.0e15),
1142  largestPrimalError_(0.0),
1143  largestDualError_(0.0),
1144  largestSolutionError_(0.0),
1145  dualBound_(1.0e10),
1146  alpha_(0.0),
1147  theta_(0.0),
1148  lowerIn_(0.0),
1149  valueIn_(0.0),
1150  upperIn_(0.0),
1151  dualIn_(0.0),
1152  lowerOut_(-1),
1153  valueOut_(-1),
1154  upperOut_(-1),
1155  dualOut_(-1),
1156  dualTolerance_(0.0),
1157  primalTolerance_(0.0),
1158  sumDualInfeasibilities_(0.0),
1159  sumPrimalInfeasibilities_(0.0),
1160  infeasibilityCost_(1.0e10),
1161  sumOfRelaxedDualInfeasibilities_(0.0),
1162  sumOfRelaxedPrimalInfeasibilities_(0.0),
1163  lower_(NULL),
1164  rowLowerWork_(NULL),
1165  columnLowerWork_(NULL),
1166  upper_(NULL),
1167  rowUpperWork_(NULL),
1168  columnUpperWork_(NULL),
1169  cost_(NULL),
1170  rowObjectiveWork_(NULL),
1171  objectiveWork_(NULL),
1172  sequenceIn_(-1),
1173  directionIn_(-1),
1174  sequenceOut_(-1),
1175  directionOut_(-1),
1176  pivotRow_(-1),
1177  lastGoodIteration_(-100),
1178  dj_(NULL),
1179  rowReducedCost_(NULL),
1180  reducedCostWork_(NULL),
1181  solution_(NULL),
1182  rowActivityWork_(NULL),
1183  columnActivityWork_(NULL),
1184  numberDualInfeasibilities_(0),
1185  numberDualInfeasibilitiesWithoutFree_(0),
1186  numberPrimalInfeasibilities_(100),
1187  numberRefinements_(0),
1188  pivotVariable_(NULL),
1189  factorization_(NULL),
1190  rowScale_(NULL),
1191  savedSolution_(NULL),
1192  columnScale_(NULL),
1193  scalingFlag_(1),
1194  numberTimesOptimal_(0),
1195  changeMade_(1),
1196  algorithm_(0),
1197  forceFactorization_(-1),
1198  perturbation_(100),
1199  nonLinearCost_(NULL),
1200  specialOptions_(0),
1201  lastBadIteration_(-999999),
1202  numberFake_(0),
1203  progressFlag_(0),
1204  firstFree_(-1)
1205{
1206  int i;
1207  for (i=0;i<6;i++) {
1208    rowArray_[i]=NULL;
1209    columnArray_[i]=NULL;
1210  }
1211  saveStatus_=NULL;
1212  factorization_ = NULL;
1213  dualRowPivot_ = NULL;
1214  primalColumnPivot_ = NULL;
1215  gutsOfDelete(0);
1216  specialOptions_ =0;
1217  delete nonLinearCost_;
1218  nonLinearCost_ = NULL;
1219  gutsOfCopy(rhs);
1220  solveType_=1; // say simplex based life form
1221}
1222// Copy constructor from model
1223ClpSimplex::ClpSimplex(const ClpModel &rhs) :
1224  ClpModel(rhs),
1225  columnPrimalInfeasibility_(0.0),
1226  rowPrimalInfeasibility_(0.0),
1227  columnPrimalSequence_(-2),
1228  rowPrimalSequence_(-2), 
1229  columnDualInfeasibility_(0.0),
1230  rowDualInfeasibility_(0.0),
1231  columnDualSequence_(-2),
1232  rowDualSequence_(-2),
1233  primalToleranceToGetOptimal_(-1.0),
1234  remainingDualInfeasibility_(0.0),
1235  largeValue_(1.0e15),
1236  largestPrimalError_(0.0),
1237  largestDualError_(0.0),
1238  largestSolutionError_(0.0),
1239  dualBound_(1.0e10),
1240  alpha_(0.0),
1241  theta_(0.0),
1242  lowerIn_(0.0),
1243  valueIn_(0.0),
1244  upperIn_(0.0),
1245  dualIn_(0.0),
1246  lowerOut_(-1),
1247  valueOut_(-1),
1248  upperOut_(-1),
1249  dualOut_(-1),
1250  dualTolerance_(0.0),
1251  primalTolerance_(0.0),
1252  sumDualInfeasibilities_(0.0),
1253  sumPrimalInfeasibilities_(0.0),
1254  infeasibilityCost_(1.0e10),
1255  sumOfRelaxedDualInfeasibilities_(0.0),
1256  sumOfRelaxedPrimalInfeasibilities_(0.0),
1257  lower_(NULL),
1258  rowLowerWork_(NULL),
1259  columnLowerWork_(NULL),
1260  upper_(NULL),
1261  rowUpperWork_(NULL),
1262  columnUpperWork_(NULL),
1263  cost_(NULL),
1264  rowObjectiveWork_(NULL),
1265  objectiveWork_(NULL),
1266  sequenceIn_(-1),
1267  directionIn_(-1),
1268  sequenceOut_(-1),
1269  directionOut_(-1),
1270  pivotRow_(-1),
1271  lastGoodIteration_(-100),
1272  dj_(NULL),
1273  rowReducedCost_(NULL),
1274  reducedCostWork_(NULL),
1275  solution_(NULL),
1276  rowActivityWork_(NULL),
1277  columnActivityWork_(NULL),
1278  numberDualInfeasibilities_(0),
1279  numberDualInfeasibilitiesWithoutFree_(0),
1280  numberPrimalInfeasibilities_(100),
1281  numberRefinements_(0),
1282  pivotVariable_(NULL),
1283  factorization_(NULL),
1284  rowScale_(NULL),
1285  savedSolution_(NULL),
1286  columnScale_(NULL),
1287  scalingFlag_(1),
1288  numberTimesOptimal_(0),
1289  changeMade_(1),
1290  algorithm_(0),
1291  forceFactorization_(-1),
1292  perturbation_(100),
1293  nonLinearCost_(NULL),
1294  specialOptions_(0),
1295  lastBadIteration_(-999999),
1296  numberFake_(0),
1297  progressFlag_(0),
1298  firstFree_(-1)
1299{
1300  int i;
1301  for (i=0;i<6;i++) {
1302    rowArray_[i]=NULL;
1303    columnArray_[i]=NULL;
1304  }
1305  saveStatus_=NULL;
1306  // get an empty factorization so we can set tolerances etc
1307  factorization_ = new ClpFactorization();
1308  // say Steepest pricing
1309  dualRowPivot_ = new ClpDualRowSteepest();
1310  // say Dantzig pricing
1311  primalColumnPivot_ = new ClpPrimalColumnDantzig();
1312  solveType_=1; // say simplex based life form
1313 
1314}
1315// Assignment operator. This copies the data
1316ClpSimplex & 
1317ClpSimplex::operator=(const ClpSimplex & rhs)
1318{
1319  if (this != &rhs) {
1320    gutsOfDelete(0);
1321    specialOptions_=0;
1322    delete nonLinearCost_;
1323    nonLinearCost_ = NULL;
1324    ClpModel::operator=(rhs);
1325    gutsOfCopy(rhs);
1326  }
1327  return *this;
1328}
1329void 
1330ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1331{
1332  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
1333  rowLowerWork_ = lower_+numberColumns_;
1334  columnLowerWork_ = lower_;
1335  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
1336  rowUpperWork_ = upper_+numberColumns_;
1337  columnUpperWork_ = upper_;
1338  cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
1339  objectiveWork_ = cost_;
1340  rowObjectiveWork_ = cost_+numberColumns_;
1341  dj_ = ClpCopyOfArray(rhs.dj_,numberRows_+numberColumns_);
1342  if (dj_) {
1343    reducedCostWork_ = dj_;
1344    rowReducedCost_ = dj_+numberColumns_;
1345  }
1346  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
1347  if (solution_) {
1348    columnActivityWork_ = solution_;
1349    rowActivityWork_ = solution_+numberColumns_;
1350  }
1351  if (rhs.pivotVariable_) {
1352    pivotVariable_ = new int[numberRows_];
1353    ClpDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
1354  } else {
1355    pivotVariable_=NULL;
1356  }
1357  if (rhs.factorization_) {
1358    factorization_ = new ClpFactorization(*rhs.factorization_);
1359  } else {
1360    factorization_=NULL;
1361  }
1362  rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_);
1363  savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1364  columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_);
1365  int i;
1366  for (i=0;i<6;i++) {
1367    rowArray_[i]=NULL;
1368    if (rhs.rowArray_[i]) 
1369      rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
1370    columnArray_[i]=NULL;
1371    if (rhs.columnArray_[i]) 
1372      columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
1373  }
1374  if (rhs.saveStatus_) {
1375    saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1376  }
1377  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1378  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1379  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1380  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1381  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1382  columnDualSequence_ = rhs.columnDualSequence_;
1383  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1384  rowDualSequence_ = rhs.rowDualSequence_;
1385  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1386  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1387  largeValue_ = rhs.largeValue_;
1388  largestPrimalError_ = rhs.largestPrimalError_;
1389  largestDualError_ = rhs.largestDualError_;
1390  largestSolutionError_ = rhs.largestSolutionError_;
1391  dualBound_ = rhs.dualBound_;
1392  alpha_ = rhs.alpha_;
1393  theta_ = rhs.theta_;
1394  lowerIn_ = rhs.lowerIn_;
1395  valueIn_ = rhs.valueIn_;
1396  upperIn_ = rhs.upperIn_;
1397  dualIn_ = rhs.dualIn_;
1398  sequenceIn_ = rhs.sequenceIn_;
1399  directionIn_ = rhs.directionIn_;
1400  lowerOut_ = rhs.lowerOut_;
1401  valueOut_ = rhs.valueOut_;
1402  upperOut_ = rhs.upperOut_;
1403  dualOut_ = rhs.dualOut_;
1404  sequenceOut_ = rhs.sequenceOut_;
1405  directionOut_ = rhs.directionOut_;
1406  pivotRow_ = rhs.pivotRow_;
1407  lastGoodIteration_ = rhs.lastGoodIteration_;
1408  numberRefinements_ = rhs.numberRefinements_;
1409  dualTolerance_ = rhs.dualTolerance_;
1410  primalTolerance_ = rhs.primalTolerance_;
1411  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1412  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1413  numberDualInfeasibilitiesWithoutFree_ = 
1414    rhs.numberDualInfeasibilitiesWithoutFree_;
1415  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1416  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1417  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1418  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1419  scalingFlag_ = rhs.scalingFlag_;
1420  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1421  changeMade_ = rhs.changeMade_;
1422  algorithm_ = rhs.algorithm_;
1423  forceFactorization_ = rhs.forceFactorization_;
1424  perturbation_ = rhs.perturbation_;
1425  infeasibilityCost_ = rhs.infeasibilityCost_;
1426  specialOptions_ = rhs.specialOptions_;
1427  lastBadIteration_ = rhs.lastBadIteration_;
1428  numberFake_ = rhs.numberFake_;
1429  progressFlag_ = rhs.progressFlag_;
1430  firstFree_ = rhs.firstFree_;
1431  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
1432  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
1433  if (rhs.nonLinearCost_!=NULL)
1434    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1435  else
1436    nonLinearCost_=NULL;
1437  solveType_=rhs.solveType_;
1438}
1439// type == 0 do everything, most + pivot data, 2 factorization data as well
1440void 
1441ClpSimplex::gutsOfDelete(int type)
1442{
1443  delete [] lower_;
1444  lower_=NULL;
1445  rowLowerWork_=NULL;
1446  columnLowerWork_=NULL;
1447  delete [] upper_;
1448  upper_=NULL;
1449  rowUpperWork_=NULL;
1450  columnUpperWork_=NULL;
1451  delete [] cost_;
1452  cost_=NULL;
1453  objectiveWork_=NULL;
1454  rowObjectiveWork_=NULL;
1455  delete [] dj_;
1456  dj_=NULL;
1457  reducedCostWork_=NULL;
1458  rowReducedCost_=NULL;
1459  delete [] solution_;
1460  solution_=NULL;
1461  rowActivityWork_=NULL;
1462  columnActivityWork_=NULL;
1463  delete [] rowScale_;
1464  rowScale_ = NULL;
1465  delete [] savedSolution_;
1466  savedSolution_ = NULL;
1467  delete [] columnScale_;
1468  columnScale_ = NULL;
1469  if ((specialOptions_&2)==0) {
1470    delete nonLinearCost_;
1471    nonLinearCost_ = NULL;
1472  }
1473  int i;
1474  for (i=0;i<6;i++) {
1475    delete rowArray_[i];
1476    rowArray_[i]=NULL;
1477    delete columnArray_[i];
1478    columnArray_[i]=NULL;
1479  }
1480  delete rowCopy_;
1481  rowCopy_=NULL;
1482  delete [] saveStatus_;
1483  saveStatus_=NULL;
1484  if (!type) {
1485    // delete everything
1486    delete factorization_;
1487    factorization_ = NULL;
1488    delete [] pivotVariable_;
1489    pivotVariable_=NULL;
1490    delete dualRowPivot_;
1491    dualRowPivot_ = NULL;
1492    delete primalColumnPivot_;
1493    primalColumnPivot_ = NULL;
1494  } else {
1495    // delete any size information in methods
1496    if (type>1) {
1497      factorization_->clearArrays();
1498      delete [] pivotVariable_;
1499      pivotVariable_=NULL;
1500    }
1501    dualRowPivot_->clearArrays();
1502    primalColumnPivot_->clearArrays();
1503  }
1504}
1505// This sets largest infeasibility and most infeasible
1506void 
1507ClpSimplex::checkPrimalSolution(const double * rowActivities,
1508                                        const double * columnActivities)
1509{
1510  double * solution;
1511  int iRow,iColumn;
1512
1513  objectiveValue_ = 0.0;
1514  firstFree_ = -1;
1515  // now look at primal solution
1516  columnPrimalInfeasibility_=0.0;
1517  columnPrimalSequence_=-1;
1518  rowPrimalInfeasibility_=0.0;
1519  rowPrimalSequence_=-1;
1520  largestSolutionError_=0.0;
1521  solution = rowActivityWork_;
1522  sumPrimalInfeasibilities_=0.0;
1523  numberPrimalInfeasibilities_=0;
1524  double primalTolerance = primalTolerance_;
1525  double relaxedTolerance=primalTolerance_;
1526  // we can't really trust infeasibilities if there is primal error
1527  double error = min(1.0e-3,largestPrimalError_);
1528  // allow tolerance at least slightly bigger than standard
1529  relaxedTolerance = relaxedTolerance +  error;
1530  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
1531
1532  for (iRow=0;iRow<numberRows_;iRow++) {
1533    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
1534    double infeasibility=0.0;
1535    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1536    if (solution[iRow]>rowUpperWork_[iRow]) {
1537      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1538    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1539      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1540    }
1541    if (infeasibility>primalTolerance) {
1542      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1543      if (infeasibility>relaxedTolerance) 
1544        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1545      numberPrimalInfeasibilities_ ++;
1546    }
1547    if (infeasibility>rowPrimalInfeasibility_) {
1548      rowPrimalInfeasibility_=infeasibility;
1549      rowPrimalSequence_=iRow;
1550    }
1551    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1552    if (infeasibility>largestSolutionError_)
1553      largestSolutionError_=infeasibility;
1554  }
1555  solution = columnActivityWork_;
1556  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1557    //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
1558    double infeasibility=0.0;
1559    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1560    if (solution[iColumn]>columnUpperWork_[iColumn]) {
1561      infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1562    } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1563      infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1564    }
1565    if (infeasibility>columnPrimalInfeasibility_) {
1566      columnPrimalInfeasibility_=infeasibility;
1567      columnPrimalSequence_=iColumn;
1568    }
1569    if (infeasibility>primalTolerance) {
1570      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1571      if (infeasibility>relaxedTolerance) 
1572        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1573      numberPrimalInfeasibilities_ ++;
1574    }
1575    infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1576    if (infeasibility>largestSolutionError_)
1577      largestSolutionError_=infeasibility;
1578  }
1579}
1580void 
1581ClpSimplex::checkDualSolution()
1582{
1583
1584  int iRow,iColumn;
1585  sumDualInfeasibilities_=0.0;
1586  numberDualInfeasibilities_=0;
1587  numberDualInfeasibilitiesWithoutFree_=0;
1588  columnDualInfeasibility_=0.0;
1589  columnDualSequence_=-1;
1590  rowDualInfeasibility_=0.0;
1591  rowDualSequence_=-1;
1592  int firstFree = -1;
1593  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
1594                                   columnPrimalInfeasibility_);
1595  remainingDualInfeasibility_=0.0;
1596  double relaxedTolerance=dualTolerance_;
1597  // we can't really trust infeasibilities if there is dual error
1598  double error = min(1.0e-3,largestDualError_);
1599  // allow tolerance at least slightly bigger than standard
1600  relaxedTolerance = relaxedTolerance +  error;
1601  sumOfRelaxedDualInfeasibilities_ = 0.0;
1602
1603  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1604    if (getColumnStatus(iColumn) != basic) {
1605      // not basic
1606      double distanceUp = columnUpperWork_[iColumn]-
1607        columnActivityWork_[iColumn];
1608      double distanceDown = columnActivityWork_[iColumn] -
1609        columnLowerWork_[iColumn];
1610      if (distanceUp>primalTolerance_) {
1611        double value = reducedCostWork_[iColumn];
1612        if (!flagged(iColumn)) {
1613          // Check if "free"
1614          if (firstFree<0&&distanceDown>primalTolerance_) {
1615            if (algorithm_>0) {
1616              // primal
1617              firstFree = iColumn;
1618            } else if (fabs(value)>1.0e2*relaxedTolerance) {
1619              // dual with dj
1620              firstFree = iColumn;
1621            }
1622          }
1623        }
1624        // should not be negative
1625        if (value<0.0) {
1626          value = - value;
1627          if (value>columnDualInfeasibility_) {
1628            columnDualInfeasibility_=value;
1629            columnDualSequence_=iColumn;
1630          }
1631          if (value>dualTolerance_) {
1632            if (getColumnStatus(iColumn) != isFree) {
1633              numberDualInfeasibilitiesWithoutFree_ ++;
1634              sumDualInfeasibilities_ += value-dualTolerance_;
1635              if (value>relaxedTolerance) 
1636                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1637              numberDualInfeasibilities_ ++;
1638            } else {
1639              // free so relax a lot
1640              value *= 0.01;
1641              if (value>dualTolerance_) {
1642                sumDualInfeasibilities_ += value-dualTolerance_;
1643                if (value>relaxedTolerance) 
1644                  sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1645                numberDualInfeasibilities_ ++;
1646              }
1647            }
1648            // maybe we can make feasible by increasing tolerance
1649            if (distanceUp<largeValue_) {
1650              if (distanceUp>primalToleranceToGetOptimal_)
1651                primalToleranceToGetOptimal_=distanceUp;
1652            } else {
1653              //gap too big for any tolerance
1654              remainingDualInfeasibility_=
1655                max(remainingDualInfeasibility_,value);
1656            }
1657          }
1658        }
1659      }
1660      if (distanceDown>primalTolerance_) {
1661        double value = reducedCostWork_[iColumn];
1662        // should not be positive
1663        if (value>0.0) {
1664          if (value>columnDualInfeasibility_) {
1665            columnDualInfeasibility_=value;
1666            columnDualSequence_=iColumn;
1667          }
1668          if (value>dualTolerance_) {
1669            sumDualInfeasibilities_ += value-dualTolerance_;
1670            if (value>relaxedTolerance) 
1671              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1672            numberDualInfeasibilities_ ++;
1673            if (getColumnStatus(iColumn) != isFree) 
1674              numberDualInfeasibilitiesWithoutFree_ ++;
1675            // maybe we can make feasible by increasing tolerance
1676            if (distanceDown<largeValue_&&
1677                distanceDown>primalToleranceToGetOptimal_)
1678              primalToleranceToGetOptimal_=distanceDown;
1679          }
1680        }
1681      }
1682    }
1683  }
1684  for (iRow=0;iRow<numberRows_;iRow++) {
1685    if (getRowStatus(iRow) != basic) {
1686      // not basic
1687      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
1688      double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
1689      if (distanceUp>primalTolerance_) {
1690        double value = rowReducedCost_[iRow];
1691        if (!flagged(iRow+numberColumns_)) {
1692          // Check if "free"
1693          if (firstFree<0&&distanceDown>primalTolerance_) {
1694            if (algorithm_>0) {
1695              // primal
1696              firstFree = iRow+numberColumns_;
1697            } else if (fabs(value)>1.0e2*relaxedTolerance) {
1698              // dual with dj
1699              firstFree = iRow+numberColumns_;
1700            }
1701          }
1702        }
1703        // should not be negative
1704        if (value<0.0) {
1705          value = - value;
1706          if (value>rowDualInfeasibility_) {
1707            rowDualInfeasibility_=value;
1708            rowDualSequence_=iRow;
1709          }
1710          if (value>dualTolerance_) {
1711            sumDualInfeasibilities_ += value-dualTolerance_;
1712            if (value>relaxedTolerance) 
1713              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1714            numberDualInfeasibilities_ ++;
1715            if (getRowStatus(iRow) != isFree) 
1716              numberDualInfeasibilitiesWithoutFree_ ++;
1717            // maybe we can make feasible by increasing tolerance
1718            if (distanceUp<largeValue_) {
1719              if (distanceUp>primalToleranceToGetOptimal_)
1720                primalToleranceToGetOptimal_=distanceUp;
1721            } else {
1722              //gap too big for any tolerance
1723              remainingDualInfeasibility_=
1724                max(remainingDualInfeasibility_,value);
1725            }
1726          }
1727        }
1728      }
1729      if (distanceDown>primalTolerance_) {
1730        double value = rowReducedCost_[iRow];
1731        // should not be positive
1732        if (value>0.0) {
1733          if (value>rowDualInfeasibility_) {
1734            rowDualInfeasibility_=value;
1735            rowDualSequence_=iRow;
1736          }
1737          if (value>dualTolerance_) {
1738            sumDualInfeasibilities_ += value-dualTolerance_;
1739            if (value>relaxedTolerance) 
1740              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1741            numberDualInfeasibilities_ ++;
1742            if (getRowStatus(iRow) != isFree) 
1743              numberDualInfeasibilitiesWithoutFree_ ++;
1744            // maybe we can make feasible by increasing tolerance
1745            if (distanceDown<largeValue_&&
1746                distanceDown>primalToleranceToGetOptimal_)
1747              primalToleranceToGetOptimal_=distanceDown;
1748          }
1749        }
1750      }
1751    }
1752  }
1753  if (firstFree>=0) {
1754    if (algorithm_<0||!numberIterations_)
1755      firstFree_=firstFree;
1756  }
1757}
1758/*
1759  Unpacks one column of the matrix into indexed array
1760*/
1761void 
1762ClpSimplex::unpack(CoinIndexedVector * rowArray) const
1763{
1764  rowArray->clear();
1765  if (sequenceIn_>=numberColumns_) {
1766    //slack
1767    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
1768  } else {
1769    // column
1770    matrix_->unpack(this,rowArray,sequenceIn_);
1771  }
1772}
1773void 
1774ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence) const
1775{
1776  rowArray->clear();
1777  if (sequence>=numberColumns_) {
1778    //slack
1779    rowArray->insert(sequence-numberColumns_,-1.0);
1780  } else {
1781    // column
1782    matrix_->unpack(this,rowArray,sequence);
1783  }
1784}
1785bool
1786ClpSimplex::createRim(int what,bool makeRowCopy)
1787{
1788  bool goodMatrix=true;
1789  int i;
1790  if ((what&(16+32))!=0) {
1791    // move information to work arrays
1792    if (optimizationDirection_<0.0) {
1793      // reverse all dual signs
1794      for (i=0;i<numberColumns_;i++) 
1795        reducedCost_[i] = -reducedCost_[i];
1796      for (i=0;i<numberRows_;i++) 
1797        dual_[i] = -dual_[i];
1798    }
1799    // row reduced costs
1800    if (!dj_) {
1801      dj_ = new double[numberRows_+numberColumns_];
1802      reducedCostWork_ = dj_;
1803      rowReducedCost_ = dj_+numberColumns_;
1804      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1805      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
1806    }
1807    if (!solution_||(what&32)!=0) {
1808      if (!solution_)
1809        solution_ = new double[numberRows_+numberColumns_];
1810      columnActivityWork_ = solution_;
1811      rowActivityWork_ = solution_+numberColumns_;
1812      memcpy(columnActivityWork_,columnActivity_,
1813             numberColumns_*sizeof(double));
1814      memcpy(rowActivityWork_,rowActivity_,
1815             numberRows_*sizeof(double));
1816    }
1817  }
1818  if ((what&16)!=0) {
1819    //check matrix
1820    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20)) {
1821      problemStatus_=4;
1822      goodMatrix= false;
1823    }
1824    if (makeRowCopy) {
1825      delete rowCopy_;
1826      // may return NULL if can't give row copy
1827      rowCopy_ = matrix_->reverseOrderedCopy();
1828#ifdef TAKEOUT
1829      {
1830
1831        ClpPackedMatrix* rowCopy =
1832          dynamic_cast< ClpPackedMatrix*>(rowCopy_);
1833        const int * column = rowCopy->getIndices();
1834        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1835        const double * element = rowCopy->getElements();
1836        int i;
1837        for (i=133;i<numberRows_;i++) {
1838          if (rowStart[i+1]-rowStart[i]==10||rowStart[i+1]-rowStart[i]==15)
1839            printf("Row %d has %d elements\n",i,rowStart[i+1]-rowStart[i]);
1840        }
1841      } 
1842#endif
1843    }
1844  }
1845  if ((what&4)!=0) {
1846    delete [] cost_;
1847    cost_ = new double[numberColumns_+numberRows_];
1848    objectiveWork_ = cost_;
1849    rowObjectiveWork_ = cost_+numberColumns_;
1850    memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double));
1851    if (rowObjective_)
1852      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1853    else
1854      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1855  }
1856  if ((what&1)!=0) {
1857    delete [] lower_;
1858    delete [] upper_;
1859    lower_ = new double[numberColumns_+numberRows_];
1860    upper_ = new double[numberColumns_+numberRows_];
1861    rowLowerWork_ = lower_+numberColumns_;
1862    columnLowerWork_ = lower_;
1863    rowUpperWork_ = upper_+numberColumns_;
1864    columnUpperWork_ = upper_;
1865    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1866    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1867    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1868    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1869    // clean up any mismatches on infinity
1870    for (i=0;i<numberColumns_;i++) {
1871      if (columnLowerWork_[i]<-1.0e30)
1872        columnLowerWork_[i] = -COIN_DBL_MAX;
1873      if (columnUpperWork_[i]>1.0e30)
1874        columnUpperWork_[i] = COIN_DBL_MAX;
1875    }
1876    // clean up any mismatches on infinity
1877    for (i=0;i<numberRows_;i++) {
1878      if (rowLowerWork_[i]<-1.0e30)
1879        rowLowerWork_[i] = -COIN_DBL_MAX;
1880      if (rowUpperWork_[i]>1.0e30)
1881        rowUpperWork_[i] = COIN_DBL_MAX;
1882    }
1883  }
1884  // do scaling if needed
1885  if (scalingFlag_>0&&!rowScale_) {
1886    if (matrix_->scale(this))
1887      scalingFlag_=-scalingFlag_; // not scaled after all
1888  }
1889  if ((what&4)!=0) {
1890    double direction = optimizationDirection_;
1891    // direction is actually scale out not scale in
1892    if (direction)
1893      direction = 1.0/direction;
1894    // but also scale by scale factors
1895    // not really sure about row scaling
1896    if (!rowScale_) {
1897      if (direction!=1.0) {
1898        for (i=0;i<numberRows_;i++)
1899          rowObjectiveWork_[i] *= direction;
1900        for (i=0;i<numberColumns_;i++)
1901          objectiveWork_[i] *= direction;
1902      }
1903    } else {
1904      for (i=0;i<numberRows_;i++)
1905        rowObjectiveWork_[i] *= direction/rowScale_[i];
1906      for (i=0;i<numberColumns_;i++)
1907        objectiveWork_[i] *= direction*columnScale_[i];
1908    }
1909  }
1910  if ((what&(1+32))!=0&&rowScale_) {
1911    for (i=0;i<numberColumns_;i++) {
1912      double multiplier = 1.0/columnScale_[i];
1913      if (columnLowerWork_[i]>-1.0e50)
1914        columnLowerWork_[i] *= multiplier;
1915      if (columnUpperWork_[i]<1.0e50)
1916        columnUpperWork_[i] *= multiplier;
1917    }
1918    for (i=0;i<numberRows_;i++) {
1919      if (rowLowerWork_[i]>-1.0e50)
1920        rowLowerWork_[i] *= rowScale_[i];
1921      if (rowUpperWork_[i]<1.0e50)
1922        rowUpperWork_[i] *= rowScale_[i];
1923    }
1924  }
1925  if ((what&(8+32))!=0&&rowScale_) {
1926    // on entry
1927    for (i=0;i<numberColumns_;i++) {
1928      columnActivityWork_[i] /= columnScale_[i];
1929      reducedCostWork_[i] *= columnScale_[i];
1930    }
1931    for (i=0;i<numberRows_;i++) {
1932      rowActivityWork_[i] *= rowScale_[i];
1933      dual_[i] /= rowScale_[i];
1934      rowReducedCost_[i] = dual_[i];
1935    }
1936  }
1937 
1938  if ((what&16)!=0) {
1939    // check rim of problem okay
1940    if (!sanityCheck())
1941      goodMatrix=false;
1942  } 
1943  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1944  // maybe we need to move scales to SimplexModel for factorization?
1945  if ((what&8)!=0&&!pivotVariable_) {
1946    pivotVariable_=new int[numberRows_];
1947  }
1948
1949  if ((what&16)!=0&&!rowArray_[2]) {
1950    // get some arrays
1951    int iRow,iColumn;
1952    // these are "indexed" arrays so we always know where nonzeros are
1953    /**********************************************************
1954      rowArray_[3] is long enough for rows+columns
1955    *********************************************************/
1956    for (iRow=0;iRow<4;iRow++) {
1957      if (!rowArray_[iRow]) {
1958        rowArray_[iRow]=new CoinIndexedVector();
1959        int length =numberRows_+factorization_->maximumPivots();
1960        if (iRow==3)
1961          length += numberColumns_;
1962        rowArray_[iRow]->reserve(length);
1963      }
1964    }
1965   
1966    for (iColumn=0;iColumn<2;iColumn++) {
1967      if (!columnArray_[iColumn]) {
1968        columnArray_[iColumn]=new CoinIndexedVector();
1969        columnArray_[iColumn]->reserve(numberColumns_);
1970      }
1971    }   
1972  }
1973  return goodMatrix;
1974}
1975void
1976ClpSimplex::deleteRim(int getRidOfFactorizationData)
1977{
1978  int i;
1979  if (problemStatus_!=1&&problemStatus_!=2) {
1980    delete [] ray_;
1981    ray_=NULL;
1982  }
1983  // ray may be null if in branch and bound
1984  if (rowScale_) {
1985    for (i=0;i<numberColumns_;i++) {
1986      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1987      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1988    }
1989    for (i=0;i<numberRows_;i++) {
1990      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1991      dual_[i] *= rowScale_[i];
1992    }
1993    if (problemStatus_==2) {
1994      for (i=0;i<numberColumns_;i++) {
1995        ray_[i] *= columnScale_[i];
1996      }
1997    } else if (problemStatus_==1&&ray_) {
1998      for (i=0;i<numberRows_;i++) {
1999        ray_[i] *= rowScale_[i];
2000      }
2001    }
2002  } else {
2003    for (i=0;i<numberColumns_;i++) {
2004      columnActivity_[i] = columnActivityWork_[i];
2005      reducedCost_[i] = reducedCostWork_[i];
2006    }
2007    for (i=0;i<numberRows_;i++) {
2008      rowActivity_[i] = rowActivityWork_[i];
2009    }
2010  }
2011  // direction may have been modified by scaling - clean up
2012  if (optimizationDirection_>0.0) {
2013    optimizationDirection_ = 1.0;
2014  }  else if (optimizationDirection_<0.0) {
2015    optimizationDirection_ = -1.0;
2016    // and reverse all dual signs
2017    for (i=0;i<numberColumns_;i++) 
2018      reducedCost_[i] = -reducedCost_[i];
2019    for (i=0;i<numberRows_;i++) 
2020      dual_[i] = -dual_[i];
2021  }
2022  // scaling may have been turned off
2023  scalingFlag_ = abs(scalingFlag_);
2024  if(getRidOfFactorizationData>=0)
2025    gutsOfDelete(getRidOfFactorizationData+1);
2026}
2027void 
2028ClpSimplex::setDualBound(double value)
2029{
2030  if (value>0.0)
2031    dualBound_=value;
2032}
2033void 
2034ClpSimplex::setInfeasibilityCost(double value)
2035{
2036  if (value>0.0)
2037    infeasibilityCost_=value;
2038}
2039void ClpSimplex::setnumberRefinements( int value) 
2040{
2041  if (value>=0&&value<10)
2042    numberRefinements_=value;
2043}
2044// Sets row pivot choice algorithm in dual
2045void 
2046ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
2047{
2048  delete dualRowPivot_;
2049  dualRowPivot_ = choice.clone(true);
2050}
2051// Sets row pivot choice algorithm in dual
2052void 
2053ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
2054{
2055  delete primalColumnPivot_;
2056  primalColumnPivot_ = choice.clone(true);
2057}
2058// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
2059void 
2060ClpSimplex::scaling(int mode)
2061{
2062  if (mode>0&&mode<3) {
2063    scalingFlag_=mode;
2064  } else if (!mode) {
2065    scalingFlag_=0;
2066    delete [] rowScale_;
2067    rowScale_ = NULL;
2068    delete [] columnScale_;
2069    columnScale_ = NULL;
2070  }
2071}
2072// Passes in factorization
2073void 
2074ClpSimplex::setFactorization( ClpFactorization & factorization)
2075{
2076  delete factorization_;
2077  factorization_= new ClpFactorization(factorization);
2078}
2079void 
2080ClpSimplex::times(double scalar,
2081                  const double * x, double * y) const
2082{
2083  if (rowScale_)
2084    matrix_->times(scalar,x,y,rowScale_,columnScale_);
2085  else
2086    matrix_->times(scalar,x,y);
2087}
2088void 
2089ClpSimplex::transposeTimes(double scalar,
2090                           const double * x, double * y) const 
2091{
2092  if (rowScale_)
2093    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
2094  else
2095    matrix_->transposeTimes(scalar,x,y);
2096}
2097/* Perturbation:
2098   -50 to +50 - perturb by this power of ten (-6 sounds good)
2099   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
2100   101 - we are perturbed
2101   102 - don't try perturbing again
2102   default is 100
2103*/
2104void 
2105ClpSimplex::setPerturbation(int value)
2106{
2107  if(value<=100&&value >=-50) {
2108    perturbation_=value;
2109  } 
2110}
2111// Sparsity on or off
2112bool 
2113ClpSimplex::sparseFactorization() const
2114{
2115  return factorization_->sparseThreshold()!=0;
2116}
2117void 
2118ClpSimplex::setSparseFactorization(bool value)
2119{
2120  if (value) {
2121    if (!factorization_->sparseThreshold())
2122      factorization_->goSparse();
2123  } else {
2124    factorization_->sparseThreshold(0);
2125  }
2126}
2127/* Tightens primal bounds to make dual faster.  Unless
2128   fixed, bounds are slightly looser than they could be.
2129   This is to make dual go faster and is probably not needed
2130   with a presolve.  Returns non-zero if problem infeasible
2131
2132   Fudge for branch and bound - put bounds on columns of factor *
2133   largest value (at continuous) - should improve stability
2134   in branch and bound on infeasible branches (0.0 is off)
2135*/
2136int 
2137ClpSimplex::tightenPrimalBounds(double factor)
2138{
2139 
2140  // Get a row copy in standard format
2141  CoinPackedMatrix copy;
2142  copy.reverseOrderedCopyOf(*matrix());
2143  // get matrix data pointers
2144  const int * column = copy.getIndices();
2145  const CoinBigIndex * rowStart = copy.getVectorStarts();
2146  const int * rowLength = copy.getVectorLengths(); 
2147  const double * element = copy.getElements();
2148  int numberChanged=1,iPass=0;
2149  double large = largeValue(); // treat bounds > this as infinite
2150  int numberInfeasible=0;
2151  int totalTightened = 0;
2152
2153  double tolerance = primalTolerance();
2154
2155
2156  // Save column bounds
2157  double * saveLower = new double [numberColumns_];
2158  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
2159  double * saveUpper = new double [numberColumns_];
2160  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
2161
2162  int iRow, iColumn;
2163
2164  // If wanted - tighten column bounds using solution
2165  if (factor) {
2166    assert (factor>1.0);
2167    double largest=0.0;
2168    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2169      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
2170        largest = max(largest,fabs(columnActivity_[iColumn]));
2171      }
2172    }
2173    largest *= factor;
2174    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2175      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
2176        columnUpper_[iColumn] = min(columnUpper_[iColumn],largest);
2177        columnLower_[iColumn] = max(columnLower_[iColumn],-largest);
2178      }
2179    }
2180  }
2181#define MAXPASS 10 
2182
2183  // Loop round seeing if we can tighten bounds
2184  // Would be faster to have a stack of possible rows
2185  // and we put altered rows back on stack
2186
2187  while(numberChanged) {
2188
2189    numberChanged = 0; // Bounds tightened this pass
2190   
2191    if (iPass==MAXPASS) break;
2192    iPass++;
2193   
2194    for (iRow = 0; iRow < numberRows_; iRow++) {
2195
2196      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
2197
2198        // possible row
2199        int infiniteUpper = 0;
2200        int infiniteLower = 0;
2201        double maximumUp = 0.0;
2202        double maximumDown = 0.0;
2203        double newBound;
2204        CoinBigIndex rStart = rowStart[iRow];
2205        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2206        CoinBigIndex j;
2207
2208        // Compute possible lower and upper ranges
2209
2210        for (j = rStart; j < rEnd; ++j) {
2211          double value=element[j];
2212          iColumn = column[j];
2213          if (value > 0.0) {
2214            if (columnUpper_[iColumn] >= large) {
2215              maximumUp = 1e31;
2216              ++infiniteUpper;
2217            } else {
2218              maximumUp += columnUpper_[iColumn] * value;
2219            }
2220            if (columnLower_[iColumn] <= -large) {
2221              maximumDown = -1e31;
2222              ++infiniteLower;
2223            } else {
2224              maximumDown += columnLower_[iColumn] * value;
2225            }
2226          } else if (value<0.0) {
2227            if (columnUpper_[iColumn] >= large) {
2228              maximumDown = -1e31;
2229              ++infiniteLower;
2230            } else {
2231              maximumDown += columnUpper_[iColumn] * value;
2232            }
2233            if (columnLower_[iColumn] <= -large) {
2234              maximumUp = 1e31;
2235              ++infiniteUpper;
2236            } else {
2237              maximumUp += columnLower_[iColumn] * value;
2238            }
2239          }
2240        }
2241        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2242            maximumDown >= rowLower_[iRow] - tolerance) {
2243
2244          // Row is redundant - make totally free
2245          rowLower_[iRow]=-COIN_DBL_MAX;
2246          rowUpper_[iRow]=COIN_DBL_MAX;
2247
2248        } else {
2249          if (maximumUp < rowLower_[iRow] -tolerance ||
2250              maximumDown > rowUpper_[iRow]+tolerance) {
2251            // problem is infeasible - exit at once
2252            numberInfeasible++;
2253            if (handler_->logLevel()>3)
2254              printf("a %g %g %g %g\n",maximumUp,rowLower_[iRow],
2255                     maximumDown,rowUpper_[iRow]);
2256            break;
2257          }
2258
2259          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2260            for (j = rStart; j < rEnd; ++j) {
2261              double value=element[j];
2262              iColumn = column[j];
2263              if (value > 0.0) {
2264                if (columnUpper_[iColumn] < large) {
2265                  newBound = columnUpper_[iColumn] + 
2266                    (rowLower_[iRow] - maximumUp) / value;
2267                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2268                    // Tighten the lower bound
2269
2270                    columnLower_[iColumn] = newBound;
2271                    ++numberChanged;
2272
2273                    // check infeasible (relaxed)
2274                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2275                        -100.0*tolerance) {
2276                      numberInfeasible++;
2277                      if (handler_->logLevel()>3)
2278                        printf("b %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
2279                    }
2280                    infiniteLower=1; // skip looking at other bound
2281                  }
2282                }
2283              } else {
2284                if (columnLower_[iColumn] > -large) {
2285                  newBound = columnLower_[iColumn] + 
2286                    (rowLower_[iRow] - maximumUp) / value;
2287                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2288                    // Tighten the upper bound
2289
2290                    columnUpper_[iColumn] = newBound;
2291                    ++numberChanged;
2292
2293                    // check infeasible (relaxed)
2294                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2295                        -100.0*tolerance) {
2296                      numberInfeasible++;
2297                      if (handler_->logLevel()>3)
2298                        printf("c %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
2299                    }
2300                    infiniteLower=1; // skip looking at other bound
2301                  }
2302                }
2303              }
2304            }
2305          }
2306         
2307          // Try other way
2308          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2309            for (j = rStart; j < rEnd; ++j) {
2310              double value=element[j];
2311              iColumn = column[j];
2312              if (value < 0.0) {
2313                if (columnUpper_[iColumn] < large) {
2314                  newBound = columnUpper_[iColumn] + 
2315                    (rowUpper_[iRow] - maximumDown) / value;
2316                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2317                    // Tighten the lower bound
2318
2319                    columnLower_[iColumn] = newBound;
2320                    ++numberChanged;
2321
2322                    // check infeasible (relaxed)
2323                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2324                        -100.0*tolerance) {
2325                      numberInfeasible++;
2326                      if (handler_->logLevel()>3)
2327                        printf("d %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
2328                    }
2329                  }
2330                } 
2331              } else {
2332                if (columnLower_[iColumn] > -large) {
2333                  newBound = columnLower_[iColumn] + 
2334                    (rowUpper_[iRow] - maximumDown) / value;
2335                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2336                    // Tighten the upper bound
2337
2338                    columnUpper_[iColumn] = newBound;
2339                    ++numberChanged;
2340
2341                    // check infeasible (relaxed)
2342                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2343                        -100.0*tolerance) {
2344                      numberInfeasible++;
2345                      if (handler_->logLevel()>3)
2346                        printf("e %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
2347                    }
2348                  }
2349                }
2350              }
2351            }
2352          }
2353        }
2354      }
2355    }
2356    totalTightened += numberChanged;
2357    if (numberInfeasible) break;
2358  }
2359  if (!numberInfeasible) {
2360    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2361      <<totalTightened
2362      <<CoinMessageEol;
2363    // Set bounds slightly loose
2364    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2365      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2366        // fix
2367        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2368          columnLower_[iColumn]=columnUpper_[iColumn];
2369        else
2370          columnUpper_[iColumn]=columnLower_[iColumn];
2371      } else {
2372        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2373          // relax a bit
2374          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2375                                      saveUpper[iColumn]);
2376        }
2377        if (columnLower_[iColumn]>saveLower[iColumn]) {
2378          // relax a bit
2379          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2380                                      saveLower[iColumn]);
2381        }
2382      }
2383    }
2384  } else {
2385    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2386      <<numberInfeasible
2387      <<CoinMessageEol;
2388    // restore column bounds
2389    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2390    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2391  }
2392  delete [] saveLower;
2393  delete [] saveUpper;
2394  return (numberInfeasible);
2395}
2396// dual
2397#include "ClpSimplexDual.hpp"
2398int ClpSimplex::dual (int ifValuesPass )
2399{
2400  return ((ClpSimplexDual *) this)->dual(ifValuesPass);
2401}
2402// primal
2403#include "ClpSimplexPrimal.hpp"
2404int ClpSimplex::primal (int ifValuesPass )
2405{
2406  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2407}
2408#ifdef QUADRATIC
2409#include "ClpSimplexPrimalQuadratic.hpp"
2410/* Solves quadratic problem using SLP - may be used as crash
2411   for other algorithms when number of iterations small
2412*/
2413int 
2414ClpSimplex::quadraticSLP(int numberPasses, double deltaTolerance)
2415{
2416  return ((ClpSimplexPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
2417}
2418// Solves quadratic using Wolfe's algorithm - primal
2419int 
2420ClpSimplex::quadraticWolfe()
2421{
2422  return ((ClpSimplexPrimalQuadratic *) this)->primalWolfe();
2423}
2424#endif
2425/* For strong branching.  On input lower and upper are new bounds
2426   while on output they are objective function values (>1.0e50 infeasible).
2427   Return code is 0 if nothing interesting, -1 if infeasible both
2428   ways and +1 if infeasible one way (check values to see which one(s))
2429*/
2430int ClpSimplex::strongBranching(int numberVariables,const int * variables,
2431                                double * newLower, double * newUpper,
2432                                bool stopOnFirstInfeasible,
2433                                bool alwaysFinish)
2434{
2435  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
2436                                                    newLower,  newUpper,
2437                                                    stopOnFirstInfeasible,
2438                                                    alwaysFinish);
2439}
2440/* Borrow model.  This is so we dont have to copy large amounts
2441   of data around.  It assumes a derived class wants to overwrite
2442   an empty model with a real one - while it does an algorithm.
2443   This is same as ClpModel one, but sets scaling on etc. */
2444void 
2445ClpSimplex::borrowModel(ClpModel & otherModel) 
2446{
2447  ClpModel::borrowModel(otherModel);
2448  createStatus();
2449  ClpDualRowSteepest steep1;
2450  setDualRowPivotAlgorithm(steep1);
2451  ClpPrimalColumnSteepest steep2;
2452  setPrimalColumnPivotAlgorithm(steep2);
2453}
2454typedef struct {
2455  double optimizationDirection;
2456  double dblParam[ClpLastDblParam];
2457  double objectiveValue;
2458  double dualBound;
2459  double dualTolerance;
2460  double primalTolerance;
2461  double sumDualInfeasibilities;
2462  double sumPrimalInfeasibilities;
2463  double infeasibilityCost;
2464  int numberRows;
2465  int numberColumns;
2466  int intParam[ClpLastIntParam];
2467  int numberIterations;
2468  int problemStatus;
2469  int maximumIterations;
2470  int lengthNames;
2471  int numberDualInfeasibilities;
2472  int numberDualInfeasibilitiesWithoutFree;
2473  int numberPrimalInfeasibilities;
2474  int numberRefinements;
2475  int scalingFlag;
2476  int algorithm;
2477  unsigned int specialOptions;
2478  int dualPivotChoice;
2479  int primalPivotChoice;
2480  int matrixStorageChoice;
2481} Clp_scalars;
2482
2483int outDoubleArray(double * array, int length, FILE * fp)
2484{
2485  int numberWritten;
2486  if (array&&length) {
2487    numberWritten = fwrite(&length,sizeof(int),1,fp);
2488    if (numberWritten!=1)
2489      return 1;
2490    numberWritten = fwrite(array,sizeof(double),length,fp);
2491    if (numberWritten!=length)
2492      return 1;
2493  } else {
2494    length = 0;
2495    numberWritten = fwrite(&length,sizeof(int),1,fp);
2496    if (numberWritten!=1)
2497      return 1;
2498  }
2499  return 0;
2500}
2501// Save model to file, returns 0 if success
2502int
2503ClpSimplex::saveModel(const char * fileName)
2504{
2505  FILE * fp = fopen(fileName,"wb");
2506  if (fp) {
2507    Clp_scalars scalars;
2508    int i;
2509    CoinBigIndex numberWritten;
2510    // Fill in scalars
2511    scalars.optimizationDirection = optimizationDirection_;
2512    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2513    scalars.objectiveValue = objectiveValue_;
2514    scalars.dualBound = dualBound_;
2515    scalars.dualTolerance = dualTolerance_;
2516    scalars.primalTolerance = primalTolerance_;
2517    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2518    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2519    scalars.infeasibilityCost = infeasibilityCost_;
2520    scalars.numberRows = numberRows_;
2521    scalars.numberColumns = numberColumns_;
2522    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2523    scalars.numberIterations = numberIterations_;
2524    scalars.problemStatus = problemStatus_;
2525    scalars.maximumIterations = maximumIterations();
2526    scalars.lengthNames = lengthNames_;
2527    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2528    scalars.numberDualInfeasibilitiesWithoutFree
2529      = numberDualInfeasibilitiesWithoutFree_;
2530    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2531    scalars.numberRefinements = numberRefinements_;
2532    scalars.scalingFlag = scalingFlag_;
2533    scalars.algorithm = algorithm_;
2534    scalars.specialOptions = specialOptions_;
2535    scalars.dualPivotChoice = dualRowPivot_->type();
2536    scalars.primalPivotChoice = primalColumnPivot_->type();
2537    scalars.matrixStorageChoice = matrix_->type();
2538
2539    // put out scalars
2540    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2541    if (numberWritten!=1)
2542      return 1;
2543    // strings
2544    CoinBigIndex length;
2545    for (i=0;i<ClpLastStrParam;i++) {
2546      length = strParam_[i].size();
2547      numberWritten = fwrite(&length,sizeof(int),1,fp);
2548      if (numberWritten!=1)
2549        return 1;
2550      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2551      if (numberWritten!=1)
2552        return 1;
2553    }
2554    // arrays - in no particular order
2555    if (outDoubleArray(rowActivity_,numberRows_,fp))
2556        return 1;
2557    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2558        return 1;
2559    if (outDoubleArray(dual_,numberRows_,fp))
2560        return 1;
2561    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2562        return 1;
2563    if (outDoubleArray(rowLower_,numberRows_,fp))
2564        return 1;
2565    if (outDoubleArray(rowUpper_,numberRows_,fp))
2566        return 1;
2567    if (outDoubleArray(objective(),numberColumns_,fp))
2568        return 1;
2569    if (outDoubleArray(rowObjective_,numberRows_,fp))
2570        return 1;
2571    if (outDoubleArray(columnLower_,numberColumns_,fp))
2572        return 1;
2573    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2574        return 1;
2575    if (ray_) {
2576      if (problemStatus_==1)
2577        if (outDoubleArray(ray_,numberRows_,fp))
2578          return 1;
2579      else if (problemStatus_==2)
2580        if (outDoubleArray(ray_,numberColumns_,fp))
2581          return 1;
2582      else
2583        if (outDoubleArray(NULL,0,fp))
2584          return 1;
2585    } else {
2586      if (outDoubleArray(NULL,0,fp))
2587        return 1;
2588    }
2589    if (status_&&(numberRows_+numberColumns_)>0) {
2590      length = numberRows_+numberColumns_;
2591      numberWritten = fwrite(&length,sizeof(int),1,fp);
2592      if (numberWritten!=1)
2593        return 1;
2594      numberWritten = fwrite(status_,sizeof(char),length, fp);
2595      if (numberWritten!=length)
2596        return 1;
2597    } else {
2598      length = 0;
2599      numberWritten = fwrite(&length,sizeof(int),1,fp);
2600      if (numberWritten!=1)
2601        return 1;
2602    }
2603    if (lengthNames_) {
2604      assert (numberRows_ == (int) rowNames_.size());
2605      for (i=0;i<numberRows_;i++) {
2606        length = rowNames_[i].size();
2607        numberWritten = fwrite(&length,sizeof(int),1,fp);
2608        if (numberWritten!=1)
2609          return 1;
2610        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
2611        if (numberWritten!=1)
2612          return 1;
2613      }
2614      assert (numberColumns_ == (int) columnNames_.size());
2615      for (i=0;i<numberColumns_;i++) {
2616        length = columnNames_[i].size();
2617        numberWritten = fwrite(&length,sizeof(int),1,fp);
2618        if (numberWritten!=1)
2619          return 1;
2620        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
2621        if (numberWritten!=1)
2622          return 1;
2623      }
2624    }
2625    // just standard type at present
2626    assert (matrix_->type()==1);
2627    assert (matrix_->getNumCols() == numberColumns_);
2628    assert (matrix_->getNumRows() == numberRows_);
2629    // we are going to save with gaps
2630    length = matrix_->getVectorStarts()[numberColumns_-1]
2631      + matrix_->getVectorLengths()[numberColumns_-1];
2632    numberWritten = fwrite(&length,sizeof(int),1,fp);
2633    if (numberWritten!=1)
2634      return 1;
2635    numberWritten = fwrite(matrix_->getElements(),
2636                           sizeof(double),length,fp);
2637    if (numberWritten!=length)
2638      return 1;
2639    numberWritten = fwrite(matrix_->getIndices(),
2640                           sizeof(int),length,fp);
2641    if (numberWritten!=length)
2642      return 1;
2643    numberWritten = fwrite(matrix_->getVectorStarts(),
2644                           sizeof(int),numberColumns_,fp);
2645    if (numberWritten!=numberColumns_)
2646      return 1;
2647    numberWritten = fwrite(matrix_->getVectorLengths(),
2648                           sizeof(int),numberColumns_,fp);
2649    if (numberWritten!=numberColumns_)
2650      return 1;
2651    // finished
2652    fclose(fp);
2653    return 0;
2654  } else {
2655    return -1;
2656  }
2657}
2658
2659int inDoubleArray(double * &array, int length, FILE * fp)
2660{
2661  int numberRead;
2662  int length2;
2663  numberRead = fread(&length2,sizeof(int),1,fp);
2664  if (numberRead!=1)
2665    return 1;
2666  if (length2) {
2667    // lengths must match
2668    if (length!=length2)
2669      return 2;
2670    array = new double[length];
2671    numberRead = fread(array,sizeof(double),length,fp);
2672    if (numberRead!=length)
2673      return 1;
2674  } 
2675  return 0;
2676}
2677/* Restore model from file, returns 0 if success,
2678   deletes current model */
2679int 
2680ClpSimplex::restoreModel(const char * fileName)
2681{
2682  FILE * fp = fopen(fileName,"rb");
2683  if (fp) {
2684    // Get rid of current model
2685    ClpModel::gutsOfDelete();
2686    gutsOfDelete(0);
2687    int i;
2688    for (i=0;i<6;i++) {
2689      rowArray_[i]=NULL;
2690      columnArray_[i]=NULL;
2691    }
2692    // get an empty factorization so we can set tolerances etc
2693    factorization_ = new ClpFactorization();
2694    // Say sparse
2695    factorization_->sparseThreshold(1);
2696    Clp_scalars scalars;
2697    CoinBigIndex numberRead;
2698
2699    // get scalars
2700    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
2701    if (numberRead!=1)
2702      return 1;
2703    // Fill in scalars
2704    optimizationDirection_ = scalars.optimizationDirection;
2705    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
2706    objectiveValue_ = scalars.objectiveValue;
2707    dualBound_ = scalars.dualBound;
2708    dualTolerance_ = scalars.dualTolerance;
2709    primalTolerance_ = scalars.primalTolerance;
2710    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
2711    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
2712    infeasibilityCost_ = scalars.infeasibilityCost;
2713    numberRows_ = scalars.numberRows;
2714    numberColumns_ = scalars.numberColumns;
2715    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
2716    numberIterations_ = scalars.numberIterations;
2717    problemStatus_ = scalars.problemStatus;
2718    setMaximumIterations(scalars.maximumIterations);
2719    lengthNames_ = scalars.lengthNames;
2720    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
2721    numberDualInfeasibilitiesWithoutFree_
2722      = scalars.numberDualInfeasibilitiesWithoutFree;
2723    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
2724    numberRefinements_ = scalars.numberRefinements;
2725    scalingFlag_ = scalars.scalingFlag;
2726    algorithm_ = scalars.algorithm;
2727    specialOptions_ = scalars.specialOptions;
2728    // strings
2729    CoinBigIndex length;
2730    for (i=0;i<ClpLastStrParam;i++) {
2731      numberRead = fread(&length,sizeof(int),1,fp);
2732      if (numberRead!=1)
2733        return 1;
2734      char * array = new char[length+1];
2735      numberRead = fread(array,length,1,fp);
2736      if (numberRead!=1)
2737        return 1;
2738      array[length]='\0';
2739      strParam_[i]=array;
2740      delete [] array;
2741    }
2742    // arrays - in no particular order
2743    if (inDoubleArray(rowActivity_,numberRows_,fp))
2744        return 1;
2745    if (inDoubleArray(columnActivity_,numberColumns_,fp))
2746        return 1;
2747    if (inDoubleArray(dual_,numberRows_,fp))
2748        return 1;
2749    if (inDoubleArray(reducedCost_,numberColumns_,fp))
2750        return 1;
2751    if (inDoubleArray(rowLower_,numberRows_,fp))
2752        return 1;
2753    if (inDoubleArray(rowUpper_,numberRows_,fp))
2754        return 1;
2755    double * objective;
2756    if (inDoubleArray(objective,numberColumns_,fp))
2757        return 1;
2758    delete objective_;
2759    objective_ = new ClpLinearObjective(objective,numberColumns_);
2760    delete [] objective;
2761    if (inDoubleArray(rowObjective_,numberRows_,fp))
2762        return 1;
2763    if (inDoubleArray(columnLower_,numberColumns_,fp))
2764        return 1;
2765    if (inDoubleArray(columnUpper_,numberColumns_,fp))
2766        return 1;
2767    if (problemStatus_==1) {
2768      if (inDoubleArray(ray_,numberRows_,fp))
2769        return 1;
2770    } else if (problemStatus_==2) {
2771      if (inDoubleArray(ray_,numberColumns_,fp))
2772        return 1;
2773    } else {
2774      // ray should be null
2775      numberRead = fread(&length,sizeof(int),1,fp);
2776      if (numberRead!=1)
2777        return 1;
2778      if (length)
2779        return 2;
2780    }
2781    delete [] status_;
2782    status_=NULL;
2783    // status region
2784    numberRead = fread(&length,sizeof(int),1,fp);
2785    if (numberRead!=1)
2786        return 1;
2787    if (length) {
2788      if (length!=numberRows_+numberColumns_)
2789        return 1;
2790      status_ = new char unsigned[length];
2791      numberRead = fread(status_,sizeof(char),length, fp);
2792      if (numberRead!=length)
2793        return 1;
2794    }
2795    if (lengthNames_) {
2796      char * array = new char[lengthNames_+1];
2797      rowNames_.resize(0);
2798      for (i=0;i<numberRows_;i++) {
2799        numberRead = fread(&length,sizeof(int),1,fp);
2800        if (numberRead!=1)
2801          return 1;
2802        numberRead = fread(array,length,1,fp);
2803        if (numberRead!=1)
2804          return 1;
2805        rowNames_[i]=array;
2806      }
2807      columnNames_.resize(0);
2808      for (i=0;i<numberColumns_;i++) {
2809        numberRead = fread(&length,sizeof(int),1,fp);
2810        if (numberRead!=1)
2811          return 1;
2812        numberRead = fread(array,length,1,fp);
2813        if (numberRead!=1)
2814          return 1;
2815        columnNames_[i]=array;
2816      }
2817      delete [] array;
2818    }
2819    // Pivot choices
2820    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
2821    delete dualRowPivot_;
2822    switch ((scalars.dualPivotChoice&63)) {
2823    case 1:
2824      // Dantzig
2825      dualRowPivot_ = new ClpDualRowDantzig();
2826      break;
2827    case 2:
2828      // Steepest - use mode
2829      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
2830      break;
2831    default:
2832      abort();
2833    }
2834    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
2835    delete primalColumnPivot_;
2836    switch ((scalars.primalPivotChoice&63)) {
2837    case 1:
2838      // Dantzig
2839      primalColumnPivot_ = new ClpPrimalColumnDantzig();
2840      break;
2841    case 2:
2842      // Steepest - use mode
2843      primalColumnPivot_
2844        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
2845      break;
2846    default:
2847      abort();
2848    }
2849    assert(scalars.matrixStorageChoice==1);
2850    delete matrix_;
2851    // get arrays
2852    numberRead = fread(&length,sizeof(int),1,fp);
2853    if (numberRead!=1)
2854      return 1;
2855    double * elements = new double[length];
2856    int * indices = new int[length];
2857    CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
2858    int * lengths = new int[numberColumns_];
2859    numberRead = fread(elements, sizeof(double),length,fp);
2860    if (numberRead!=length)
2861      return 1;
2862    numberRead = fread(indices, sizeof(int),length,fp);
2863    if (numberRead!=length)
2864      return 1;
2865    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
2866    if (numberRead!=numberColumns_)
2867      return 1;
2868    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
2869    if (numberRead!=numberColumns_)
2870      return 1;
2871    // assign matrix
2872    CoinPackedMatrix * matrix = new CoinPackedMatrix();
2873    matrix->assignMatrix(true, numberRows_, numberColumns_,
2874                         length, elements, indices, starts, lengths);
2875    // and transfer to Clp
2876    matrix_ = new ClpPackedMatrix(matrix);
2877    // finished
2878    fclose(fp);
2879    return 0;
2880  } else {
2881    return -1;
2882  }
2883  return 0;
2884}
2885// value of incoming variable (in Dual)
2886double 
2887ClpSimplex::valueIncomingDual() const
2888{
2889  // Need value of incoming for list of infeasibilities as may be infeasible
2890  double valueIncoming = (dualOut_/alpha_)*directionOut_;
2891  if (directionIn_==-1)
2892    valueIncoming = upperIn_-valueIncoming;
2893  else
2894    valueIncoming = lowerIn_-valueIncoming;
2895  return valueIncoming;
2896}
2897//#############################################################################
2898
2899ClpSimplexProgress::ClpSimplexProgress () 
2900{
2901  int i;
2902  for (i=0;i<CLP_PROGRESS;i++) {
2903    objective_[i] = 0.0;
2904    infeasibility_[i] = -1.0; // set to an impossible value
2905    numberInfeasibilities_[i]=-1; 
2906  }
2907  numberTimes_ = 0;
2908  numberBadTimes_ = 0;
2909  model_ = NULL;
2910}
2911
2912
2913//-----------------------------------------------------------------------------
2914
2915ClpSimplexProgress::~ClpSimplexProgress ()
2916{
2917}
2918// Copy constructor.
2919ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2920{
2921  int i;
2922  for (i=0;i<CLP_PROGRESS;i++) {
2923    objective_[i] = rhs.objective_[i];
2924    infeasibility_[i] = rhs.infeasibility_[i];
2925    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2926  }
2927  numberTimes_ = rhs.numberTimes_;
2928  numberBadTimes_ = rhs.numberBadTimes_;
2929  model_ = rhs.model_;
2930}
2931// Copy constructor.from model
2932ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2933{
2934  int i;
2935  for (i=0;i<CLP_PROGRESS;i++) {
2936    objective_[i] = 0.0;
2937    infeasibility_[i] = -1.0; // set to an impossible value
2938    numberInfeasibilities_[i]=-1; 
2939  }
2940  numberTimes_ = 0;
2941  numberBadTimes_ = 0;
2942  model_ = model;
2943}
2944// Assignment operator. This copies the data
2945ClpSimplexProgress & 
2946ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2947{
2948  if (this != &rhs) {
2949    int i;
2950    for (i=0;i<CLP_PROGRESS;i++) {
2951      objective_[i] = rhs.objective_[i];
2952      infeasibility_[i] = rhs.infeasibility_[i];
2953      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2954    }
2955    numberTimes_ = rhs.numberTimes_;
2956    numberBadTimes_ = rhs.numberBadTimes_;
2957    model_ = rhs.model_;
2958  }
2959  return *this;
2960}
2961// Seems to be something odd about exact comparison of doubles on linux
2962static bool equalDouble(double value1, double value2)
2963{
2964  assert(sizeof(unsigned int)*2==sizeof(double));
2965  unsigned int *i1 = (unsigned int *) &value1;
2966  unsigned int *i2 = (unsigned int *) &value2;
2967  return (i1[0]==i2[0]&&i1[1]==i2[1]);
2968}
2969int
2970ClpSimplexProgress::looping()
2971{
2972  if (!model_)
2973    return -1;
2974  double objective = model_->objectiveValue();
2975  double infeasibility;
2976  int numberInfeasibilities;
2977  if (model_->algorithm()<0) {
2978    // dual
2979    infeasibility = model_->sumPrimalInfeasibilities();
2980    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2981  } else {
2982    //primal
2983    infeasibility = model_->sumDualInfeasibilities();
2984    numberInfeasibilities = model_->numberDualInfeasibilities();
2985  }
2986  int i;
2987  int numberMatched=0;
2988  int matched=0;
2989
2990  for (i=0;i<CLP_PROGRESS;i++) {
2991    bool matchedOnObjective = equalDouble(objective,objective_[i]);
2992    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
2993    bool matchedOnInfeasibilities = 
2994      (numberInfeasibilities==numberInfeasibilities_[i]);
2995
2996    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
2997      matched |= (1<<i);
2998      numberMatched++;
2999      // here mainly to get over compiler bug?
3000      if (model_->messageHandler()->logLevel()>10)
3001        printf("%d %d %d %d %d loop check\n",i,numberMatched,
3002               matchedOnObjective, matchedOnInfeasibility, 
3003               matchedOnInfeasibilities);
3004    }
3005    if (i) {
3006      objective_[i-1] = objective_[i];
3007      infeasibility_[i-1] = infeasibility_[i];
3008      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
3009    }
3010  }
3011  objective_[CLP_PROGRESS-1] = objective;
3012  infeasibility_[CLP_PROGRESS-1] = infeasibility;
3013  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
3014  if (model_->progressFlag())
3015    numberMatched=0;
3016  numberTimes_++;
3017  if (numberTimes_<10)
3018    numberMatched=0;
3019  // skip if just last time as may be checking something
3020  if (matched==(1<<(CLP_PROGRESS-1)))
3021    numberMatched=0;
3022  if (numberMatched) {
3023    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
3024      <<numberMatched
3025      <<matched
3026      <<numberTimes_
3027      <<CoinMessageEol;
3028    numberBadTimes_++;
3029    if (numberBadTimes_<10) {
3030      // make factorize every iteration
3031      model_->forceFactorization(1);
3032      if (model_->algorithm()<0) {
3033        // dual - change tolerance
3034        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
3035        // if infeasible increase dual bound
3036        if (model_->dualBound()<1.0e19) {
3037          model_->setDualBound(model_->dualBound()*1.1);
3038        }
3039      } else {
3040        // primal - change tolerance
3041        model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
3042        // if infeasible increase infeasibility cost
3043        if (model_->infeasibilityCost()<1.0e19) {
3044          model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
3045        }
3046      }
3047      return -2;
3048    } else {
3049      model_->messageHandler()->message(CLP_LOOP,model_->messages())
3050        <<CoinMessageEol;
3051      // look at solution and maybe declare victory
3052      if (infeasibility<1.0e-4)
3053        return 0;
3054      else
3055        return 3;
3056    }
3057  }
3058  return -1;
3059}
3060// Sanity check on input data - returns true if okay
3061bool 
3062ClpSimplex::sanityCheck()
3063{
3064  // bad if empty
3065  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
3066    handler_->message(CLP_EMPTY_PROBLEM,messages_)
3067      <<numberRows_
3068      <<numberColumns_
3069      <<matrix_->getNumElements()
3070      <<CoinMessageEol;
3071    problemStatus_=4;
3072    return false;
3073  }
3074  int numberBad ;
3075  double largestBound, smallestBound, minimumGap;
3076  double smallestObj, largestObj;
3077  int firstBad;
3078  int modifiedBounds=0;
3079  int i;
3080  numberBad=0;
3081  firstBad=-1;
3082  minimumGap=1.0e100;
3083  smallestBound=1.0e100;
3084  largestBound=0.0;
3085  smallestObj=1.0e100;
3086  largestObj=0.0;
3087  // If bounds are too close - fix
3088  double fixTolerance = 1.1*primalTolerance_;
3089  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
3090    double value;
3091    value = fabs(cost_[i]);
3092    if (value>1.0e50) {
3093      numberBad++;
3094      if (firstBad<0)
3095        firstBad=i;
3096    } else if (value) {
3097      if (value>largestObj)
3098        largestObj=value;
3099      if (value<smallestObj)
3100        smallestObj=value;
3101    }
3102    value=upper_[i]-lower_[i];
3103    if (value<-primalTolerance_) {
3104      numberBad++;
3105      if (firstBad<0)
3106        firstBad=i;
3107    } else if (value<=fixTolerance) {
3108      if (value) {
3109        // modify
3110        upper_[i] = lower_[i];
3111        modifiedBounds++;
3112      }
3113    } else {
3114      if (value<minimumGap)
3115        minimumGap=value;
3116    }
3117    if (lower_[i]>-1.0e100&&lower_[i]) {
3118      value = fabs(lower_[i]);
3119      if (value>largestBound)
3120        largestBound=value;
3121      if (value<smallestBound)
3122        smallestBound=value;
3123    }
3124    if (upper_[i]<1.0e100&&upper_[i]) {
3125      value = fabs(upper_[i]);
3126      if (value>largestBound)
3127        largestBound=value;
3128      if (value<smallestBound)
3129        smallestBound=value;
3130    }
3131  }
3132  if (largestBound)
3133    handler_->message(CLP_RIMSTATISTICS3,messages_)
3134      <<smallestBound
3135      <<largestBound
3136      <<minimumGap
3137      <<CoinMessageEol;
3138  minimumGap=1.0e100;
3139  smallestBound=1.0e100;
3140  largestBound=0.0;
3141  for (i=0;i<numberColumns_;i++) {
3142    double value;
3143    value = fabs(cost_[i]);
3144    if (value>1.0e50) {
3145      numberBad++;
3146      if (firstBad<0)
3147        firstBad=i;
3148    } else if (value) {
3149      if (value>largestObj)
3150        largestObj=value;
3151      if (value<smallestObj)
3152        smallestObj=value;
3153    }
3154    value=upper_[i]-lower_[i];
3155    if (value<-primalTolerance_) {
3156      numberBad++;
3157      if (firstBad<0)
3158        firstBad=i;
3159    } else if (value<=fixTolerance) {
3160      if (value) {
3161        // modify
3162        upper_[i] = lower_[i];
3163        modifiedBounds++;
3164      }
3165    } else {
3166      if (value<minimumGap)
3167        minimumGap=value;
3168    }
3169    if (lower_[i]>-1.0e100&&lower_[i]) {
3170      value = fabs(lower_[i]);
3171      if (value>largestBound)
3172        largestBound=value;
3173      if (value<smallestBound)
3174        smallestBound=value;
3175    }
3176    if (upper_[i]<1.0e100&&upper_[i]) {
3177      value = fabs(upper_[i]);
3178      if (value>largestBound)
3179        largestBound=value;
3180      if (value<smallestBound)
3181        smallestBound=value;
3182    }
3183  }
3184  char rowcol[]={'R','C'};
3185  if (numberBad) {
3186    handler_->message(CLP_BAD_BOUNDS,messages_)
3187      <<numberBad
3188      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
3189      <<CoinMessageEol;
3190    problemStatus_=4;
3191    return false;
3192  }
3193  if (modifiedBounds)
3194    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
3195      <<modifiedBounds
3196      <<CoinMessageEol;
3197  handler_->message(CLP_RIMSTATISTICS1,messages_)
3198    <<smallestObj
3199    <<largestObj
3200    <<CoinMessageEol;
3201  if (largestBound)
3202    handler_->message(CLP_RIMSTATISTICS2,messages_)
3203      <<smallestBound
3204      <<largestBound
3205      <<minimumGap
3206      <<CoinMessageEol;
3207  return true;
3208}
3209// Set up status array (for OsiClp)
3210void 
3211ClpSimplex::createStatus() 
3212{
3213  if(!status_)
3214    status_ = new unsigned char [numberColumns_+numberRows_];
3215  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3216  int i;
3217  // set column status to one nearest zero
3218  for (i=0;i<numberColumns_;i++) {
3219#if 0
3220    if (columnLower_[i]>=0.0) {
3221      setColumnStatus(i,atLowerBound);
3222    } else if (columnUpper_[i]<=0.0) {
3223      setColumnStatus(i,atUpperBound);
3224    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
3225      // free
3226      setColumnStatus(i,isFree);
3227    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
3228      setColumnStatus(i,atLowerBound);
3229    } else {
3230      setColumnStatus(i,atUpperBound);
3231    }
3232#else
3233    setColumnStatus(i,atLowerBound);
3234#endif
3235  }
3236  for (i=0;i<numberRows_;i++) {
3237    setRowStatus(i,basic);
3238  }
3239}
3240/* Loads a problem (the constraints on the
3241   rows are given by lower and upper bounds). If a pointer is 0 then the
3242   following values are the default:
3243   <ul>
3244   <li> <code>colub</code>: all columns have upper bound infinity
3245   <li> <code>collb</code>: all columns have lower bound 0
3246   <li> <code>rowub</code>: all rows have upper bound infinity
3247   <li> <code>rowlb</code>: all rows have lower bound -infinity
3248   <li> <code>obj</code>: all variables have 0 objective coefficient
3249   </ul>
3250*/
3251void 
3252ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
3253                    const double* collb, const double* colub,   
3254                    const double* obj,
3255                    const double* rowlb, const double* rowub,
3256                    const double * rowObjective)
3257{
3258  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3259                        rowObjective);
3260  createStatus();
3261}
3262void 
3263ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
3264                    const double* collb, const double* colub,   
3265                    const double* obj,
3266                    const double* rowlb, const double* rowub,
3267                    const double * rowObjective)
3268{
3269  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3270                        rowObjective);
3271  createStatus();
3272}
3273
3274/* Just like the other loadProblem() method except that the matrix is
3275   given in a standard column major ordered format (without gaps). */
3276void 
3277ClpSimplex::loadProblem (  const int numcols, const int numrows,
3278                    const CoinBigIndex* start, const int* index,
3279                    const double* value,
3280                    const double* collb, const double* colub,   
3281                    const double* obj,
3282                    const double* rowlb, const double* rowub,
3283                    const double * rowObjective)
3284{
3285  ClpModel::loadProblem(numcols, numrows, start, index, value,
3286                          collb, colub, obj, rowlb, rowub,
3287                          rowObjective);
3288  createStatus();
3289}
3290void 
3291ClpSimplex::loadProblem (  const int numcols, const int numrows,
3292                           const CoinBigIndex* start, const int* index,
3293                           const double* value,const int * length,
3294                           const double* collb, const double* colub,   
3295                           const double* obj,
3296                           const double* rowlb, const double* rowub,
3297                           const double * rowObjective)
3298{
3299  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3300                          collb, colub, obj, rowlb, rowub,
3301                          rowObjective);
3302  createStatus();
3303}
3304// Read an mps file from the given filename
3305int 
3306ClpSimplex::readMps(const char *filename,
3307            bool keepNames,
3308            bool ignoreErrors)
3309{
3310  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3311  createStatus();
3312  return status;
3313}
3314// Just check solution (for external use)
3315void 
3316ClpSimplex::checkSolution()
3317{
3318  // put in standard form
3319  createRim(7+8+16);
3320  dualTolerance_=dblParam_[ClpDualTolerance];
3321  primalTolerance_=dblParam_[ClpPrimalTolerance];
3322  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
3323  checkDualSolution();
3324#ifdef CLP_DEBUG
3325  int i;
3326  double value=0.0;
3327  for (i=0;i<numberRows_+numberColumns_;i++)
3328    value += dj_[i]*solution_[i];
3329  printf("dual value %g, primal %g\n",value,objectiveValue());
3330#endif
3331  // release extra memory
3332  deleteRim(0);
3333}
3334/* Crash - at present just aimed at dual, returns
3335   -2 if dual preferred and crash basis created
3336   -1 if dual preferred and all slack basis preferred
3337   0 if basis going in was not all slack
3338   1 if primal preferred and all slack basis preferred
3339   2 if primal preferred and crash basis created.
3340   
3341   if gap between bounds <="gap" variables can be flipped
3342   
3343   If "pivot" is
3344   0 No pivoting (so will just be choice of algorithm)
3345   1 Simple pivoting e.g. gub
3346   2 Mini iterations
3347*/
3348int 
3349ClpSimplex::crash(double gap,int pivot)
3350{
3351  assert(!rowObjective_); // not coded
3352  int iColumn;
3353  int numberBad=0;
3354  int numberBasic=0;
3355  double dualTolerance=dblParam_[ClpDualTolerance];
3356  //double primalTolerance=dblParam_[ClpPrimalTolerance];
3357
3358  // If no basis then make all slack one
3359  if (!status_)
3360    createStatus();
3361 
3362  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3363    if (getColumnStatus(iColumn)==basic)
3364      numberBasic++;
3365  }
3366  if (numberBasic) {
3367    // not all slack
3368    return 0;
3369  } else {
3370    double * dj = new double [numberColumns_];
3371    double * solution = columnActivity_;
3372    const double * linearObjective = objective();
3373    //double objectiveValue=0.0;
3374    int iColumn;
3375    for (iColumn=0;iColumn<numberColumns_;iColumn++)
3376      dj[iColumn] = optimizationDirection_*linearObjective[iColumn];
3377    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3378      // assume natural place is closest to zero
3379      double lowerBound = columnLower_[iColumn];
3380      double upperBound = columnUpper_[iColumn];
3381      if (lowerBound>-1.0e20||upperBound<1.0e20) {
3382        bool atLower;
3383        if (fabs(upperBound)<fabs(lowerBound)) {
3384          atLower=false;
3385          setColumnStatus(iColumn,atUpperBound);
3386          solution[iColumn]=upperBound;
3387        } else {
3388          atLower=true;
3389          setColumnStatus(iColumn,atLowerBound);
3390          solution[iColumn]=lowerBound;
3391        }
3392        if (dj[iColumn]<0.0) {
3393          // should be at upper bound
3394          if (atLower) {
3395            // can we flip
3396            if (upperBound-lowerBound<=gap) {
3397              columnActivity_[iColumn]=upperBound;
3398              setColumnStatus(iColumn,atUpperBound);
3399            } else if (dj[iColumn]<-dualTolerance) {
3400              numberBad++;
3401            }
3402          }
3403        } else if (dj[iColumn]>0.0) {
3404          // should be at lower bound
3405          if (!atLower) {
3406            // can we flip
3407            if (upperBound-lowerBound<=gap) {
3408              columnActivity_[iColumn]=lowerBound;
3409              setColumnStatus(iColumn,atLowerBound);
3410            } else if (dj[iColumn]>dualTolerance) {
3411              numberBad++;
3412            }
3413          }
3414        }
3415      } else {
3416        // free
3417        setColumnStatus(iColumn,isFree);
3418        if (fabs(dj[iColumn])>dualTolerance) 
3419          numberBad++;
3420      }
3421    }
3422    if (numberBad||pivot) {
3423      if (!pivot) {
3424        delete [] dj;
3425        return 1;
3426      } else {
3427        // see if can be made dual feasible with gubs etc
3428        double * pi = new double[numberRows_];
3429        memset (pi,0,numberRows_*sizeof(double));
3430        int * way = new int[numberColumns_];
3431        int numberIn = 0;
3432       
3433        // Get column copy
3434        CoinPackedMatrix * columnCopy = matrix();
3435        // Get a row copy in standard format
3436        CoinPackedMatrix copy;
3437        copy.reverseOrderedCopyOf(*columnCopy);
3438        // get matrix data pointers
3439        const int * column = copy.getIndices();
3440        const CoinBigIndex * rowStart = copy.getVectorStarts();
3441        const int * rowLength = copy.getVectorLengths(); 
3442        const double * elementByRow = copy.getElements();
3443        //const int * row = columnCopy->getIndices();
3444        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3445        //const int * columnLength = columnCopy->getVectorLengths();
3446        //const double * element = columnCopy->getElements();
3447       
3448       
3449        // if equality row and bounds mean artificial in basis bad
3450        // then do anyway
3451       
3452        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3453          // - if we want to reduce dj, + if we want to increase
3454          int thisWay = 100;
3455          double lowerBound = columnLower_[iColumn];
3456          double upperBound = columnUpper_[iColumn];
3457          if (upperBound>lowerBound) {
3458            switch(getColumnStatus(iColumn)) {
3459             
3460            case basic:
3461              thisWay=0;
3462            case ClpSimplex::isFixed:
3463              break;
3464            case isFree:
3465            case superBasic:
3466              if (dj[iColumn]<-dualTolerance) 
3467                thisWay = 1;
3468              else if (dj[iColumn]>dualTolerance) 
3469                thisWay = -1;
3470              else
3471                thisWay =0;
3472              break;
3473            case atUpperBound:
3474              if (dj[iColumn]>dualTolerance) 
3475                thisWay = -1;
3476              else if (dj[iColumn]<-dualTolerance) 
3477                thisWay = -3;
3478              else
3479                thisWay = -2;
3480              break;
3481            case atLowerBound:
3482              if (dj[iColumn]<-dualTolerance) 
3483                thisWay = 1;
3484              else if (dj[iColumn]>dualTolerance) 
3485                thisWay = 3;
3486              else
3487                thisWay = 2;
3488              break;
3489            }
3490          }
3491          way[iColumn] = thisWay;
3492        }
3493        /*if (!numberBad)
3494          printf("Was dual feasible before passes - rows %d\n",
3495          numberRows_);*/
3496        int lastNumberIn = -100000;
3497        int numberPasses=5;
3498        while (numberIn>lastNumberIn+numberRows_/100) {
3499          lastNumberIn = numberIn;
3500          // we need to maximize chance of doing good
3501          int iRow;
3502          for (iRow=0;iRow<numberRows_;iRow++) {
3503            double lowerBound = rowLower_[iRow];
3504            double upperBound = rowUpper_[iRow];
3505            if (getRowStatus(iRow)==basic) {
3506              // see if we can find a column to pivot on
3507              int j;
3508              // down is amount pi can go down
3509              double maximumDown = COIN_DBL_MAX;
3510              double maximumUp = COIN_DBL_MAX;
3511              double minimumDown =0.0;
3512              double minimumUp =0.0;
3513              int iUp=-1;
3514              int iDown=-1;
3515              int iUpB=-1;
3516              int iDownB=-1;
3517              if (lowerBound<-1.0e20)
3518                maximumUp = -1.0;
3519              if (upperBound>1.0e20)
3520                maximumDown = -1.0;
3521              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3522                int iColumn = column[j];
3523                double value = elementByRow[j];
3524                double djValue = dj[iColumn];
3525                /* way -
3526                   -3 - okay at upper bound with negative dj
3527                   -2 - marginal at upper bound with zero dj - can only decrease
3528                   -1 - bad at upper bound
3529                   0 - we can never pivot on this row
3530                   1 - bad at lower bound
3531                   2 - marginal at lower bound with zero dj - can only increase
3532                   3 - okay at lower bound with positive dj
3533                   100 - fine we can just ignore
3534                */
3535                if (way[iColumn]!=100) {
3536                  switch(way[iColumn]) {
3537                   
3538                  case -3:
3539                    if (value>0.0) {
3540                      if (maximumDown*value>-djValue) {
3541                        maximumDown = -djValue/value;
3542                        iDown = iColumn;
3543                      }
3544                    } else {
3545                      if (-maximumUp*value>-djValue) {
3546                        maximumUp = djValue/value;
3547                        iUp = iColumn;
3548                      }
3549                    }
3550                    break;
3551                  case -2:
3552                    if (value>0.0) {
3553                      maximumDown = 0.0;
3554                    } else {
3555                      maximumUp = 0.0;
3556                    }
3557                    break;
3558                  case -1:
3559                    // see if could be satisfied
3560                    // dj value > 0
3561                    if (value>0.0) {
3562                      maximumDown=0.0;
3563                      if (maximumUp*value<djValue-dualTolerance) {
3564                        maximumUp = 0.0; // would improve but not enough
3565                      } else {
3566                        if (minimumUp*value<djValue) {
3567                          minimumUp = djValue/value;
3568                          iUpB = iColumn;
3569                        }
3570                      }
3571                    } else {
3572                      maximumUp=0.0;
3573                      if (-maximumDown*value<djValue-dualTolerance) {
3574                        maximumDown = 0.0; // would improve but not enough
3575                      } else {
3576                        if (-minimumDown*value<djValue) {
3577                          minimumDown = -djValue/value;
3578                          iDownB = iColumn;
3579                        }
3580                      }
3581                    }
3582                   
3583                    break;
3584                  case 0:
3585                    maximumDown = -1.0;
3586                    maximumUp=-1.0;
3587                    break;
3588                  case 1:
3589                    // see if could be satisfied
3590                    // dj value < 0
3591                    if (value>0.0) {
3592                      maximumUp=0.0;
3593                      if (maximumDown*value<-djValue-dualTolerance) {
3594                        maximumDown = 0.0; // would improve but not enough
3595                      } else {
3596                        if (minimumDown*value<-djValue) {
3597                          minimumDown = -djValue/value;
3598                          iDownB = iColumn;
3599                        }
3600                      }
3601                    } else {
3602                      maximumDown=0.0;
3603                      if (-maximumUp*value<-djValue-dualTolerance) {
3604                        maximumUp = 0.0; // would improve but not enough
3605                      } else {
3606                        if (-minimumUp*value<-djValue) {
3607                          minimumUp = djValue/value;
3608                          iUpB = iColumn;
3609                        }
3610                      }
3611                    }
3612                   
3613                    break;
3614                  case 2:
3615                    if (value>0.0) {
3616                      maximumUp = 0.0;
3617                    } else {
3618                      maximumDown = 0.0;
3619                    }
3620                   
3621                    break;
3622                  case 3:
3623                    if (value>0.0) {
3624                      if (maximumUp*value>djValue) {
3625                        maximumUp = djValue/value;
3626                        iUp = iColumn;
3627                      }
3628                    } else {
3629                      if (-maximumDown*value>djValue) {
3630                        maximumDown = -djValue/value;
3631                        iDown = iColumn;
3632                      }
3633                    }
3634                   
3635                    break;
3636                  default:
3637                    break;
3638                  }
3639                }
3640              }
3641              if (iUpB>=0)
3642                iUp=iUpB;
3643              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
3644                iUp=-1;
3645              if (iDownB>=0)
3646                iDown=iDownB;
3647              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
3648                iDown=-1;
3649              if (iUp>=0||iDown>=0) {
3650                // do something
3651                if (iUp>=0&&iDown>=0) {
3652                  if (maximumDown>maximumUp)
3653                    iUp=-1;
3654                }
3655                double change;
3656                int kColumn;
3657                if (iUp>=0) {
3658                  kColumn=iUp;
3659                  change=maximumUp;
3660                  // just do minimum if was dual infeasible
3661                  // ? only if maximum large?
3662                  if (minimumUp>0.0)
3663                    change=minimumUp;
3664                  setRowStatus(iRow,atUpperBound);
3665                } else {
3666                  kColumn=iDown;
3667                  change=-maximumDown;
3668                  // just do minimum if was dual infeasible
3669                  // ? only if maximum large?
3670                  if (minimumDown>0.0)
3671                    change=-minimumDown;
3672                  setRowStatus(iRow,atLowerBound);
3673                }
3674                assert (fabs(change)<1.0e20);
3675                setColumnStatus(kColumn,basic);
3676                numberIn++;
3677                pi[iRow]=change;
3678                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3679                  int iColumn = column[j];
3680                  double value = elementByRow[j];
3681                  double djValue = dj[iColumn]-change*value;
3682                  dj[iColumn]=djValue;
3683                  if (abs(way[iColumn])==1) {
3684                    numberBad--;
3685                    /*if (!numberBad)
3686                      printf("Became dual feasible at row %d out of %d\n",
3687                      iRow, numberRows_);*/
3688                    lastNumberIn=-1000000;
3689                  }
3690                  int thisWay = 100;
3691                  double lowerBound = columnLower_[iColumn];
3692                  double upperBound = columnUpper_[iColumn];
3693                  if (upperBound>lowerBound) {
3694                    switch(getColumnStatus(iColumn)) {
3695                     
3696                    case basic:
3697                      thisWay=0;
3698                    case isFixed:
3699                      break;
3700                    case isFree:
3701                    case superBasic:
3702                      if (djValue<-dualTolerance) 
3703                        thisWay = 1;
3704                      else if (djValue>dualTolerance) 
3705                        thisWay = -1;
3706                      else
3707                        { thisWay =0; abort();}
3708                      break;
3709                    case atUpperBound:
3710                      if (djValue>dualTolerance) 
3711                        { thisWay =-1; abort();}
3712                      else if (djValue<-dualTolerance) 
3713                        thisWay = -3;
3714                      else
3715                        thisWay = -2;
3716                      break;
3717                    case atLowerBound:
3718                      if (djValue<-dualTolerance) 
3719                        { thisWay =1; abort();}
3720                      else if (djValue>dualTolerance) 
3721                        thisWay = 3;
3722                      else
3723                        thisWay = 2;
3724                      break;
3725                    }
3726                  }
3727                  way[iColumn] = thisWay;
3728                }
3729              }
3730            }
3731          }
3732          if (numberIn==lastNumberIn||numberBad||pivot<2)
3733            break;
3734          if (!(--numberPasses))
3735            break;
3736          //printf("%d put in so far\n",numberIn);
3737        }
3738        // last attempt to flip
3739        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3740          double lowerBound = columnLower_[iColumn];
3741          double upperBound = columnUpper_[iColumn];
3742          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
3743            double djValue=dj[iColumn];
3744            switch(getColumnStatus(iColumn)) {
3745             
3746            case basic:
3747            case ClpSimplex::isFixed:
3748              break;
3749            case isFree:
3750            case superBasic:
3751              break;
3752            case atUpperBound:
3753              if (djValue>dualTolerance) {
3754                setColumnStatus(iColumn,atUpperBound);
3755                solution[iColumn]=upperBound;
3756              } 
3757              break;
3758            case atLowerBound:
3759              if (djValue<-dualTolerance) {
3760                setColumnStatus(iColumn,atUpperBound);
3761                solution[iColumn]=upperBound;
3762              }
3763              break;
3764            }
3765          }
3766        }
3767        delete [] pi;
3768        delete [] dj;
3769        delete [] way;
3770        handler_->message(CLP_CRASH,messages_)
3771          <<numberIn
3772          <<numberBad
3773          <<CoinMessageEol;
3774        return -1;
3775      }
3776    } else {
3777      delete [] dj;
3778      return -1;
3779    }
3780  }
3781}
3782int
3783ClpSimplex::nextSuperBasic()
3784{
3785  if (firstFree_>=0) {
3786    int returnValue=firstFree_;
3787    int iColumn=firstFree_+1;
3788    if (algorithm_>0) {
3789      // primal
3790      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
3791        if (getStatus(iColumn)==superBasic) {
3792          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
3793            solution_[iColumn]=lower_[iColumn];
3794            setStatus(iColumn,atLowerBound);
3795          } else if (fabs(solution_[iColumn]-upper_[iColumn])
3796                     <=primalTolerance_) {
3797            solution_[iColumn]=upper_[iColumn];
3798            setStatus(iColumn,atUpperBound);
3799          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
3800            setStatus(iColumn,isFree);
3801            break;
3802          } else {
3803            break;
3804          }
3805        }
3806      }
3807    } else {
3808      // dual
3809      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
3810        if (getStatus(iColumn)==isFree) 
3811          if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
3812            break;
3813      }
3814    }
3815    firstFree_ = iColumn;
3816    if (firstFree_==numberRows_+numberColumns_)
3817      firstFree_=-1;
3818    return returnValue;
3819  } else {
3820    return -1;
3821  }
3822}
3823/* Pivot in a variable and out a variable.  Returns 0 if okay,
3824   1 if inaccuracy forced re-factorization, -1 if would be singular.
3825   Also updates primal/dual infeasibilities.
3826   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
3827*/
3828int ClpSimplex::pivot()
3829{
3830  // scaling not allowed
3831  assert (!scalingFlag_);
3832  // assume In_ and Out_ are correct and directionOut_ set
3833  // (or In_ if flip
3834  lowerIn_ = lower_[sequenceIn_];
3835  valueIn_ = solution_[sequenceIn_];
3836  upperIn_ = upper_[sequenceIn_];
3837  dualIn_ = dj_[sequenceIn_];
3838  if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) {
3839    assert (pivotRow_>=0&&pivotRow_<numberRows_);
3840    assert (pivotVariable_[pivotRow_]==sequenceOut_);
3841    lowerOut_ = lower_[sequenceOut_];
3842    valueOut_ = solution_[sequenceOut_];
3843    upperOut_ = upper_[sequenceOut_];
3844    // for now assume primal is feasible (or in dual)
3845    dualOut_ = dj_[sequenceOut_];
3846    assert(fabs(dualOut_)<1.0e-6);
3847  } else {
3848    assert (pivotRow_<0);
3849  }
3850  bool roundAgain = true;
3851  int returnCode=0;
3852  while (roundAgain) {
3853    roundAgain=false;
3854    unpack(rowArray_[1]);
3855    factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
3856    // we are going to subtract movement from current basic
3857    double movement;
3858    // see where incoming will go to
3859    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
3860      // flip so go to bound
3861      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
3862    } else {
3863      // get where outgoing needs to get to
3864      double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
3865      // solutionOut_ - movement*alpha_ == outValue
3866      movement = (outValue-valueOut_)/alpha_;
3867      // set directionIn_ correctly
3868      directionIn_ = (movement>0) ? 1 :-1;
3869    }
3870    // update primal solution
3871    {
3872      int i;
3873      int * index = rowArray_[1]->getIndices();
3874      int number = rowArray_[1]->getNumElements();
3875      double * element = rowArray_[1]->denseVector();
3876      for (i=0;i<number;i++) {
3877        int ii = index[i];
3878        // get column
3879        ii = pivotVariable_[ii];
3880        solution_[ii] -= movement*element[i];
3881        element[i]=0.0;
3882      }
3883      // see where something went to
3884      if (sequenceOut_<0) {
3885        if (directionIn_<0) {
3886          assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
3887          solution_[sequenceIn_]=upperIn_;
3888        } else {
3889          assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
3890          solution_[sequenceIn_]=lowerIn_;
3891        }
3892      } else {
3893        if (directionOut_<0) {
3894          assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
3895          solution_[sequenceOut_]=upperOut_;
3896        } else {
3897          assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
3898          solution_[sequenceOut_]=lowerOut_;
3899        }
3900        solution_[sequenceIn_]=valueIn_+movement;
3901      }
3902    }   
3903    double objectiveChange = dualIn_*movement;
3904    // update duals
3905    if (pivotRow_>=0) {
3906      alpha_ = rowArray_[1]->denseVector()[pivotRow_];
3907      assert (fabs(alpha_)>1.0e-8);
3908      double multiplier = dualIn_/alpha_;
3909      rowArray_[0]->insert(pivotRow_,multiplier);
3910      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
3911      // put row of tableau in rowArray[0] and columnArray[0]
3912      matrix_->transposeTimes(this,-1.0,
3913                              rowArray_[0],columnArray_[1],columnArray_[0]);
3914      // update column djs
3915      int i;
3916      int * index = columnArray_[0]->getIndices();
3917      int number = columnArray_[0]->getNumElements();
3918      double * element = columnArray_[0]->denseVector();
3919      for (i=0;i<number;i++) {
3920        int ii = index[i];
3921        dj_[ii] += element[ii];
3922        element[ii]=0.0;
3923      }
3924      columnArray_[0]->setNumElements(0);
3925      // and row djs
3926      index = rowArray_[0]->getIndices();
3927      number = rowArray_[0]->getNumElements();
3928      element = rowArray_[0]->denseVector();
3929      for (i=0;i<number;i++) {
3930        int ii = index[i];
3931        dj_[ii+numberColumns_] += element[ii];
3932        dual_[ii] = dj_[ii+numberColumns_];
3933        element[ii]=0.0;
3934      }
3935      rowArray_[0]->setNumElements(0);
3936      // check incoming
3937      assert (fabs(dj_[sequenceIn_])<1.0e-6);
3938    }
3939   
3940    // if stable replace in basis
3941    int updateStatus = factorization_->replaceColumn(rowArray_[2],
3942                                                   pivotRow_,
3943                                                     alpha_);
3944    bool takePivot=true;
3945    // if no pivots, bad update but reasonable alpha - take and invert
3946    if (updateStatus==2&&
3947        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
3948      updateStatus=4;
3949    if (updateStatus==1||updateStatus==4) {
3950      // slight error
3951      if (factorization_->pivots()>5||updateStatus==4) {
3952        returnCode=-1;
3953      }
3954    } else if (updateStatus==2) {
3955      // major error
3956      rowArray_[1]->clear();
3957      takePivot=false;
3958      if (factorization_->pivots()) {
3959        // refactorize here
3960        statusOfProblem();
3961        roundAgain=true;
3962      } else {
3963        returnCode=1;
3964      }
3965    } else if (updateStatus==3) {
3966      // out of memory
3967      // increase space if not many iterations
3968      if (factorization_->pivots()<
3969          0.5*factorization_->maximumPivots()&&
3970          factorization_->pivots()<200)
3971        factorization_->areaFactor(
3972                                   factorization_->areaFactor() * 1.1);
3973      returnCode =-1; // factorize now
3974    }
3975    if (takePivot) {
3976      int save = algorithm_;
3977      // make simple so always primal
3978      algorithm_=1;
3979      housekeeping(objectiveChange);
3980      algorithm_=save;
3981    }
3982  }
3983  if (returnCode == -1) {
3984    // refactorize here
3985    statusOfProblem();
3986  } else {
3987    // just for now - recompute anyway
3988    gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
3989  }
3990  return returnCode;
3991}
3992
3993/* Pivot in a variable and choose an outgoing one.  Assumes primal
3994   feasible - will not go through a bound.  Returns step length in theta
3995   Returns ray in ray_ (or NULL if no pivot)
3996   Return codes as before but -1 means no acceptable pivot
3997*/
3998int ClpSimplex::primalPivotResult()
3999{
4000  assert (sequenceIn_>=0);
4001  valueIn_=solution_[sequenceIn_];
4002  lowerIn_=lower_[sequenceIn_];
4003  upperIn_=upper_[sequenceIn_];
4004  dualIn_=dj_[sequenceIn_];
4005
4006  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
4007  if (returnCode<0&&returnCode>-4) {
4008    return 0;
4009  } else {
4010    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
4011           returnCode);
4012    return -1;
4013  }
4014}
4015 
4016/* Pivot out a variable and choose an incoing one.  Assumes dual
4017   feasible - will not go through a reduced cost. 
4018   Returns step length in theta
4019   Returns ray in ray_ (or NULL if no pivot)
4020   Return codes as before but -1 means no acceptable pivot
4021*/
4022int 
4023ClpSimplex::dualPivotResult()
4024{
4025  return ((ClpSimplexDual *) this)->pivotResult();
4026}
4027// Common bits of coding for dual and primal
4028int 
4029ClpSimplex::startup(int ifValuesPass)
4030{
4031  // sanity check
4032  assert (numberRows_==matrix_->getNumRows());
4033  assert (numberColumns_==matrix_->getNumCols());
4034  // for moment all arrays must exist
4035  assert(columnLower_);
4036  assert(columnUpper_);
4037  assert(rowLower_);
4038  assert(rowUpper_);
4039
4040
4041  primalTolerance_=dblParam_[ClpPrimalTolerance];
4042  dualTolerance_=dblParam_[ClpDualTolerance];
4043
4044  // put in standard form (and make row copy)
4045  // create modifiable copies of model rim and do optional scaling
4046  bool goodMatrix=createRim(7+8+16,true);
4047
4048  if (goodMatrix) {
4049    // Model looks okay
4050    // Do initial factorization
4051    // and set certain stuff
4052    // We can either set increasing rows so ...IsBasic gives pivot row
4053    // or we can just increment iBasic one by one
4054    // for now let ...iBasic give pivot row
4055    factorization_->increasingRows(2);
4056    // row activities have negative sign
4057    factorization_->slackValue(-1.0);
4058    factorization_->zeroTolerance(1.0e-13);
4059    // Switch off dense
4060    int saveThreshold = factorization_->denseThreshold();
4061    factorization_->setDenseThreshold(0);
4062    numberIterations_=0;
4063
4064    // do perturbation if asked for
4065
4066    if (perturbation_<100) {
4067      if (algorithm_>0) {
4068        ((ClpSimplexPrimal *) this)->perturb();
4069      } else if (algorithm_<0) {
4070        ((ClpSimplexDual *) this)->perturb();
4071      }
4072    }
4073    // for primal we will change bounds using infeasibilityCost_
4074    if (nonLinearCost_==NULL&&algorithm_>0) {
4075      // get a valid nonlinear cost function
4076      delete nonLinearCost_;
4077      nonLinearCost_= new ClpNonLinearCost(this);
4078    }
4079   
4080    // loop round to clean up solution if values pass
4081    int numberThrownOut = -1;
4082    int totalNumberThrownOut=0;
4083    while(numberThrownOut) {
4084      int status = internalFactorize(0+10*ifValuesPass);
4085      if (status<0)
4086        return 1; // some error
4087      else
4088        totalNumberThrownOut+= status;
4089     
4090      // for this we need clean basis so it is after factorize
4091      numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
4092                                     NULL,NULL,
4093                                     ifValuesPass!=0);
4094      totalNumberThrownOut+= numberThrownOut;
4095     
4096    }
4097   
4098    if (totalNumberThrownOut)
4099      handler_->message(CLP_SINGULARITIES,messages_)
4100        <<totalNumberThrownOut
4101        <<CoinMessageEol;
4102    // Switch back dense
4103    factorization_->setDenseThreshold(saveThreshold);
4104   
4105    problemStatus_ = -1;
4106   
4107    // number of times we have declared optimality
4108    numberTimesOptimal_=0;
4109
4110    return 0;
4111  } else {
4112    // bad matrix
4113    return 2;
4114  }
4115   
4116}
4117
4118
4119void 
4120ClpSimplex::finish()
4121{
4122  // Get rid of some arrays and empty factorization
4123  deleteRim();
4124  assert(problemStatus_>=0&&problemStatus_<5);
4125  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
4126    <<objectiveValue()
4127    <<CoinMessageEol;
4128  factorization_->relaxAccuracyCheck(1.0);
4129  // get rid of any network stuff - could do more
4130  factorization_->cleanUp();
4131}
4132// Save data
4133ClpDataSave
4134ClpSimplex::saveData() const
4135{
4136  ClpDataSave saved;
4137  saved.dualBound_ = dualBound_;
4138  saved.infeasibilityCost_ = infeasibilityCost_;
4139  saved.sparseThreshold_ = factorization_->sparseThreshold();
4140  saved.perturbation_ = perturbation_;
4141  return saved;
4142}
4143// Restore data
4144void 
4145ClpSimplex::restoreData(ClpDataSave saved)
4146{
4147  factorization_->sparseThreshold(saved.sparseThreshold_);
4148  perturbation_ = saved.perturbation_;
4149  infeasibilityCost_ = saved.infeasibilityCost_;
4150  dualBound_ = saved.dualBound_;
4151}
4152/* Factorizes and returns true if optimal.  Used by user */
4153bool
4154ClpSimplex::statusOfProblem()
4155{
4156  // is factorization okay?
4157  assert (internalFactorize(1)==0);
4158  // put back original costs and then check
4159  // also move to work arrays
4160  createRim(4+32);
4161  abort();
4162  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
4163  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
4164  gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
4165  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
4166  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
4167  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
4168  deleteRim(-1);
4169  return (primalFeasible()&&dualFeasible());
4170}
4171/* Return model - updates any scalars */
4172void 
4173ClpSimplex::returnModel(ClpSimplex & otherModel)
4174{
4175  ClpModel::returnModel(otherModel);
4176  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
4177  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
4178  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
4179  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
4180  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
4181  otherModel.columnDualSequence_ = columnDualSequence_;
4182  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
4183  otherModel.rowDualSequence_ = rowDualSequence_;
4184  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
4185  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
4186  otherModel.largestPrimalError_ = largestPrimalError_;
4187  otherModel.largestDualError_ = largestDualError_;
4188  otherModel.largestSolutionError_ = largestSolutionError_;
4189  otherModel.alpha_ = alpha_;
4190  otherModel.theta_ = theta_;
4191  otherModel.lowerIn_ = lowerIn_;
4192  otherModel.valueIn_ = valueIn_;
4193  otherModel.upperIn_ = upperIn_;
4194  otherModel.dualIn_ = dualIn_;
4195  otherModel.sequenceIn_ = sequenceIn_;
4196  otherModel.directionIn_ = directionIn_;
4197  otherModel.lowerOut_ = lowerOut_;
4198  otherModel.valueOut_ = valueOut_;
4199  otherModel.upperOut_ = upperOut_;
4200  otherModel.dualOut_ = dualOut_;
4201  otherModel.sequenceOut_ = sequenceOut_;
4202  otherModel.directionOut_ = directionOut_;
4203  otherModel.pivotRow_ = pivotRow_;
4204  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
4205  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
4206  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
4207    numberDualInfeasibilitiesWithoutFree_;
4208  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
4209  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
4210  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
4211  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
4212  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
4213}
4214/* Constructs a non linear cost from list of non-linearities (columns only)
4215   First lower of each column is taken as real lower
4216   Last lower is taken as real upper and cost ignored
4217   
4218   Returns nonzero if bad data e.g. lowers not monotonic
4219*/
4220int 
4221ClpSimplex::createPiecewiseLinearCosts(const int * starts,
4222                                       const double * lower, const double * gradient)
4223{
4224  delete nonLinearCost_;
4225  // Set up feasible bounds and check monotonicity
4226  int iColumn;
4227  int returnCode=0;
4228
4229  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4230    int iIndex = starts[iColumn];
4231    int end = starts[iColumn+1]-1;
4232    columnLower_[iColumn] = lower[iIndex];
4233    columnUpper_[iColumn] = lower[end];
4234    double value = columnLower_[iColumn];
4235    iIndex++;
4236    for (;iIndex<end;iIndex++) {
4237      if (lower[iIndex]<value)
4238        returnCode++; // not monotonic
4239      value=lower[iIndex];
4240    }
4241  }
4242  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
4243  specialOptions_ |= 2; // say keep
4244  return returnCode;
4245}
Note: See TracBrowser for help on using the repository browser.