source:
trunk/Clp/src/ClpSimplex.cpp
@
1722
Last change on this file since 1722 was 1722, checked in by stefan, 9 years ago  



File size: 495.4 KB 
<
Rev  Line  

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