source: tags/clp-0-94-0/ClpSimplex.cpp @ 1814

Last change on this file since 1814 was 98, checked in by forrest, 17 years ago

Relax a tolerance

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