source: trunk/Clp/src/ClpSimplexOther.cpp @ 1869

Last change on this file since 1869 was 1869, checked in by forrest, 7 years ago

fix stupid new length error in mini presolve (thanks Valentin)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 355.3 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1869 2012-07-19 14:53:33Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9
10#include "CoinHelperFunctions.hpp"
11#include "ClpSimplexOther.hpp"
12#include "ClpSimplexDual.hpp"
13#include "ClpSimplexPrimal.hpp"
14#include "ClpEventHandler.hpp"
15#include "ClpHelperFunctions.hpp"
16#include "ClpFactorization.hpp"
17#include "ClpDualRowDantzig.hpp"
18#include "ClpNonLinearCost.hpp"
19#include "ClpDynamicMatrix.hpp"
20#include "CoinPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
22#include "CoinBuild.hpp"
23#include "CoinMpsIO.hpp"
24#include "CoinFloatEqual.hpp"
25#include "ClpMessage.hpp"
26#include <cfloat>
27#include <cassert>
28#include <string>
29#include <stdio.h>
30#include <iostream>
31#ifdef HAS_CILK
32#include <cilk/cilk.h>
33#else
34#define cilk_for for
35#define cilk_spawn
36#define cilk_sync
37#endif
38#ifdef INT_IS_8
39#define COIN_ANY_BITS_PER_INT 64
40#define COIN_ANY_SHIFT_PER_INT 6
41#define COIN_ANY_MASK_PER_INT 0x3f
42#else
43#define COIN_ANY_BITS_PER_INT 32
44#define COIN_ANY_SHIFT_PER_INT 5
45#define COIN_ANY_MASK_PER_INT 0x1f
46#endif
47/* Dual ranging.
48   This computes increase/decrease in cost for each given variable and corresponding
49   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
50   and numberColumns.. for artificials/slacks.
51   For non-basic variables the sequence number will be that of the non-basic variables.
52
53   Up to user to provide correct length arrays.
54
55*/
56void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
57                                  double * costIncreased, int * sequenceIncreased,
58                                  double * costDecreased, int * sequenceDecreased,
59                                  double * valueIncrease, double * valueDecrease)
60{
61     rowArray_[1]->clear();
62     columnArray_[1]->clear();
63     // long enough for rows+columns
64     assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_);
65     rowArray_[3]->clear();
66     int * backPivot = rowArray_[3]->getIndices();
67     int i;
68     for ( i = 0; i < numberRows_ + numberColumns_; i++) {
69          backPivot[i] = -1;
70     }
71     for (i = 0; i < numberRows_; i++) {
72          int iSequence = pivotVariable_[i];
73          backPivot[iSequence] = i;
74     }
75     // dualTolerance may be zero if from CBC.  In fact use that fact
76     bool inCBC = !dualTolerance_;
77     if (inCBC)
78          assert (integerType_);
79     dualTolerance_ = dblParam_[ClpDualTolerance];
80     double * arrayX = rowArray_[0]->denseVector();
81     for ( i = 0; i < numberCheck; i++) {
82          rowArray_[0]->clear();
83          //rowArray_[0]->checkClear();
84          //rowArray_[1]->checkClear();
85          //columnArray_[1]->checkClear();
86          columnArray_[0]->clear();
87          //columnArray_[0]->checkClear();
88          int iSequence = which[i];
89          if (iSequence < 0) {
90               costIncreased[i] = 0.0;
91               sequenceIncreased[i] = -1;
92               costDecreased[i] = 0.0;
93               sequenceDecreased[i] = -1;
94               continue;
95          }
96          double costIncrease = COIN_DBL_MAX;
97          double costDecrease = COIN_DBL_MAX;
98          int sequenceIncrease = -1;
99          int sequenceDecrease = -1;
100          if (valueIncrease) {
101               assert (valueDecrease);
102               valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
103               valueDecrease[i] = valueIncrease[i];
104          }
105
106          switch(getStatus(iSequence)) {
107
108          case basic: {
109               // non-trvial
110               // Get pivot row
111               int iRow = backPivot[iSequence];
112               assert (iRow >= 0);
113#ifndef COIN_FAC_NEW
114               double plusOne = 1.0;
115               rowArray_[0]->createPacked(1, &iRow, &plusOne);
116#else
117               rowArray_[0]->createOneUnpackedElement( iRow, 1.0);
118#endif
119               factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
120               // put row of tableau in rowArray[0] and columnArray[0]
121               matrix_->transposeTimes(this, -1.0,
122                                       rowArray_[0], columnArray_[1], columnArray_[0]);
123#if COIN_FAC_NEW
124               assert (!rowArray_[0]->packedMode());
125#endif
126               double alphaIncrease;
127               double alphaDecrease;
128               // do ratio test up and down
129               checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease,
130                               costDecrease, sequenceDecrease, alphaDecrease);
131               if (!inCBC) {
132                    if (valueIncrease) {
133                         if (sequenceIncrease >= 0)
134                              valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence);
135                         if (sequenceDecrease >= 0)
136                              valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence);
137                    }
138               } else {
139                    int number = rowArray_[0]->getNumElements();
140#if COIN_FAC_NEW
141                    const int * index = rowArray_[0]->getIndices();
142#endif
143                    double scale2 = 0.0;
144                    int j;
145                    for (j = 0; j < number; j++) {
146#ifndef COIN_FAC_NEW
147                         scale2 += arrayX[j] * arrayX[j];
148#else
149                         int iRow=index[j];
150                         scale2 += arrayX[iRow] * arrayX[iRow];
151#endif
152                    }
153                    scale2 = 1.0 / sqrt(scale2);
154                    //valueIncrease[i] = scale2;
155                    if (sequenceIncrease >= 0) {
156                         double djValue = dj_[sequenceIncrease];
157                         if (fabs(djValue) > 10.0 * dualTolerance_) {
158                              // we are going to use for cutoff so be exact
159                              costIncrease = fabs(djValue / alphaIncrease);
160                              /* Not sure this is good idea as I don't think correct e.g.
161                                 suppose a continuous variable has dj slightly greater. */
162                              if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) {
163                                   // can improve
164                                   double movement = (columnScale_ == NULL) ? 1.0 :
165                                                     rhsScale_ * inverseColumnScale_[sequenceIncrease];
166                                   costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
167                              }
168                         } else {
169                              costIncrease = 0.0;
170                         }
171                    }
172                    if (sequenceDecrease >= 0) {
173                         double djValue = dj_[sequenceDecrease];
174                         if (fabs(djValue) > 10.0 * dualTolerance_) {
175                              // we are going to use for cutoff so be exact
176                              costDecrease = fabs(djValue / alphaDecrease);
177                              if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
178                                   // can improve
179                                   double movement = (columnScale_ == NULL) ? 1.0 :
180                                                     rhsScale_ * inverseColumnScale_[sequenceDecrease];
181                                   costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
182                              }
183                         } else {
184                              costDecrease = 0.0;
185                         }
186                    }
187                    costIncrease *= scale2;
188                    costDecrease *= scale2;
189               }
190          }
191          break;
192          case isFixed:
193               break;
194          case isFree:
195          case superBasic:
196               costIncrease = 0.0;
197               costDecrease = 0.0;
198               sequenceIncrease = iSequence;
199               sequenceDecrease = iSequence;
200               break;
201          case atUpperBound:
202               costIncrease = CoinMax(0.0, -dj_[iSequence]);
203               sequenceIncrease = iSequence;
204               if (valueIncrease)
205                    valueIncrease[i] = primalRanging1(iSequence, iSequence);
206               break;
207          case atLowerBound:
208               costDecrease = CoinMax(0.0, dj_[iSequence]);
209               sequenceDecrease = iSequence;
210               if (valueIncrease)
211                    valueDecrease[i] = primalRanging1(iSequence, iSequence);
212               break;
213          }
214          double scaleFactor;
215          if (rowScale_) {
216               if (iSequence < numberColumns_)
217                    scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
218               else
219                    scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
220          } else {
221               scaleFactor = 1.0 / objectiveScale_;
222          }
223          if (costIncrease < 1.0e30)
224               costIncrease *= scaleFactor;
225          if (costDecrease < 1.0e30)
226               costDecrease *= scaleFactor;
227          if (optimizationDirection_ == 1.0) {
228               costIncreased[i] = costIncrease;
229               sequenceIncreased[i] = sequenceIncrease;
230               costDecreased[i] = costDecrease;
231               sequenceDecreased[i] = sequenceDecrease;
232          } else if (optimizationDirection_ == -1.0) {
233               costIncreased[i] = costDecrease;
234               sequenceIncreased[i] = sequenceDecrease;
235               costDecreased[i] = costIncrease;
236               sequenceDecreased[i] = sequenceIncrease;
237               if (valueIncrease) {
238                    double temp = valueIncrease[i];
239                    valueIncrease[i] = valueDecrease[i];
240                    valueDecrease[i] = temp;
241               }
242          } else if (optimizationDirection_ == 0.0) {
243               // !!!!!! ???
244               costIncreased[i] = COIN_DBL_MAX;
245               sequenceIncreased[i] = -1;
246               costDecreased[i] = COIN_DBL_MAX;
247               sequenceDecreased[i] = -1;
248          } else {
249               abort();
250          }
251     }
252     rowArray_[0]->clear();
253     //rowArray_[1]->clear();
254     //columnArray_[1]->clear();
255     columnArray_[0]->clear();
256     //rowArray_[3]->clear();
257     if (!optimizationDirection_)
258          printf("*** ????? Ranging with zero optimization costs\n");
259}
260/*
261   Row array has row part of pivot row
262   Column array has column part.
263   This is used in dual ranging
264*/
265void
266ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
267                                 CoinIndexedVector * columnArray,
268                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
269                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
270{
271     double acceptablePivot = 1.0e-9;
272     double * work;
273     int number;
274     int * which;
275     int iSection;
276
277     double thetaDown = 1.0e31;
278     double thetaUp = 1.0e31;
279     int sequenceDown = -1;
280     int sequenceUp = -1;
281     double alphaDown = 0.0;
282     double alphaUp = 0.0;
283
284     int addSequence;
285
286     for (iSection = 0; iSection < 2; iSection++) {
287
288          int i;
289          if (!iSection) {
290               work = rowArray->denseVector();
291               number = rowArray->getNumElements();
292               which = rowArray->getIndices();
293               addSequence = numberColumns_;
294          } else {
295               work = columnArray->denseVector();
296               number = columnArray->getNumElements();
297               which = columnArray->getIndices();
298               addSequence = 0;
299          }
300
301          for (i = 0; i < number; i++) {
302               int iSequence = which[i];
303               int iSequence2 = iSequence + addSequence;
304#ifndef COIN_FAC_NEW
305               double alpha = work[i];
306#else
307               double alpha = !addSequence ? work[i] : work[iSequence];
308#endif
309               if (fabs(alpha) < acceptablePivot)
310                    continue;
311               double oldValue = dj_[iSequence2];
312
313               switch(getStatus(iSequence2)) {
314
315               case basic:
316                    break;
317               case ClpSimplex::isFixed:
318                    break;
319               case isFree:
320               case superBasic:
321                    // treat dj as if zero
322                    thetaDown = 0.0;
323                    thetaUp = 0.0;
324                    sequenceDown = iSequence2;
325                    sequenceUp = iSequence2;
326                    break;
327               case atUpperBound:
328                    if (alpha > 0.0) {
329                         // test up
330                         if (oldValue + thetaUp * alpha > dualTolerance_) {
331                              thetaUp = (dualTolerance_ - oldValue) / alpha;
332                              sequenceUp = iSequence2;
333                              alphaUp = alpha;
334                         }
335                    } else {
336                         // test down
337                         if (oldValue - thetaDown * alpha > dualTolerance_) {
338                              thetaDown = -(dualTolerance_ - oldValue) / alpha;
339                              sequenceDown = iSequence2;
340                              alphaDown = alpha;
341                         }
342                    }
343                    break;
344               case atLowerBound:
345                    if (alpha < 0.0) {
346                         // test up
347                         if (oldValue + thetaUp * alpha < - dualTolerance_) {
348                              thetaUp = -(dualTolerance_ + oldValue) / alpha;
349                              sequenceUp = iSequence2;
350                              alphaUp = alpha;
351                         }
352                    } else {
353                         // test down
354                         if (oldValue - thetaDown * alpha < -dualTolerance_) {
355                              thetaDown = (dualTolerance_ + oldValue) / alpha;
356                              sequenceDown = iSequence2;
357                              alphaDown = alpha;
358                         }
359                    }
360                    break;
361               }
362          }
363     }
364     if (sequenceUp >= 0) {
365          costIncrease = thetaUp;
366          sequenceIncrease = sequenceUp;
367          alphaIncrease = alphaUp;
368     }
369     if (sequenceDown >= 0) {
370          costDecrease = thetaDown;
371          sequenceDecrease = sequenceDown;
372          alphaDecrease = alphaDown;
373     }
374}
375/** Primal ranging.
376    This computes increase/decrease in value for each given variable and corresponding
377    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
378    and numberColumns.. for artificials/slacks.
379    For basic variables the sequence number will be that of the basic variables.
380
381    Up to user to provide correct length arrays.
382
383    When here - guaranteed optimal
384*/
385void
386ClpSimplexOther::primalRanging(int numberCheck, const int * which,
387                               double * valueIncreased, int * sequenceIncreased,
388                               double * valueDecreased, int * sequenceDecreased)
389{
390     rowArray_[0]->clear();
391     rowArray_[1]->clear();
392     lowerIn_ = -COIN_DBL_MAX;
393     upperIn_ = COIN_DBL_MAX;
394     valueIn_ = 0.0;
395     for ( int i = 0; i < numberCheck; i++) {
396          int iSequence = which[i];
397          double valueIncrease = COIN_DBL_MAX;
398          double valueDecrease = COIN_DBL_MAX;
399          int sequenceIncrease = -1;
400          int sequenceDecrease = -1;
401
402          switch(getStatus(iSequence)) {
403
404          case basic:
405          case isFree:
406          case superBasic:
407               // Easy
408               valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
409               valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
410               sequenceDecrease = iSequence;
411               sequenceIncrease = iSequence;
412               break;
413          case isFixed:
414          case atUpperBound:
415          case atLowerBound: {
416               // Non trivial
417               // Other bound is ignored
418#ifndef COIN_FAC_NEW
419               unpackPacked(rowArray_[1], iSequence);
420#else
421               unpack(rowArray_[1], iSequence);
422#endif
423               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
424               // Get extra rows
425               matrix_->extendUpdated(this, rowArray_[1], 0);
426               // do ratio test
427               checkPrimalRatios(rowArray_[1], 1);
428               if (pivotRow_ >= 0) {
429                    valueIncrease = theta_;
430                    sequenceIncrease = pivotVariable_[pivotRow_];
431               }
432               checkPrimalRatios(rowArray_[1], -1);
433               if (pivotRow_ >= 0) {
434                    valueDecrease = theta_;
435                    sequenceDecrease = pivotVariable_[pivotRow_];
436               }
437               rowArray_[1]->clear();
438          }
439          break;
440          }
441          double scaleFactor;
442          if (rowScale_) {
443               if (iSequence < numberColumns_)
444                    scaleFactor = columnScale_[iSequence] / rhsScale_;
445               else
446                    scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
447          } else {
448               scaleFactor = 1.0 / rhsScale_;
449          }
450          if (valueIncrease < 1.0e30)
451               valueIncrease *= scaleFactor;
452          else
453               valueIncrease = COIN_DBL_MAX;
454          if (valueDecrease < 1.0e30)
455               valueDecrease *= scaleFactor;
456          else
457               valueDecrease = COIN_DBL_MAX;
458          valueIncreased[i] = valueIncrease;
459          sequenceIncreased[i] = sequenceIncrease;
460          valueDecreased[i] = valueDecrease;
461          sequenceDecreased[i] = sequenceDecrease;
462     }
463}
464// Returns new value of whichOther when whichIn enters basis
465double
466ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
467{
468     rowArray_[0]->clear();
469     rowArray_[1]->clear();
470     int iSequence = whichIn;
471     double newValue = solution_[whichOther];
472     double alphaOther = 0.0;
473     Status status = getStatus(iSequence);
474     assert (status == atLowerBound || status == atUpperBound);
475     int wayIn = (status == atLowerBound) ? 1 : -1;
476
477     switch(getStatus(iSequence)) {
478
479     case basic:
480     case isFree:
481     case superBasic:
482          assert (whichIn == whichOther);
483          // Easy
484          newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
485          break;
486     case isFixed:
487     case atUpperBound:
488     case atLowerBound:
489          // Non trivial
490     {
491          // Other bound is ignored
492#ifndef COIN_FAC_NEW
493          unpackPacked(rowArray_[1], iSequence);
494#else
495          unpack(rowArray_[1], iSequence);
496#endif
497          factorization_->updateColumn(rowArray_[2], rowArray_[1]);
498          // Get extra rows
499          matrix_->extendUpdated(this, rowArray_[1], 0);
500          // do ratio test
501          double acceptablePivot = 1.0e-7;
502          double * work = rowArray_[1]->denseVector();
503          int number = rowArray_[1]->getNumElements();
504          int * which = rowArray_[1]->getIndices();
505
506          // we may need to swap sign
507          double way = wayIn;
508          double theta = 1.0e30;
509          for (int iIndex = 0; iIndex < number; iIndex++) {
510
511               int iRow = which[iIndex];
512#ifndef COIN_FAC_NEW
513               double alpha = work[iIndex] * way;
514#else
515               double alpha = work[iRow] * way;
516#endif
517               int iPivot = pivotVariable_[iRow];
518               if (iPivot == whichOther) {
519                    alphaOther = alpha;
520                    continue;
521               }
522               double oldValue = solution_[iPivot];
523               if (fabs(alpha) > acceptablePivot) {
524                    if (alpha > 0.0) {
525                         // basic variable going towards lower bound
526                         double bound = lower_[iPivot];
527                         oldValue -= bound;
528                         if (oldValue - theta * alpha < 0.0) {
529                              theta = CoinMax(0.0, oldValue / alpha);
530                         }
531                    } else {
532                         // basic variable going towards upper bound
533                         double bound = upper_[iPivot];
534                         oldValue = oldValue - bound;
535                         if (oldValue - theta * alpha > 0.0) {
536                              theta = CoinMax(0.0, oldValue / alpha);
537                         }
538                    }
539               }
540          }
541          if (whichIn != whichOther) {
542               if (theta < 1.0e30)
543                    newValue -= theta * alphaOther;
544               else
545                    newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
546          } else {
547               newValue += theta * wayIn;
548          }
549     }
550     rowArray_[1]->clear();
551     break;
552     }
553     double scaleFactor;
554     if (rowScale_) {
555          if (whichOther < numberColumns_)
556               scaleFactor = columnScale_[whichOther] / rhsScale_;
557          else
558               scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
559     } else {
560          scaleFactor = 1.0 / rhsScale_;
561     }
562     if (newValue < 1.0e29)
563          if (newValue > -1.0e29)
564               newValue *= scaleFactor;
565          else
566               newValue = -COIN_DBL_MAX;
567     else
568          newValue = COIN_DBL_MAX;
569     return newValue;
570}
571/*
572   Row array has pivot column
573   This is used in primal ranging
574*/
575void
576ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
577                                   int direction)
578{
579     // sequence stays as row number until end
580     pivotRow_ = -1;
581     double acceptablePivot = 1.0e-7;
582     double * work = rowArray->denseVector();
583     int number = rowArray->getNumElements();
584     int * which = rowArray->getIndices();
585
586     // we need to swap sign if going down
587     double way = direction;
588     theta_ = 1.0e30;
589     for (int iIndex = 0; iIndex < number; iIndex++) {
590
591          int iRow = which[iIndex];
592#ifndef COIN_FAC_NEW
593          double alpha = work[iIndex] * way;
594#else
595          double alpha = work[iRow] * way;
596#endif
597          int iPivot = pivotVariable_[iRow];
598          double oldValue = solution_[iPivot];
599          if (fabs(alpha) > acceptablePivot) {
600               if (alpha > 0.0) {
601                    // basic variable going towards lower bound
602                    double bound = lower_[iPivot];
603                    oldValue -= bound;
604                    if (oldValue - theta_ * alpha < 0.0) {
605                         pivotRow_ = iRow;
606                         theta_ = CoinMax(0.0, oldValue / alpha);
607                    }
608               } else {
609                    // basic variable going towards upper bound
610                    double bound = upper_[iPivot];
611                    oldValue = oldValue - bound;
612                    if (oldValue - theta_ * alpha > 0.0) {
613                         pivotRow_ = iRow;
614                         theta_ = CoinMax(0.0, oldValue / alpha);
615                    }
616               }
617          }
618     }
619}
620/* Write the basis in MPS format to the specified file.
621   If writeValues true writes values of structurals
622   (and adds VALUES to end of NAME card)
623
624   Row and column names may be null.
625   formatType is
626   <ul>
627   <li> 0 - normal
628   <li> 1 - extra accuracy
629   <li> 2 - IEEE hex (later)
630   </ul>
631
632   Returns non-zero on I/O error
633
634   This is based on code contributed by Thorsten Koch
635*/
636int
637ClpSimplexOther::writeBasis(const char *filename,
638                            bool writeValues,
639                            int formatType) const
640{
641     formatType = CoinMax(0, formatType);
642     formatType = CoinMin(2, formatType);
643     if (!writeValues)
644          formatType = 0;
645     // See if INTEL if IEEE
646     if (formatType == 2) {
647          // test intel here and add 1 if not intel
648          double value = 1.0;
649          char x[8];
650          memcpy(x, &value, 8);
651          if (x[0] == 63) {
652               formatType ++; // not intel
653          } else {
654               assert (x[0] == 0);
655          }
656     }
657
658     char number[20];
659     FILE * fp = fopen(filename, "w");
660     if (!fp)
661          return -1;
662
663     // NAME card
664
665     if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
666          fprintf(fp, "NAME          BLANK      ");
667     } else {
668          fprintf(fp, "NAME          %s       ", strParam_[ClpProbName].c_str());
669     }
670     if (formatType >= 2)
671          fprintf(fp, "FREEIEEE");
672     else if (writeValues)
673          fprintf(fp, "VALUES");
674     // finish off name
675     fprintf(fp, "\n");
676     int iRow = 0;
677     for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
678          bool printit = false;
679          if( getColumnStatus(iColumn) == ClpSimplex::basic) {
680               printit = true;
681               // Find non basic row
682               for(; iRow < numberRows_; iRow++) {
683                    if (getRowStatus(iRow) != ClpSimplex::basic)
684                         break;
685               }
686               if (lengthNames_) {
687                    if (iRow != numberRows_) {
688                         fprintf(fp, " %s %-8s       %s",
689                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
690                                 columnNames_[iColumn].c_str(),
691                                 rowNames_[iRow].c_str());
692                         iRow++;
693                    } else {
694                         // Allow for too many basics!
695                         fprintf(fp, " BS %-8s       ",
696                                 columnNames_[iColumn].c_str());
697                         // Dummy row name if values
698                         if (writeValues)
699                              fprintf(fp, "      _dummy_");
700                    }
701               } else {
702                    // no names
703                    if (iRow != numberRows_) {
704                         fprintf(fp, " %s C%7.7d     R%7.7d",
705                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
706                                 iColumn, iRow);
707                         iRow++;
708                    } else {
709                         // Allow for too many basics!
710                         fprintf(fp, " BS C%7.7d", iColumn);
711                         // Dummy row name if values
712                         if (writeValues)
713                              fprintf(fp, "      _dummy_");
714                    }
715               }
716          } else  {
717               if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
718                    printit = true;
719                    if (lengthNames_)
720                         fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
721                    else
722                         fprintf(fp, " UL C%7.7d", iColumn);
723                    // Dummy row name if values
724                    if (writeValues)
725                         fprintf(fp, "      _dummy_");
726               }
727          }
728          if (printit && writeValues) {
729               // add value
730               CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
731               fprintf(fp, "     %s", number);
732          }
733          if (printit)
734               fprintf(fp, "\n");
735     }
736     fprintf(fp, "ENDATA\n");
737     fclose(fp);
738     return 0;
739}
740// Read a basis from the given filename
741int
742ClpSimplexOther::readBasis(const char *fileName)
743{
744     int status = 0;
745     bool canOpen = false;
746     if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
747          // stdin
748          canOpen = true;
749     } else {
750          FILE *fp = fopen(fileName, "r");
751          if (fp) {
752               // can open - lets go for it
753               fclose(fp);
754               canOpen = true;
755          } else {
756               handler_->message(CLP_UNABLE_OPEN, messages_)
757                         << fileName << CoinMessageEol;
758               return -1;
759          }
760     }
761     CoinMpsIO m;
762     m.passInMessageHandler(handler_);
763     *m.messagesPointer() = coinMessages();
764     bool savePrefix = m.messageHandler()->prefix();
765     m.messageHandler()->setPrefix(handler_->prefix());
766     status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
767                          status_,
768                          columnNames_, numberColumns_,
769                          rowNames_, numberRows_);
770     m.messageHandler()->setPrefix(savePrefix);
771     if (status >= 0) {
772          if (!status) {
773               // set values
774               int iColumn, iRow;
775               for (iRow = 0; iRow < numberRows_; iRow++) {
776                    if (getRowStatus(iRow) == atLowerBound)
777                         rowActivity_[iRow] = rowLower_[iRow];
778                    else if (getRowStatus(iRow) == atUpperBound)
779                         rowActivity_[iRow] = rowUpper_[iRow];
780               }
781               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
782                    if (getColumnStatus(iColumn) == atLowerBound)
783                         columnActivity_[iColumn] = columnLower_[iColumn];
784                    else if (getColumnStatus(iColumn) == atUpperBound)
785                         columnActivity_[iColumn] = columnUpper_[iColumn];
786               }
787          } else {
788               memset(rowActivity_, 0, numberRows_ * sizeof(double));
789               matrix_->times(-1.0, columnActivity_, rowActivity_);
790          }
791     } else {
792          // errors
793          handler_->message(CLP_IMPORT_ERRORS, messages_)
794                    << status << fileName << CoinMessageEol;
795     }
796     return status;
797}
798/* Creates dual of a problem if looks plausible
799   (defaults will always create model)
800   fractionRowRanges is fraction of rows allowed to have ranges
801   fractionColumnRanges is fraction of columns allowed to have ranges
802*/
803ClpSimplex *
804ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
805{
806     const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
807     bool changed = false;
808     int numberChanged = 0;
809     int iColumn;
810     // check if we need to change bounds to rows
811     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
812          if (columnUpper_[iColumn] < 1.0e20 &&
813                    columnLower_[iColumn] > -1.0e20) {
814               changed = true;
815               numberChanged++;
816          }
817     }
818     int iRow;
819     int numberExtraRows = 0;
820     if (numberChanged <= fractionColumnRanges * numberColumns_) {
821          for (iRow = 0; iRow < numberRows_; iRow++) {
822               if (rowLower_[iRow] > -1.0e20 &&
823                         rowUpper_[iRow] < 1.0e20) {
824                    if (rowUpper_[iRow] != rowLower_[iRow])
825                         numberExtraRows++;
826               }
827          }
828          if (numberExtraRows > fractionRowRanges * numberRows_)
829               return NULL;
830     } else {
831          return NULL;
832     }
833     if (changed) {
834          ClpSimplex * model3 = new ClpSimplex(*model2);
835          CoinBuild build;
836          double one = 1.0;
837          int numberColumns = model3->numberColumns();
838          const double * columnLower = model3->columnLower();
839          const double * columnUpper = model3->columnUpper();
840          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
841               if (columnUpper[iColumn] < 1.0e20 &&
842                         columnLower[iColumn] > -1.0e20) {
843                    if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
844                         double value = columnUpper[iColumn];
845                         model3->setColumnUpper(iColumn, COIN_DBL_MAX);
846                         build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
847                    } else {
848                         double value = columnLower[iColumn];
849                         model3->setColumnLower(iColumn, -COIN_DBL_MAX);
850                         build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
851                    }
852               }
853          }
854          model3->addRows(build);
855          model2 = model3;
856     }
857     int numberColumns = model2->numberColumns();
858     const double * columnLower = model2->columnLower();
859     const double * columnUpper = model2->columnUpper();
860     int numberRows = model2->numberRows();
861     double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
862     double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
863
864     const double * objective = model2->objective();
865     CoinPackedMatrix * matrix = model2->matrix();
866     // get transpose
867     CoinPackedMatrix rowCopy = *matrix;
868     const int * row = matrix->getIndices();
869     const int * columnLength = matrix->getVectorLengths();
870     const CoinBigIndex * columnStart = matrix->getVectorStarts();
871     const double * elementByColumn = matrix->getElements();
872     double objOffset = 0.0;
873     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
874          double offset = 0.0;
875          double objValue = optimizationDirection_ * objective[iColumn];
876          if (columnUpper[iColumn] > 1.0e20) {
877               if (columnLower[iColumn] > -1.0e20)
878                    offset = columnLower[iColumn];
879          } else if (columnLower[iColumn] < -1.0e20) {
880               offset = columnUpper[iColumn];
881          } else {
882               // taken care of before
883               abort();
884          }
885          if (offset) {
886               objOffset += offset * objValue;
887               for (CoinBigIndex j = columnStart[iColumn];
888                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
889                    int iRow = row[j];
890                    if (rowLower[iRow] > -1.0e20)
891                         rowLower[iRow] -= offset * elementByColumn[j];
892                    if (rowUpper[iRow] < 1.0e20)
893                         rowUpper[iRow] -= offset * elementByColumn[j];
894               }
895          }
896     }
897     int * which = new int[numberRows+numberExtraRows];
898     rowCopy.reverseOrdering();
899     rowCopy.transpose();
900     double * fromRowsLower = new double[numberRows+numberExtraRows];
901     double * fromRowsUpper = new double[numberRows+numberExtraRows];
902     double * newObjective = new double[numberRows+numberExtraRows];
903     double * fromColumnsLower = new double[numberColumns];
904     double * fromColumnsUpper = new double[numberColumns];
905     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
906          double objValue = optimizationDirection_ * objective[iColumn];
907          // Offset is already in
908          if (columnUpper[iColumn] > 1.0e20) {
909               if (columnLower[iColumn] > -1.0e20) {
910                    fromColumnsLower[iColumn] = -COIN_DBL_MAX;
911                    fromColumnsUpper[iColumn] = objValue;
912               } else {
913                    // free
914                    fromColumnsLower[iColumn] = objValue;
915                    fromColumnsUpper[iColumn] = objValue;
916               }
917          } else if (columnLower[iColumn] < -1.0e20) {
918               fromColumnsLower[iColumn] = objValue;
919               fromColumnsUpper[iColumn] = COIN_DBL_MAX;
920          } else {
921               abort();
922          }
923     }
924     int kRow = 0;
925     int kExtraRow = numberRows;
926     for (iRow = 0; iRow < numberRows; iRow++) {
927          if (rowLower[iRow] < -1.0e20) {
928               assert (rowUpper[iRow] < 1.0e20);
929               newObjective[kRow] = -rowUpper[iRow];
930               fromRowsLower[kRow] = -COIN_DBL_MAX;
931               fromRowsUpper[kRow] = 0.0;
932               which[kRow] = iRow;
933               kRow++;
934          } else if (rowUpper[iRow] > 1.0e20) {
935               newObjective[kRow] = -rowLower[iRow];
936               fromRowsLower[kRow] = 0.0;
937               fromRowsUpper[kRow] = COIN_DBL_MAX;
938               which[kRow] = iRow;
939               kRow++;
940          } else {
941               if (rowUpper[iRow] == rowLower[iRow]) {
942                    newObjective[kRow] = -rowLower[iRow];
943                    fromRowsLower[kRow] = -COIN_DBL_MAX;;
944                    fromRowsUpper[kRow] = COIN_DBL_MAX;
945                    which[kRow] = iRow;
946                    kRow++;
947               } else {
948                    // range
949                    newObjective[kRow] = -rowUpper[iRow];
950                    fromRowsLower[kRow] = -COIN_DBL_MAX;
951                    fromRowsUpper[kRow] = 0.0;
952                    which[kRow] = iRow;
953                    kRow++;
954                    newObjective[kExtraRow] = -rowLower[iRow];
955                    fromRowsLower[kExtraRow] = 0.0;
956                    fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
957                    which[kExtraRow] = iRow;
958                    kExtraRow++;
959               }
960          }
961     }
962     if (numberExtraRows) {
963          CoinPackedMatrix newCopy;
964          newCopy.setExtraGap(0.0);
965          newCopy.setExtraMajor(0.0);
966          newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
967          rowCopy = newCopy;
968     }
969     ClpSimplex * modelDual = new ClpSimplex();
970     modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
971                            fromColumnsLower, fromColumnsUpper);
972     modelDual->setObjectiveOffset(objOffset);
973     modelDual->setDualBound(model2->dualBound());
974     modelDual->setInfeasibilityCost(model2->infeasibilityCost());
975     modelDual->setDualTolerance(model2->dualTolerance());
976     modelDual->setPrimalTolerance(model2->primalTolerance());
977     modelDual->setPerturbation(model2->perturbation());
978     modelDual->setSpecialOptions(model2->specialOptions());
979     modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
980     delete [] fromRowsLower;
981     delete [] fromRowsUpper;
982     delete [] fromColumnsLower;
983     delete [] fromColumnsUpper;
984     delete [] newObjective;
985     delete [] which;
986     delete [] rowLower;
987     delete [] rowUpper;
988     if (changed)
989          delete model2;
990     modelDual->createStatus();
991     return modelDual;
992}
993// Restores solution from dualized problem
994int
995ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
996                                 bool checkAccuracy)
997{
998     int returnCode = 0;;
999     createStatus();
1000     // Number of rows in dual problem was original number of columns
1001     assert (numberColumns_ == dualProblem->numberRows());
1002     // If slack on d-row basic then column at bound otherwise column basic
1003     // If d-column basic then rhs tight
1004     int numberBasic = 0;
1005     int iRow, iColumn = 0;
1006     // Get number of extra rows from ranges
1007     int numberExtraRows = 0;
1008     for (iRow = 0; iRow < numberRows_; iRow++) {
1009          if (rowLower_[iRow] > -1.0e20 &&
1010                    rowUpper_[iRow] < 1.0e20) {
1011               if (rowUpper_[iRow] != rowLower_[iRow])
1012                    numberExtraRows++;
1013          }
1014     }
1015     const double * objective = this->objective();
1016     const double * dualDual = dualProblem->dualRowSolution();
1017     const double * dualDj = dualProblem->dualColumnSolution();
1018     const double * dualSol = dualProblem->primalColumnSolution();
1019     const double * dualActs = dualProblem->primalRowSolution();
1020#if 0
1021     ClpSimplex thisCopy = *this;
1022     thisCopy.dual(); // for testing
1023     const double * primalDual = thisCopy.dualRowSolution();
1024     const double * primalDj = thisCopy.dualColumnSolution();
1025     const double * primalSol = thisCopy.primalColumnSolution();
1026     const double * primalActs = thisCopy.primalRowSolution();
1027     char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
1028     printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
1029     for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
1030          printf("%d at %c primal %g dual %g\n",
1031                 iRow, ss[dualProblem->getRowStatus(iRow)],
1032                 dualActs[iRow], dualDual[iRow]);
1033     printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
1034     for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
1035          printf("%d at %c primal %g dual %g\n",
1036                 iColumn, ss[dualProblem->getColumnStatus(iColumn)],
1037                 dualSol[iColumn], dualDj[iColumn]);
1038     printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
1039     for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
1040          printf("%d at %c primal %g dual %g\n",
1041                 iRow, ss[thisCopy.getRowStatus(iRow)],
1042                 primalActs[iRow], primalDual[iRow]);
1043     printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
1044     for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
1045          printf("%d at %c primal %g dual %g\n",
1046                 iColumn, ss[thisCopy.getColumnStatus(iColumn)],
1047                 primalSol[iColumn], primalDj[iColumn]);
1048#endif
1049     // position at bound information
1050     int jColumn = numberRows_;
1051     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1052          double objValue = optimizationDirection_ * objective[iColumn];
1053          Status status = dualProblem->getRowStatus(iColumn);
1054          double otherValue = COIN_DBL_MAX;
1055          if (columnUpper_[iColumn] < 1.0e20 &&
1056                    columnLower_[iColumn] > -1.0e20) {
1057               if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1058                    otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1059               } else {
1060                    otherValue = columnLower_[iColumn] + dualDj[jColumn];
1061               }
1062               jColumn++;
1063          }
1064          if (status == basic) {
1065               // column is at bound
1066               if (otherValue == COIN_DBL_MAX) {
1067                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1068                    if (columnUpper_[iColumn] > 1.0e20) {
1069                         if (columnLower_[iColumn] > -1.0e20) {
1070                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1071                                   setColumnStatus(iColumn, atLowerBound);
1072                              else
1073                                   setColumnStatus(iColumn, isFixed);
1074                              columnActivity_[iColumn] = columnLower_[iColumn];
1075                         } else {
1076                              // free
1077                              setColumnStatus(iColumn, isFree);
1078                              columnActivity_[iColumn] = 0.0;
1079                         }
1080                    } else {
1081                         setColumnStatus(iColumn, atUpperBound);
1082                         columnActivity_[iColumn] = columnUpper_[iColumn];
1083                    }
1084               } else {
1085                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1086                    //printf("other dual sol %g\n",otherValue);
1087                    if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1088                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1089                              setColumnStatus(iColumn, atLowerBound);
1090                         else
1091                              setColumnStatus(iColumn, isFixed);
1092                         columnActivity_[iColumn] = columnLower_[iColumn];
1093                    } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1094                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1095                              setColumnStatus(iColumn, atUpperBound);
1096                         else
1097                              setColumnStatus(iColumn, isFixed);
1098                         columnActivity_[iColumn] = columnUpper_[iColumn];
1099                    } else {
1100                         abort();
1101                    }
1102               }
1103          } else {
1104               if (otherValue == COIN_DBL_MAX) {
1105                    // column basic
1106                    setColumnStatus(iColumn, basic);
1107                    numberBasic++;
1108                    if (columnLower_[iColumn] > -1.0e20) {
1109                         columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1110                    } else if (columnUpper_[iColumn] < 1.0e20) {
1111                         columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1112                    } else {
1113                         columnActivity_[iColumn] = -dualDual[iColumn];
1114                    }
1115                    reducedCost_[iColumn] = 0.0;
1116               } else {
1117                    // may be at other bound
1118                    //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1119                    if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1120                         // column basic
1121                         setColumnStatus(iColumn, basic);
1122                         numberBasic++;
1123                         //printf("Col %d otherV %g dualDual %g\n",iColumn,
1124                         // otherValue,dualDual[iColumn]);
1125                         columnActivity_[iColumn] = -dualDual[iColumn];
1126                         columnActivity_[iColumn] = otherValue;
1127                         reducedCost_[iColumn] = 0.0;
1128                    } else {
1129                         reducedCost_[iColumn] = objValue - dualActs[iColumn];
1130                         if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1131                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1132                                   setColumnStatus(iColumn, atLowerBound);
1133                              else
1134                                   setColumnStatus(iColumn, isFixed);
1135                              columnActivity_[iColumn] = columnLower_[iColumn];
1136                         } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1137                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1138                                   setColumnStatus(iColumn, atUpperBound);
1139                              else
1140                                   setColumnStatus(iColumn, isFixed);
1141                              columnActivity_[iColumn] = columnUpper_[iColumn];
1142                         } else {
1143                              abort();
1144                         }
1145                    }
1146               }
1147          }
1148     }
1149     // now rows
1150     int kExtraRow = jColumn;
1151     int numberRanges = 0;
1152     for (iRow = 0; iRow < numberRows_; iRow++) {
1153          Status status = dualProblem->getColumnStatus(iRow);
1154          if (status == basic) {
1155               // row is at bound
1156               dual_[iRow] = dualSol[iRow];;
1157          } else {
1158               // row basic
1159               setRowStatus(iRow, basic);
1160               numberBasic++;
1161               dual_[iRow] = 0.0;
1162          }
1163          if (rowLower_[iRow] < -1.0e20) {
1164               if (status == basic) {
1165                    rowActivity_[iRow] = rowUpper_[iRow];
1166                    setRowStatus(iRow, atUpperBound);
1167               } else {
1168                    assert (dualDj[iRow] < 1.0e-5);
1169                    rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1170               }
1171          } else if (rowUpper_[iRow] > 1.0e20) {
1172               if (status == basic) {
1173                    rowActivity_[iRow] = rowLower_[iRow];
1174                    setRowStatus(iRow, atLowerBound);
1175               } else {
1176                    rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1177                    assert (dualDj[iRow] > -1.0e-5);
1178               }
1179          } else {
1180               if (rowUpper_[iRow] == rowLower_[iRow]) {
1181                    rowActivity_[iRow] = rowLower_[iRow];
1182                    if (status == basic) {
1183                         setRowStatus(iRow, isFixed);
1184                    }
1185               } else {
1186                    // range
1187                    numberRanges++;
1188                    Status statusL = dualProblem->getColumnStatus(kExtraRow);
1189                    //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1190                    //     iRow,status,kExtraRow,statusL, dualSol[iRow],
1191                    //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1192                    if (status == basic) {
1193                         assert (statusL != basic);
1194                         rowActivity_[iRow] = rowUpper_[iRow];
1195                         setRowStatus(iRow, atUpperBound);
1196                    } else if (statusL == basic) {
1197                         numberBasic--; // already counted
1198                         rowActivity_[iRow] = rowLower_[iRow];
1199                         setRowStatus(iRow, atLowerBound);
1200                         dual_[iRow] = dualSol[kExtraRow];;
1201                    } else {
1202                         rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1203                         assert (dualDj[iRow] < 1.0e-5);
1204                         // row basic
1205                         //setRowStatus(iRow,basic);
1206                         //numberBasic++;
1207                         dual_[iRow] = 0.0;
1208                    }
1209                    kExtraRow++;
1210               }
1211          }
1212     }
1213     if (numberBasic != numberRows_) {
1214          printf("Bad basis - ranges - coding needed\n");
1215          assert (numberRanges);
1216          abort();
1217     }
1218     if (optimizationDirection_ < 0.0) {
1219          for (iRow = 0; iRow < numberRows_; iRow++) {
1220               dual_[iRow] = -dual_[iRow];
1221          }
1222     }
1223     // redo row activities
1224     memset(rowActivity_, 0, numberRows_ * sizeof(double));
1225     matrix_->times(1.0, columnActivity_, rowActivity_);
1226     // redo reduced costs
1227     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1228     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1229     checkSolutionInternal();
1230     if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
1231          returnCode = 1;
1232#ifdef CLP_INVESTIGATE
1233          printf("There are %d dual infeasibilities summing to %g ",
1234                 numberDualInfeasibilities_, sumDualInfeasibilities_);
1235          printf("and %d primal infeasibilities summing to %g\n",
1236                 numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
1237#endif
1238     }
1239     // Below will go to ..DEBUG later
1240#if 1 //ndef NDEBUG
1241     if (checkAccuracy) {
1242       // Check if correct
1243       double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1244       double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1245       double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1246       double * dual = CoinCopyOfArray(dual_, numberRows_);
1247       this->dual(); //primal();
1248       CoinRelFltEq eq(1.0e-5);
1249       for (iRow = 0; iRow < numberRows_; iRow++) {
1250         assert(eq(dual[iRow], dual_[iRow]));
1251       }
1252       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1253         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1254       }
1255       for (iRow = 0; iRow < numberRows_; iRow++) {
1256         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1257       }
1258       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1259         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1260       }
1261       delete [] columnActivity;
1262       delete [] rowActivity;
1263       delete [] reducedCost;
1264       delete [] dual;
1265     }
1266#endif
1267     return returnCode;
1268}
1269/* Does very cursory presolve.
1270   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1271*/
1272ClpSimplex *
1273ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1274                        int & nBound, bool moreBounds, bool tightenBounds)
1275{
1276     //#define CHECK_STATUS
1277#ifdef CHECK_STATUS
1278     {
1279          int n = 0;
1280          int i;
1281          for (i = 0; i < numberColumns_; i++)
1282               if (getColumnStatus(i) == ClpSimplex::basic)
1283                    n++;
1284          for (i = 0; i < numberRows_; i++)
1285               if (getRowStatus(i) == ClpSimplex::basic)
1286                    n++;
1287          assert (n == numberRows_);
1288     }
1289#endif
1290
1291     const double * element = matrix_->getElements();
1292     const int * row = matrix_->getIndices();
1293     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1294     const int * columnLength = matrix_->getVectorLengths();
1295
1296     CoinZeroN(rhs, numberRows_);
1297     int iColumn;
1298     int iRow;
1299     CoinZeroN(whichRow, numberRows_);
1300     int * backColumn = whichColumn + numberColumns_;
1301     int numberRows2 = 0;
1302     int numberColumns2 = 0;
1303     double offset = 0.0;
1304     const double * objective = this->objective();
1305     double * solution = columnActivity_;
1306     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1307          double lower = columnLower_[iColumn];
1308          double upper = columnUpper_[iColumn];
1309          if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
1310               backColumn[iColumn] = numberColumns2;
1311               whichColumn[numberColumns2++] = iColumn;
1312               for (CoinBigIndex j = columnStart[iColumn];
1313                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1314                    int iRow = row[j];
1315                    int n = whichRow[iRow];
1316                    if (n == 0 && element[j])
1317                         whichRow[iRow] = -iColumn - 1;
1318                    else if (n < 0)
1319                         whichRow[iRow] = 2;
1320               }
1321          } else {
1322               // fixed
1323               backColumn[iColumn] = -1;
1324               solution[iColumn] = upper;
1325               if (upper) {
1326                    offset += objective[iColumn] * upper;
1327                    for (CoinBigIndex j = columnStart[iColumn];
1328                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1329                         int iRow = row[j];
1330                         double value = element[j];
1331                         rhs[iRow] += upper * value;
1332                    }
1333               }
1334          }
1335     }
1336     int returnCode = 0;
1337     double tolerance = primalTolerance();
1338     nBound = 2 * numberRows_;
1339     for (iRow = 0; iRow < numberRows_; iRow++) {
1340          int n = whichRow[iRow];
1341          if (n > 0) {
1342               whichRow[numberRows2++] = iRow;
1343          } else if (n < 0) {
1344               //whichRow[numberRows2++]=iRow;
1345               //continue;
1346               // Can only do in certain circumstances as we don't know current value
1347               if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
1348                    // save row and column for bound
1349                    whichRow[--nBound] = iRow;
1350                    whichRow[nBound+numberRows_] = -n - 1;
1351               } else if (moreBounds) {
1352                    // save row and column for bound
1353                    whichRow[--nBound] = iRow;
1354                    whichRow[nBound+numberRows_] = -n - 1;
1355               } else {
1356                    whichRow[numberRows2++] = iRow;
1357               }
1358          } else {
1359               // empty
1360               double rhsValue = rhs[iRow];
1361               if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
1362                    returnCode = 1; // infeasible
1363               }
1364          }
1365     }
1366     ClpSimplex * small = NULL;
1367     if (!returnCode) {
1368       //printf("CRUNCH from (%d,%d) to (%d,%d)\n",
1369       //     numberRows_,numberColumns_,numberRows2,numberColumns2);
1370          small = new ClpSimplex(this, numberRows2, whichRow,
1371                                 numberColumns2, whichColumn, true, false);
1372#if 0
1373          ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1374          if (rowCopy) {
1375               assert(!small->rowCopy());
1376               small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1377                                    numberColumns2, whichColumn));
1378          }
1379#endif
1380          // Set some stuff
1381          small->setDualBound(dualBound_);
1382          small->setInfeasibilityCost(infeasibilityCost_);
1383          small->setSpecialOptions(specialOptions_);
1384          small->setPerturbation(perturbation_);
1385          small->defaultFactorizationFrequency();
1386          small->setAlphaAccuracy(alphaAccuracy_);
1387          // If no rows left then no tightening!
1388          if (!numberRows2 || !numberColumns2)
1389               tightenBounds = false;
1390
1391          int numberElements = getNumElements();
1392          int numberElements2 = small->getNumElements();
1393          small->setObjectiveOffset(objectiveOffset() - offset);
1394          handler_->message(CLP_CRUNCH_STATS, messages_)
1395                    << numberRows2 << -(numberRows_ - numberRows2)
1396                    << numberColumns2 << -(numberColumns_ - numberColumns2)
1397                    << numberElements2 << -(numberElements - numberElements2)
1398                    << CoinMessageEol;
1399          // And set objective value to match
1400          small->setObjectiveValue(this->objectiveValue());
1401          double * rowLower2 = small->rowLower();
1402          double * rowUpper2 = small->rowUpper();
1403          int jRow;
1404          for (jRow = 0; jRow < numberRows2; jRow++) {
1405               iRow = whichRow[jRow];
1406               if (rowLower2[jRow] > -1.0e20)
1407                    rowLower2[jRow] -= rhs[iRow];
1408               if (rowUpper2[jRow] < 1.0e20)
1409                    rowUpper2[jRow] -= rhs[iRow];
1410          }
1411          // and bounds
1412          double * columnLower2 = small->columnLower();
1413          double * columnUpper2 = small->columnUpper();
1414          const char * integerInformation = integerType_;
1415          for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1416               iRow = whichRow[jRow];
1417               iColumn = whichRow[jRow+numberRows_];
1418               double lowerRow = rowLower_[iRow];
1419               if (lowerRow > -1.0e20)
1420                    lowerRow -= rhs[iRow];
1421               double upperRow = rowUpper_[iRow];
1422               if (upperRow < 1.0e20)
1423                    upperRow -= rhs[iRow];
1424               int jColumn = backColumn[iColumn];
1425               double lower = columnLower2[jColumn];
1426               double upper = columnUpper2[jColumn];
1427               double value = 0.0;
1428               for (CoinBigIndex j = columnStart[iColumn];
1429                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1430                    if (iRow == row[j]) {
1431                         value = element[j];
1432                         break;
1433                    }
1434               }
1435               assert (value);
1436               // convert rowLower and Upper to implied bounds on column
1437               double newLower = -COIN_DBL_MAX;
1438               double newUpper = COIN_DBL_MAX;
1439               if (value > 0.0) {
1440                    if (lowerRow > -1.0e20)
1441                         newLower = lowerRow / value;
1442                    if (upperRow < 1.0e20)
1443                         newUpper = upperRow / value;
1444               } else {
1445                    if (upperRow < 1.0e20)
1446                         newLower = upperRow / value;
1447                    if (lowerRow > -1.0e20)
1448                         newUpper = lowerRow / value;
1449               }
1450               if (integerInformation && integerInformation[iColumn]) {
1451                    if (newLower - floor(newLower) < 10.0 * tolerance)
1452                         newLower = floor(newLower);
1453                    else
1454                         newLower = ceil(newLower);
1455                    if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1456                         newUpper = ceil(newUpper);
1457                    else
1458                         newUpper = floor(newUpper);
1459               }
1460               newLower = CoinMax(lower, newLower);
1461               newUpper = CoinMin(upper, newUpper);
1462               if (newLower > newUpper + tolerance) {
1463                    //printf("XXYY inf on bound\n");
1464                    returnCode = 1;
1465               }
1466               columnLower2[jColumn] = newLower;
1467               columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1468               if (getRowStatus(iRow) != ClpSimplex::basic) {
1469                    if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1470                         if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1471                              // can only get here if will be fixed
1472                              small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1473                         } else {
1474                              // solution is valid
1475                              if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1476                                        fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1477                                   small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1478                              else
1479                                   small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1480                         }
1481                    } else {
1482                         //printf("what now neither basic\n");
1483                    }
1484               }
1485          }
1486          if (returnCode) {
1487               delete small;
1488               small = NULL;
1489          } else if (tightenBounds && integerInformation) {
1490               // See if we can tighten any bounds
1491               // use rhs for upper and small duals for lower
1492               double * up = rhs;
1493               double * lo = small->dualRowSolution();
1494               const double * element = small->clpMatrix()->getElements();
1495               const int * row = small->clpMatrix()->getIndices();
1496               const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1497               //const int * columnLength = small->clpMatrix()->getVectorLengths();
1498               CoinZeroN(lo, numberRows2);
1499               CoinZeroN(up, numberRows2);
1500               for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1501                    double upper = columnUpper2[iColumn];
1502                    double lower = columnLower2[iColumn];
1503                    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1504                    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1505                         int iRow = row[j];
1506                         double value = element[j];
1507                         if (value > 0.0) {
1508                              if (upper < 1.0e20)
1509                                   up[iRow] += upper * value;
1510                              else
1511                                   up[iRow] = COIN_DBL_MAX;
1512                              if (lower > -1.0e20)
1513                                   lo[iRow] += lower * value;
1514                              else
1515                                   lo[iRow] = -COIN_DBL_MAX;
1516                         } else {
1517                              if (upper < 1.0e20)
1518                                   lo[iRow] += upper * value;
1519                              else
1520                                   lo[iRow] = -COIN_DBL_MAX;
1521                              if (lower > -1.0e20)
1522                                   up[iRow] += lower * value;
1523                              else
1524                                   up[iRow] = COIN_DBL_MAX;
1525                         }
1526                    }
1527               }
1528               double * rowLower2 = small->rowLower();
1529               double * rowUpper2 = small->rowUpper();
1530               bool feasible = true;
1531               // make safer
1532               for (int iRow = 0; iRow < numberRows2; iRow++) {
1533                    double lower = lo[iRow];
1534                    if (lower > rowUpper2[iRow] + tolerance) {
1535                         feasible = false;
1536                         break;
1537                    } else {
1538                         lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1539                    }
1540                    double upper = up[iRow];
1541                    if (upper < rowLower2[iRow] - tolerance) {
1542                         feasible = false;
1543                         break;
1544                    } else {
1545                         up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1546                    }
1547               }
1548               if (!feasible) {
1549                    delete small;
1550                    small = NULL;
1551               } else {
1552                    // and tighten
1553                    for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1554                         if (integerInformation[whichColumn[iColumn]]) {
1555                              double upper = columnUpper2[iColumn];
1556                              double lower = columnLower2[iColumn];
1557                              double newUpper = upper;
1558                              double newLower = lower;
1559                              double difference = upper - lower;
1560                              if (lower > -1000.0 && upper < 1000.0) {
1561                                   for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1562                                        int iRow = row[j];
1563                                        double value = element[j];
1564                                        if (value > 0.0) {
1565                                             double upWithOut = up[iRow] - value * difference;
1566                                             if (upWithOut < 0.0) {
1567                                                  newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1568                                             }
1569                                             double lowWithOut = lo[iRow] + value * difference;
1570                                             if (lowWithOut > 0.0) {
1571                                                  newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1572                                             }
1573                                        } else {
1574                                             double upWithOut = up[iRow] + value * difference;
1575                                             if (upWithOut < 0.0) {
1576                                                  newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1577                                             }
1578                                             double lowWithOut = lo[iRow] - value * difference;
1579                                             if (lowWithOut > 0.0) {
1580                                                  newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1581                                             }
1582                                        }
1583                                   }
1584                                   if (newLower > lower || newUpper < upper) {
1585                                        if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1586                                             newUpper = floor(newUpper);
1587                                        else
1588                                             newUpper = floor(newUpper + 0.5);
1589                                        if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1590                                             newLower = ceil(newLower);
1591                                        else
1592                                             newLower = ceil(newLower - 0.5);
1593                                        // change may be too small - check
1594                                        if (newLower > lower || newUpper < upper) {
1595                                             if (newUpper >= newLower) {
1596                                                  // Could also tighten in this
1597                                                  //printf("%d bounds %g %g tightened to %g %g\n",
1598                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1599                                                  //     newLower,newUpper);
1600#if 1
1601                                                  columnUpper2[iColumn] = newUpper;
1602                                                  columnLower2[iColumn] = newLower;
1603                                                  columnUpper_[whichColumn[iColumn]] = newUpper;
1604                                                  columnLower_[whichColumn[iColumn]] = newLower;
1605#endif
1606                                                  // and adjust bounds on rows
1607                                                  newUpper -= upper;
1608                                                  newLower -= lower;
1609                                                  for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1610                                                       int iRow = row[j];
1611                                                       double value = element[j];
1612                                                       if (value > 0.0) {
1613                                                            up[iRow] += newUpper * value;
1614                                                            lo[iRow] += newLower * value;
1615                                                       } else {
1616                                                            lo[iRow] += newUpper * value;
1617                                                            up[iRow] += newLower * value;
1618                                                       }
1619                                                  }
1620                                             } else {
1621                                                  // infeasible
1622                                                  //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1623                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1624                                                  //     newLower,newUpper);
1625#if 1
1626                                                  delete small;
1627                                                  small = NULL;
1628                                                  break;
1629#endif
1630                                             }
1631                                        }
1632                                   }
1633                              }
1634                         }
1635                    }
1636               }
1637          }
1638     }
1639#if 0
1640     if (small) {
1641          static int which = 0;
1642          which++;
1643          char xxxx[20];
1644          sprintf(xxxx, "bad%d.mps", which);
1645          small->writeMps(xxxx, 0, 1);
1646          sprintf(xxxx, "largebad%d.mps", which);
1647          writeMps(xxxx, 0, 1);
1648          printf("bad%d %x old size %d %d new %d %d\n", which, small,
1649                 numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
1650#if 0
1651          for (int i = 0; i < numberColumns_; i++)
1652               printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1653          for (int i = 0; i < numberRows_; i++)
1654               printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
1655#endif
1656     }
1657#endif
1658#ifdef CHECK_STATUS
1659     {
1660          int n = 0;
1661          int i;
1662          for (i = 0; i < small->numberColumns(); i++)
1663               if (small->getColumnStatus(i) == ClpSimplex::basic)
1664                    n++;
1665          for (i = 0; i < small->numberRows(); i++)
1666               if (small->getRowStatus(i) == ClpSimplex::basic)
1667                    n++;
1668          assert (n == small->numberRows());
1669     }
1670#endif
1671     return small;
1672}
1673/* After very cursory presolve.
1674   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1675*/
1676void
1677ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1678                             const int * whichRow,
1679                             const int * whichColumn, int nBound)
1680{
1681#ifndef NDEBUG
1682     for (int i = 0; i < small.numberRows(); i++)
1683          assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1684     for (int i = 0; i < small.numberColumns(); i++)
1685          assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
1686#endif
1687     getbackSolution(small, whichRow, whichColumn);
1688     // and deal with status for bounds
1689     const double * element = matrix_->getElements();
1690     const int * row = matrix_->getIndices();
1691     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1692     const int * columnLength = matrix_->getVectorLengths();
1693     double tolerance = primalTolerance();
1694     double djTolerance = dualTolerance();
1695     for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1696          int iRow = whichRow[jRow];
1697          int iColumn = whichRow[jRow+numberRows_];
1698          if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1699               double lower = columnLower_[iColumn];
1700               double upper = columnUpper_[iColumn];
1701               double value = columnActivity_[iColumn];
1702               double djValue = reducedCost_[iColumn];
1703               dual_[iRow] = 0.0;
1704               if (upper > lower) {
1705                    if (value < lower + tolerance && djValue > -djTolerance) {
1706                         setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1707                         setRowStatus(iRow, ClpSimplex::basic);
1708                    } else if (value > upper - tolerance && djValue < djTolerance) {
1709                         setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1710                         setRowStatus(iRow, ClpSimplex::basic);
1711                    } else {
1712                         // has to be basic
1713                         setColumnStatus(iColumn, ClpSimplex::basic);
1714                         reducedCost_[iColumn] = 0.0;
1715                         double value = 0.0;
1716                         for (CoinBigIndex j = columnStart[iColumn];
1717                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1718                              if (iRow == row[j]) {
1719                                   value = element[j];
1720                                   break;
1721                              }
1722                         }
1723                         dual_[iRow] = djValue / value;
1724                         if (rowUpper_[iRow] > rowLower_[iRow]) {
1725                              if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1726                                        fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1727                                   setRowStatus(iRow, ClpSimplex::atLowerBound);
1728                              else
1729                                   setRowStatus(iRow, ClpSimplex::atUpperBound);
1730                         } else {
1731                              setRowStatus(iRow, ClpSimplex::isFixed);
1732                         }
1733                    }
1734               } else {
1735                    // row can always be basic
1736                    setRowStatus(iRow, ClpSimplex::basic);
1737               }
1738          } else {
1739               // row can always be basic
1740               setRowStatus(iRow, ClpSimplex::basic);
1741          }
1742     }
1743     //#ifndef NDEBUG
1744#if 0
1745     if  (small.status() == 0) {
1746          int n = 0;
1747          int i;
1748          for (i = 0; i < numberColumns; i++)
1749               if (getColumnStatus(i) == ClpSimplex::basic)
1750                    n++;
1751          for (i = 0; i < numberRows; i++)
1752               if (getRowStatus(i) == ClpSimplex::basic)
1753                    n++;
1754          assert (n == numberRows);
1755     }
1756#endif
1757}
1758/* Tightens integer bounds - returns number tightened or -1 if infeasible
1759 */
1760int
1761ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1762{
1763     // See if we can tighten any bounds
1764     // use rhs for upper and small duals for lower
1765     double * up = rhsSpace;
1766     double * lo = dual_;
1767     const double * element = matrix_->getElements();
1768     const int * row = matrix_->getIndices();
1769     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1770     const int * columnLength = matrix_->getVectorLengths();
1771     CoinZeroN(lo, numberRows_);
1772     CoinZeroN(up, numberRows_);
1773     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1774          double upper = columnUpper_[iColumn];
1775          double lower = columnLower_[iColumn];
1776          //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1777          for (CoinBigIndex j = columnStart[iColumn];
1778                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1779               int iRow = row[j];
1780               double value = element[j];
1781               if (value > 0.0) {
1782                    if (upper < 1.0e20)
1783                         up[iRow] += upper * value;
1784                    else
1785                         up[iRow] = COIN_DBL_MAX;
1786                    if (lower > -1.0e20)
1787                         lo[iRow] += lower * value;
1788                    else
1789                         lo[iRow] = -COIN_DBL_MAX;
1790               } else {
1791                    if (upper < 1.0e20)
1792                         lo[iRow] += upper * value;
1793                    else
1794                         lo[iRow] = -COIN_DBL_MAX;
1795                    if (lower > -1.0e20)
1796                         up[iRow] += lower * value;
1797                    else
1798                         up[iRow] = COIN_DBL_MAX;
1799               }
1800          }
1801     }
1802     bool feasible = true;
1803     // make safer
1804     double tolerance = primalTolerance();
1805     for (int iRow = 0; iRow < numberRows_; iRow++) {
1806          double lower = lo[iRow];
1807          if (lower > rowUpper_[iRow] + tolerance) {
1808               feasible = false;
1809               break;
1810          } else {
1811               lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
1812          }
1813          double upper = up[iRow];
1814          if (upper < rowLower_[iRow] - tolerance) {
1815               feasible = false;
1816               break;
1817          } else {
1818               up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1819          }
1820     }
1821     int numberTightened = 0;
1822     if (!feasible) {
1823          return -1;
1824     } else if (integerType_) {
1825          // and tighten
1826          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1827               if (integerType_[iColumn]) {
1828                    double upper = columnUpper_[iColumn];
1829                    double lower = columnLower_[iColumn];
1830                    double newUpper = upper;
1831                    double newLower = lower;
1832                    double difference = upper - lower;
1833                    if (lower > -1000.0 && upper < 1000.0) {
1834                         for (CoinBigIndex j = columnStart[iColumn];
1835                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1836                              int iRow = row[j];
1837                              double value = element[j];
1838                              if (value > 0.0) {
1839                                   double upWithOut = up[iRow] - value * difference;
1840                                   if (upWithOut < 0.0) {
1841                                        newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1842                                   }
1843                                   double lowWithOut = lo[iRow] + value * difference;
1844                                   if (lowWithOut > 0.0) {
1845                                        newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1846                                   }
1847                              } else {
1848                                   double upWithOut = up[iRow] + value * difference;
1849                                   if (upWithOut < 0.0) {
1850                                        newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1851                                   }
1852                                   double lowWithOut = lo[iRow] - value * difference;
1853                                   if (lowWithOut > 0.0) {
1854                                        newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1855                                   }
1856                              }
1857                         }
1858                         if (newLower > lower || newUpper < upper) {
1859                              if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1860                                   newUpper = floor(newUpper);
1861                              else
1862                                   newUpper = floor(newUpper + 0.5);
1863                              if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1864                                   newLower = ceil(newLower);
1865                              else
1866                                   newLower = ceil(newLower - 0.5);
1867                              // change may be too small - check
1868                              if (newLower > lower || newUpper < upper) {
1869                                   if (newUpper >= newLower) {
1870                                        numberTightened++;
1871                                        //printf("%d bounds %g %g tightened to %g %g\n",
1872                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1873                                        //     newLower,newUpper);
1874                                        columnUpper_[iColumn] = newUpper;
1875                                        columnLower_[iColumn] = newLower;
1876                                        // and adjust bounds on rows
1877                                        newUpper -= upper;
1878                                        newLower -= lower;
1879                                        for (CoinBigIndex j = columnStart[iColumn];
1880                                                  j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1881                                             int iRow = row[j];
1882                                             double value = element[j];
1883                                             if (value > 0.0) {
1884                                                  up[iRow] += newUpper * value;
1885                                                  lo[iRow] += newLower * value;
1886                                             } else {
1887                                                  lo[iRow] += newUpper * value;
1888                                                  up[iRow] += newLower * value;
1889                                             }
1890                                        }
1891                                   } else {
1892                                        // infeasible
1893                                        //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1894                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1895                                        //     newLower,newUpper);
1896                                        return -1;
1897                                   }
1898                              }
1899                         }
1900                    }
1901               }
1902          }
1903     }
1904     return numberTightened;
1905}
1906/* Parametrics
1907   This is an initial slow version.
1908   The code uses current bounds + theta * change (if change array not NULL)
1909   and similarly for objective.
1910   It starts at startingTheta and returns ending theta in endingTheta.
1911   If reportIncrement 0.0 it will report on any movement
1912   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1913   If it can not reach input endingTheta return code will be 1 for infeasible,
1914   2 for unbounded, if error on ranges -1,  otherwise 0.
1915   Normal report is just theta and objective but
1916   if event handler exists it may do more
1917   On exit endingTheta is maximum reached (can be used for next startingTheta)
1918*/
1919int
1920ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
1921                             const double * lowerChangeBound, const double * upperChangeBound,
1922                             const double * lowerChangeRhs, const double * upperChangeRhs,
1923                             const double * changeObjective)
1924{
1925     bool needToDoSomething = true;
1926     bool canTryQuick = (reportIncrement) ? true : false;
1927     // Save copy of model
1928     ClpSimplex copyModel = *this;
1929     int savePerturbation = perturbation_;
1930     perturbation_ = 102; // switch off
1931     while (needToDoSomething) {
1932          needToDoSomething = false;
1933          algorithm_ = -1;
1934
1935          // save data
1936          ClpDataSave data = saveData();
1937          // Dantzig
1938          ClpDualRowPivot * savePivot = dualRowPivot_;
1939          dualRowPivot_ = new ClpDualRowDantzig();
1940          dualRowPivot_->setModel(this);
1941          int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1942          int iRow, iColumn;
1943          double * chgUpper = NULL;
1944          double * chgLower = NULL;
1945          double * chgObjective = NULL;
1946
1947
1948          if (!returnCode) {
1949               // Find theta when bounds will cross over and create arrays
1950               int numberTotal = numberRows_ + numberColumns_;
1951               chgLower = new double[numberTotal];
1952               memset(chgLower, 0, numberTotal * sizeof(double));
1953               chgUpper = new double[numberTotal];
1954               memset(chgUpper, 0, numberTotal * sizeof(double));
1955               chgObjective = new double[numberTotal];
1956               memset(chgObjective, 0, numberTotal * sizeof(double));
1957               assert (!rowScale_);
1958               double maxTheta = 1.0e50;
1959               if (lowerChangeRhs || upperChangeRhs) {
1960                    for (iRow = 0; iRow < numberRows_; iRow++) {
1961                         double lower = rowLower_[iRow];
1962                         double upper = rowUpper_[iRow];
1963                         if (lower > upper) {
1964                              maxTheta = -1.0;
1965                              break;
1966                         }
1967                         double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0;
1968                         double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0;
1969                         if (lower > -1.0e20 && upper < 1.0e20) {
1970                              if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1971                                   maxTheta = (upper - lower) / (lowerChange - upperChange);
1972                              }
1973                         }
1974                         if (lower > -1.0e20) {
1975                              lower_[numberColumns_+iRow] += startingTheta * lowerChange;
1976                              chgLower[numberColumns_+iRow] = lowerChange;
1977                         }
1978                         if (upper < 1.0e20) {
1979                              upper_[numberColumns_+iRow] += startingTheta * upperChange;
1980                              chgUpper[numberColumns_+iRow] = upperChange;
1981                         }
1982                    }
1983               }
1984               if (maxTheta > 0.0) {
1985                    if (lowerChangeBound || upperChangeBound) {
1986                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1987                              double lower = columnLower_[iColumn];
1988                              double upper = columnUpper_[iColumn];
1989                              if (lower > upper) {
1990                                   maxTheta = -1.0;
1991                                   break;
1992                              }
1993                              double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0;
1994                              double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0;
1995                              if (lower > -1.0e20 && upper < 1.0e20) {
1996                                   if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1997                                        maxTheta = (upper - lower) / (lowerChange - upperChange);
1998                                   }
1999                              }
2000                              if (lower > -1.0e20) {
2001                                   lower_[iColumn] += startingTheta * lowerChange;
2002                                   chgLower[iColumn] = lowerChange;
2003                              }
2004                              if (upper < 1.0e20) {
2005                                   upper_[iColumn] += startingTheta * upperChange;
2006                                   chgUpper[iColumn] = upperChange;
2007                              }
2008                         }
2009                    }
2010                    if (maxTheta == 1.0e50)
2011                         maxTheta = COIN_DBL_MAX;
2012               }
2013               if (maxTheta < 0.0) {
2014                    // bad ranges or initial
2015                    returnCode = -1;
2016               }
2017               if (maxTheta < endingTheta) {
2018                 char line[100];
2019                 sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
2020                         endingTheta,maxTheta);
2021                 handler_->message(CLP_GENERAL,messages_)
2022                   << line << CoinMessageEol;
2023                 endingTheta = maxTheta;
2024               }
2025               if (endingTheta < startingTheta) {
2026                    // bad initial
2027                    returnCode = -2;
2028               }
2029          }
2030          double saveEndingTheta = endingTheta;
2031          if (!returnCode) {
2032               if (changeObjective) {
2033                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2034                         chgObjective[iColumn] = changeObjective[iColumn];
2035                         cost_[iColumn] += startingTheta * changeObjective[iColumn];
2036                    }
2037               }
2038               double * saveDuals = NULL;
2039               reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
2040               assert (!problemStatus_);
2041               for (int i=0;i<numberRows_+numberColumns_;i++)
2042                 setFakeBound(i, noFake);
2043               // Now do parametrics
2044               handler_->message(CLP_PARAMETRICS_STATS, messages_)
2045                 << startingTheta << objectiveValue() << CoinMessageEol;
2046               while (!returnCode) {
2047                    //assert (reportIncrement);
2048                 parametricsData paramData;
2049                 paramData.startingTheta=startingTheta;
2050                 paramData.endingTheta=endingTheta;
2051                 paramData.maxTheta=COIN_DBL_MAX;
2052                 paramData.lowerChange = chgLower;
2053                 paramData.upperChange = chgUpper;
2054                    returnCode = parametricsLoop(paramData, reportIncrement,
2055                                                 chgLower, chgUpper, chgObjective, data,
2056                                                 canTryQuick);
2057                 startingTheta=paramData.startingTheta;
2058                 endingTheta=paramData.endingTheta;
2059                    if (!returnCode) {
2060                         //double change = endingTheta-startingTheta;
2061                         startingTheta = endingTheta;
2062                         endingTheta = saveEndingTheta;
2063                         //for (int i=0;i<numberTotal;i++) {
2064                         //lower_[i] += change*chgLower[i];
2065                         //upper_[i] += change*chgUpper[i];
2066                         //cost_[i] += change*chgObjective[i];
2067                         //}
2068                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
2069                           << startingTheta << objectiveValue() << CoinMessageEol;
2070                         if (startingTheta >= endingTheta)
2071                              break;
2072                    } else if (returnCode == -1) {
2073                         // trouble - do external solve
2074                         needToDoSomething = true;
2075                    } else if (problemStatus_==1) {
2076                      // can't move any further
2077                      if (!canTryQuick) {
2078                        handler_->message(CLP_PARAMETRICS_STATS, messages_)
2079                          << endingTheta << objectiveValue() << CoinMessageEol;
2080                        problemStatus_=0;
2081                      }
2082                    } else {
2083                         abort();
2084                    }
2085               }
2086          }
2087          reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2088
2089          delete dualRowPivot_;
2090          dualRowPivot_ = savePivot;
2091          // Restore any saved stuff
2092          restoreData(data);
2093          if (needToDoSomething) {
2094               double saveStartingTheta = startingTheta; // known to be feasible
2095               int cleanedUp = 1;
2096               while (cleanedUp) {
2097                    // tweak
2098                    if (cleanedUp == 1) {
2099                         if (!reportIncrement)
2100                              startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2101                         else
2102                              startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2103                    } else {
2104                         // restoring to go slowly
2105                         startingTheta = saveStartingTheta;
2106                    }
2107                    // only works if not scaled
2108                    int i;
2109                    const double * obj1 = objective();
2110                    double * obj2 = copyModel.objective();
2111                    const double * lower1 = columnLower_;
2112                    double * lower2 = copyModel.columnLower();
2113                    const double * upper1 = columnUpper_;
2114                    double * upper2 = copyModel.columnUpper();
2115                    for (i = 0; i < numberColumns_; i++) {
2116                         obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2117                         lower2[i] = lower1[i] + startingTheta * chgLower[i];
2118                         upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2119                    }
2120                    lower1 = rowLower_;
2121                    lower2 = copyModel.rowLower();
2122                    upper1 = rowUpper_;
2123                    upper2 = copyModel.rowUpper();
2124                    for (i = 0; i < numberRows_; i++) {
2125                         lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2126                         upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2127                    }
2128                    copyModel.dual();
2129                    if (copyModel.problemStatus()) {
2130                      char line[100];
2131                      sprintf(line,"Can not get to theta of %g\n", startingTheta);
2132                      handler_->message(CLP_GENERAL,messages_)
2133                        << line << CoinMessageEol;
2134                         canTryQuick = false; // do slowly to get exact amount
2135                         // back to last known good
2136                         if (cleanedUp == 1)
2137                              cleanedUp = 2;
2138                         else
2139                              abort();
2140                    } else {
2141                         // and move stuff back
2142                         int numberTotal = numberRows_ + numberColumns_;
2143                         CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2144                         CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2145                         CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2146                         cleanedUp = 0;
2147                    }
2148               }
2149          }
2150          delete [] chgLower;
2151          delete [] chgUpper;
2152          delete [] chgObjective;
2153     }
2154     perturbation_ = savePerturbation;
2155     char line[100];
2156     sprintf(line,"Ending theta %g\n", endingTheta);
2157     handler_->message(CLP_GENERAL,messages_)
2158       << line << CoinMessageEol;
2159     return problemStatus_;
2160}
2161/* Version of parametrics which reads from file
2162   See CbcClpParam.cpp for details of format
2163   Returns -2 if unable to open file */
2164int 
2165ClpSimplexOther::parametrics(const char * dataFile)
2166{
2167  int returnCode=-2;
2168  FILE *fp = fopen(dataFile, "r");
2169  char line[200];
2170  if (!fp) {
2171    handler_->message(CLP_UNABLE_OPEN, messages_)
2172      << dataFile << CoinMessageEol;
2173    return -2;
2174  }
2175
2176  if (!fgets(line, 200, fp)) {
2177    sprintf(line,"Empty parametrics file %s?",dataFile);
2178    handler_->message(CLP_GENERAL,messages_)
2179      << line << CoinMessageEol;
2180    fclose(fp);
2181    return -2;
2182  }
2183  char * pos = line;
2184  char * put = line;
2185  while (*pos >= ' ' && *pos != '\n') {
2186    if (*pos != ' ' && *pos != '\t') {
2187      *put = static_cast<char>(tolower(*pos));
2188      put++;
2189    }
2190    pos++;
2191  }
2192  *put = '\0';
2193  pos = line;
2194  double startTheta=0.0;
2195  double endTheta=0.0;
2196  double intervalTheta=COIN_DBL_MAX;
2197  int detail=0;
2198  bool good = true;
2199  while (good) {
2200    good=false;
2201    // check ROWS
2202    char * comma = strchr(pos, ',');
2203    if (!comma)
2204      break;
2205    *comma = '\0';
2206    if (strcmp(pos,"rows"))
2207      break;
2208    *comma = ',';
2209    pos = comma+1;
2210    // check lower theta
2211    comma = strchr(pos, ',');
2212    if (!comma)
2213      break;
2214    *comma = '\0';
2215    startTheta = atof(pos);
2216    *comma = ',';
2217    pos = comma+1;
2218    // check upper theta
2219    comma = strchr(pos, ',');
2220    good=true;
2221    if (comma)
2222      *comma = '\0';
2223    endTheta = atof(pos);
2224    if (comma) {
2225      *comma = ',';
2226      pos = comma+1;
2227      comma = strchr(pos, ',');
2228      if (comma)
2229        *comma = '\0';
2230      intervalTheta = atof(pos);
2231      if (comma) {
2232        *comma = ',';
2233        pos = comma+1;
2234        comma = strchr(pos, ',');
2235        if (comma)
2236          *comma = '\0';
2237        detail = atoi(pos);
2238        if (comma) 
2239        *comma = ',';
2240      }
2241    }
2242    break;
2243  }
2244  if (good) {
2245    if (startTheta<0.0||
2246        startTheta>endTheta||
2247        intervalTheta<0.0)
2248      good=false;
2249    if (detail<0||detail>1)
2250      good=false;
2251  }
2252  if (intervalTheta>=endTheta)
2253    intervalTheta=0.0;
2254  if (!good) {
2255    sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2256    handler_->message(CLP_GENERAL,messages_)
2257      << line << CoinMessageEol;
2258    fclose(fp);
2259    return -2;
2260  }
2261  if (!fgets(line, 200, fp)) {
2262    sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2263    handler_->message(CLP_GENERAL,messages_)
2264      << line << CoinMessageEol;
2265    fclose(fp);
2266    return -2;
2267  }
2268  double * lowerRowMove = NULL;
2269  double * upperRowMove = NULL;
2270  double * lowerColumnMove = NULL;
2271  double * upperColumnMove = NULL;
2272  double * objectiveMove = NULL;
2273  char saveLine[200];
2274  saveLine[0]='\0';
2275  std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2276  int gotRow[] = { -1, -1, -1, -1, -1};
2277  int orderRow[5];
2278  assert(sizeof(gotRow) == sizeof(orderRow));
2279  int nAcross = 0;
2280  pos = line;
2281  put = line;
2282  while (*pos >= ' ' && *pos != '\n') {
2283    if (*pos != ' ' && *pos != '\t') {
2284      *put = static_cast<char>(tolower(*pos));
2285      put++;
2286    }
2287    pos++;
2288  }
2289  *put = '\0';
2290  pos = line;
2291  int i;
2292  good = true;
2293  if (strncmp(line,"column",6)) {
2294    while (pos) {
2295      char * comma = strchr(pos, ',');
2296      if (comma)
2297        *comma = '\0';
2298      for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2299        if (headingsRow[i] == pos) {
2300          if (gotRow[i] < 0) {
2301            orderRow[nAcross] = i;
2302            gotRow[i] = nAcross++;
2303          } else {
2304            // duplicate
2305            good = false;
2306          }
2307          break;
2308        }
2309      }
2310      if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2311        good = false;
2312      if (comma) {
2313        *comma = ',';
2314        pos = comma + 1;
2315      } else {
2316        break;
2317      }
2318    }
2319    if (gotRow[0] < 0 && gotRow[1] < 0)
2320      good = false;
2321    if (gotRow[0] >= 0 && gotRow[1] >= 0)
2322      good = false;
2323    if (gotRow[0] >= 0 && !lengthNames())
2324      good = false;
2325    if (gotRow[4]<0) {
2326      if (gotRow[2] < 0 && gotRow[3] >= 0)
2327        good = false;
2328      else if (gotRow[3] < 0 && gotRow[2] >= 0)
2329        good = false;
2330    } else if (gotRow[2]>=0||gotRow[3]>=0) {
2331      good = false;
2332    }
2333    if (good) {
2334      char ** rowNames = new char * [numberRows_];
2335      int iRow;
2336      for (iRow = 0; iRow < numberRows_; iRow++) {
2337        rowNames[iRow] =
2338          CoinStrdup(rowName(iRow).c_str());
2339      }
2340      lowerRowMove = new double [numberRows_];
2341      memset(lowerRowMove,0,numberRows_*sizeof(double));
2342      upperRowMove = new double [numberRows_];
2343      memset(upperRowMove,0,numberRows_*sizeof(double));
2344      int nLine = 0;
2345      int nBadLine = 0;
2346      int nBadName = 0;
2347      bool goodLine=false;
2348      while (fgets(line, 200, fp)) {
2349        goodLine=true;
2350        if (!strncmp(line, "ENDATA", 6)||
2351            !strncmp(line, "COLUMN",6))
2352          break;
2353        goodLine=false;
2354        nLine++;
2355        iRow = -1;
2356        double upper = 0.0;
2357        double lower = 0.0;
2358        char * pos = line;
2359        char * put = line;
2360        while (*pos >= ' ' && *pos != '\n') {
2361          if (*pos != ' ' && *pos != '\t') {
2362            *put = *pos;
2363            put++;
2364          }
2365          pos++;
2366        }
2367        *put = '\0';
2368        pos = line;
2369        for (int i = 0; i < nAcross; i++) {
2370          char * comma = strchr(pos, ',');
2371          if (comma) {
2372            *comma = '\0';
2373          } else if (i < nAcross - 1) {
2374            nBadLine++;
2375            break;
2376          }
2377          switch (orderRow[i]) {
2378            // name
2379          case 0:
2380            // For large problems this could be slow
2381            for (iRow = 0; iRow < numberRows_; iRow++) {
2382              if (!strcmp(rowNames[iRow], pos))
2383                break;
2384            }
2385            if (iRow == numberRows_)
2386              iRow = -1;
2387            break;
2388            // number
2389          case 1:
2390            iRow = atoi(pos);
2391            if (iRow < 0 || iRow >= numberRows_)
2392              iRow = -1;
2393            break;
2394            // lower
2395          case 2:
2396            upper = atof(pos);
2397            break;
2398            // upper
2399          case 3:
2400            lower = atof(pos);
2401            break;
2402            // rhs
2403          case 4:
2404            lower = atof(pos);
2405            upper = lower;
2406            break;
2407          }
2408          if (comma) {
2409            *comma = ',';
2410            pos = comma + 1;
2411          }
2412        }
2413        if (iRow >= 0) {
2414          if (rowLower_[iRow]>-1.0e20)
2415            lowerRowMove[iRow] = lower;
2416          else
2417            lowerRowMove[iRow]=0.0;
2418          if (rowUpper_[iRow]<1.0e20)
2419            upperRowMove[iRow] = upper;
2420          else
2421            upperRowMove[iRow] = lower;
2422        } else {
2423          nBadName++;
2424          if(saveLine[0]=='\0')
2425            strcpy(saveLine,line);
2426        }
2427      }
2428      sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2429      handler_->message(CLP_GENERAL,messages_)
2430        << line << CoinMessageEol;
2431      if (nBadName) {
2432        sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2433        handler_->message(CLP_GENERAL,messages_)
2434          << line << CoinMessageEol;
2435        returnCode=-1;
2436        good=false;
2437      }
2438      for (iRow = 0; iRow < numberRows_; iRow++) {
2439        free(rowNames[iRow]);
2440      }
2441      delete [] rowNames;
2442    } else {
2443      sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2444      handler_->message(CLP_GENERAL,messages_)
2445        << line << CoinMessageEol;
2446      returnCode=-1;
2447      good=false;
2448    }
2449  }
2450  if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
2451    if (!fgets(line, 200, fp)) {
2452      sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2453      handler_->message(CLP_GENERAL,messages_)
2454        << line << CoinMessageEol;
2455      fclose(fp);
2456      return -2;
2457    }
2458    std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2459    saveLine[0]='\0';
2460    int gotColumn[] = { -1, -1, -1, -1, -1};
2461    int orderColumn[5];
2462    assert(sizeof(gotColumn) == sizeof(orderColumn));
2463    nAcross = 0;
2464    pos = line;
2465    put = line;
2466    while (*pos >= ' ' && *pos != '\n') {
2467      if (*pos != ' ' && *pos != '\t') {
2468        *put = static_cast<char>(tolower(*pos));
2469        put++;
2470      }
2471      pos++;
2472    }
2473    *put = '\0';
2474    pos = line;
2475    int i;
2476    if (strncmp(line,"endata",6)&&good) {
2477      while (pos) {
2478        char * comma = strchr(pos, ',');
2479        if (comma)
2480          *comma = '\0';
2481        for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2482          if (headingsColumn[i] == pos) {
2483            if (gotColumn[i] < 0) {
2484              orderColumn[nAcross] = i;
2485              gotColumn[i] = nAcross++;
2486            } else {
2487              // duplicate
2488              good = false;
2489            }
2490            break;
2491          }
2492        }
2493        if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2494          good = false;
2495        if (comma) {
2496          *comma = ',';
2497          pos = comma + 1;
2498        } else {
2499          break;
2500        }
2501      }
2502      if (gotColumn[0] < 0 && gotColumn[1] < 0)
2503        good = false;
2504      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2505        good = false;
2506      if (gotColumn[0] >= 0 && !lengthNames())
2507        good = false;
2508      if (good) {
2509        char ** columnNames = new char * [numberColumns_];
2510        int iColumn;
2511        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2512          columnNames[iColumn] =
2513            CoinStrdup(columnName(iColumn).c_str());
2514        }
2515        lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2516        memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2517        upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2518        memset(upperColumnMove,0,numberColumns_*sizeof(double));
2519        objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2520        memset(objectiveMove,0,numberColumns_*sizeof(double));
2521        int nLine = 0;
2522        int nBadLine = 0;
2523        int nBadName = 0;
2524        bool goodLine=false;
2525        while (fgets(line, 200, fp)) {
2526          goodLine=true;
2527          if (!strncmp(line, "ENDATA", 6))
2528            break;
2529          goodLine=false;
2530          nLine++;
2531          iColumn = -1;
2532          double upper = 0.0;
2533          double lower = 0.0;
2534          double obj =0.0;
2535          char * pos = line;
2536          char * put = line;
2537          while (*pos >= ' ' && *pos != '\n') {
2538            if (*pos != ' ' && *pos != '\t') {
2539              *put = *pos;
2540              put++;
2541            }
2542            pos++;
2543          }
2544          *put = '\0';
2545          pos = line;
2546          for (int i = 0; i < nAcross; i++) {
2547            char * comma = strchr(pos, ',');
2548            if (comma) {
2549              *comma = '\0';
2550            } else if (i < nAcross - 1) {
2551              nBadLine++;
2552              break;
2553            }
2554            switch (orderColumn[i]) {
2555              // name
2556            case 0:
2557              // For large problems this could be slow
2558              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2559                if (!strcmp(columnNames[iColumn], pos))
2560                  break;
2561              }
2562              if (iColumn == numberColumns_)
2563                iColumn = -1;
2564              break;
2565              // number
2566            case 1:
2567              iColumn = atoi(pos);
2568              if (iColumn < 0 || iColumn >= numberColumns_)
2569                iColumn = -1;
2570              break;
2571              // lower
2572            case 2:
2573              upper = atof(pos);
2574              break;
2575              // upper
2576            case 3:
2577              lower = atof(pos);
2578              break;
2579              // objective
2580            case 4:
2581              obj = atof(pos);
2582              upper = lower;
2583              break;
2584            }
2585            if (comma) {
2586              *comma = ',';
2587              pos = comma + 1;
2588            }
2589          }
2590          if (iColumn >= 0) {
2591            if (columnLower_[iColumn]>-1.0e20)
2592              lowerColumnMove[iColumn] = lower;
2593            else
2594              lowerColumnMove[iColumn]=0.0;
2595            if (columnUpper_[iColumn]<1.0e20)
2596              upperColumnMove[iColumn] = upper;
2597            else
2598              upperColumnMove[iColumn] = lower;
2599            objectiveMove[iColumn] = obj;
2600          } else {
2601            nBadName++;
2602            if(saveLine[0]=='\0')
2603              strcpy(saveLine,line);
2604          }
2605        }
2606        sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2607        handler_->message(CLP_GENERAL,messages_)
2608          << line << CoinMessageEol;
2609        if (nBadName) {
2610          sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2611          handler_->message(CLP_GENERAL,messages_)
2612            << line << CoinMessageEol;
2613          returnCode=-1;
2614          good=false;
2615        }
2616        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2617          free(columnNames[iColumn]);
2618        }
2619        delete [] columnNames;
2620      } else {
2621        sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2622        handler_->message(CLP_GENERAL,messages_)
2623          << line << CoinMessageEol;
2624        returnCode=-1;
2625        good=false;
2626      }
2627    }
2628  }
2629  returnCode=-1;
2630  if (good) {
2631    // clean arrays
2632    if (lowerRowMove) {
2633      bool empty=true;
2634      for (int i=0;i<numberRows_;i++) {
2635        if (lowerRowMove[i]) {
2636          empty=false;
2637        break;
2638        }
2639      }
2640      if (empty) {
2641        delete [] lowerRowMove;
2642        lowerRowMove=NULL;
2643      }
2644    }
2645    if (upperRowMove) {
2646      bool empty=true;
2647      for (int i=0;i<numberRows_;i++) {
2648        if (upperRowMove[i]) {
2649          empty=false;
2650        break;
2651        }
2652      }
2653      if (empty) {
2654        delete [] upperRowMove;
2655        upperRowMove=NULL;
2656      }
2657    }
2658    if (lowerColumnMove) {
2659      bool empty=true;
2660      for (int i=0;i<numberColumns_;i++) {
2661        if (lowerColumnMove[i]) {
2662          empty=false;
2663        break;
2664        }
2665      }
2666      if (empty) {
2667        delete [] lowerColumnMove;
2668        lowerColumnMove=NULL;
2669      }
2670    }
2671    if (upperColumnMove) {
2672      bool empty=true;
2673      for (int i=0;i<numberColumns_;i++) {
2674        if (upperColumnMove[i]) {
2675          empty=false;
2676        break;
2677        }
2678      }
2679      if (empty) {
2680        delete [] upperColumnMove;
2681        upperColumnMove=NULL;
2682      }
2683    }
2684    if (objectiveMove) {
2685      bool empty=true;
2686      for (int i=0;i<numberColumns_;i++) {
2687        if (objectiveMove[i]) {
2688          empty=false;
2689        break;
2690        }
2691      }
2692      if (empty) {
2693        delete [] objectiveMove;
2694        objectiveMove=NULL;
2695      }
2696    }
2697    int saveScaling = scalingFlag_;
2698    scalingFlag_ = 0;
2699    int saveLogLevel = handler_->logLevel();
2700    if (detail>0&&!intervalTheta)
2701      handler_->setLogLevel(3);
2702    else
2703      handler_->setLogLevel(1);
2704    returnCode = parametrics(startTheta,endTheta,intervalTheta,
2705                             lowerColumnMove,upperColumnMove,
2706                             lowerRowMove,upperRowMove,
2707                             objectiveMove);
2708    scalingFlag_ = saveScaling;
2709    handler_->setLogLevel(saveLogLevel);
2710  }
2711  delete [] lowerRowMove;
2712  delete [] upperRowMove;
2713  delete [] lowerColumnMove;
2714  delete [] upperColumnMove;
2715  delete [] objectiveMove;
2716  fclose(fp);
2717  return returnCode;
2718}
2719int
2720ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
2721                                 const double * lowerChange, const double * upperChange,
2722                                 const double * changeObjective, ClpDataSave & data,
2723                                 bool canTryQuick)
2724{
2725  double startingTheta = paramData.startingTheta;
2726  double & endingTheta = paramData.endingTheta;
2727     // stuff is already at starting
2728     // For this crude version just try and go to end
2729     double change = 0.0;
2730     if (reportIncrement && canTryQuick) {
2731          endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2732          change = endingTheta - startingTheta;
2733     }
2734     int numberTotal = numberRows_ + numberColumns_;
2735     int i;
2736     for ( i = 0; i < numberTotal; i++) {
2737          lower_[i] += change * lowerChange[i];
2738          upper_[i] += change * upperChange[i];
2739          switch(getStatus(i)) {
2740
2741          case basic:
2742          case isFree:
2743          case superBasic:
2744               break;
2745          case isFixed:
2746          case atUpperBound:
2747               solution_[i] = upper_[i];
2748               break;
2749          case atLowerBound:
2750               solution_[i] = lower_[i];
2751               break;
2752          }
2753          cost_[i] += change * changeObjective[i];
2754     }
2755     problemStatus_ = -1;
2756
2757     // This says whether to restore things etc
2758     // startup will have factorized so can skip
2759     int factorType = 0;
2760     // Start check for cycles
2761     progress_.startCheck();
2762     // Say change made on first iteration
2763     changeMade_ = 1;
2764     /*
2765       Status of problem:
2766       0 - optimal
2767       1 - infeasible
2768       2 - unbounded
2769       -1 - iterating
2770       -2 - factorization wanted
2771       -3 - redo checking without factorization
2772       -4 - looks infeasible
2773     */
2774     while (problemStatus_ < 0) {
2775          int iRow, iColumn;
2776          // clear
2777          for (iRow = 0; iRow < 4; iRow++) {
2778               rowArray_[iRow]->clear();
2779          }
2780
2781          for (iColumn = 0; iColumn < 2; iColumn++) {
2782               columnArray_[iColumn]->clear();
2783          }
2784
2785          // give matrix (and model costs and bounds a chance to be
2786          // refreshed (normally null)
2787          matrix_->refresh(this);
2788          // may factorize, checks if problem finished
2789          statusOfProblemInParametrics(factorType, data);
2790          // Say good factorization
2791          factorType = 1;
2792          if (data.sparseThreshold_) {
2793               // use default at present
2794               factorization_->sparseThreshold(0);
2795               factorization_->goSparse();
2796          }
2797
2798          // exit if victory declared
2799          if (problemStatus_ >= 0 && 
2800              (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
2801               break;
2802
2803          // test for maximum iterations
2804          if (hitMaximumIterations()) {
2805               problemStatus_ = 3;
2806               break;
2807          }
2808          // Check event
2809          {
2810               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2811               if (status >= 0) {
2812                    problemStatus_ = 5;
2813                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
2814                    break;
2815               }
2816          }
2817          // Do iterations
2818          problemStatus_=-1;
2819          if (canTryQuick) {
2820               double * saveDuals = NULL;
2821               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2822          } else {
2823               whileIterating(paramData, reportIncrement,
2824                              changeObjective);
2825               startingTheta = endingTheta;
2826          }
2827     }
2828     if (!problemStatus_) {
2829          theta_ = change + startingTheta;
2830          eventHandler_->event(ClpEventHandler::theta);
2831          return 0;
2832     } else if (problemStatus_ == 10) {
2833          return -1;
2834     } else {
2835          return problemStatus_;
2836     }
2837}
2838/* Parametrics
2839   The code uses current bounds + theta * change (if change array not NULL)
2840   It starts at startingTheta and returns ending theta in endingTheta.
2841   If it can not reach input endingTheta return code will be 1 for infeasible,
2842   2 for unbounded, if error on ranges -1,  otherwise 0.
2843   Event handler may do more
2844   On exit endingTheta is maximum reached (can be used for next startingTheta)
2845*/
2846int
2847ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
2848                             const double * lowerChangeBound, const double * upperChangeBound,
2849                             const double * lowerChangeRhs, const double * upperChangeRhs)
2850{
2851  int savePerturbation = perturbation_;
2852  perturbation_ = 102; // switch off
2853  algorithm_ = -1;
2854  // extra region
2855  int maximumPivots = factorization_->maximumPivots();
2856  int numberDense = factorization_->numberDense();
2857  int length = numberRows_ + numberDense + maximumPivots;
2858  assert (!rowArray_[4]);
2859  rowArray_[4]=new CoinIndexedVector(length);
2860  assert (!rowArray_[5]);
2861  rowArray_[5]=new CoinIndexedVector(length);
2862
2863  // save data
2864  ClpDataSave data = saveData();
2865  int numberTotal = numberRows_ + numberColumns_;
2866  int ratio = (2*sizeof(int))/sizeof(double);
2867  assert (ratio==1||ratio==2);
2868  // allow for unscaled - even if not needed
2869  int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1;
2870  /*
2871    Save information and modify
2872  */
2873  double * saveLower = new double [lengthArrays];
2874  double * saveUpper = new double [lengthArrays];
2875  double * lowerCopy = saveLower+2*numberTotal;
2876  double * upperCopy = saveUpper+2*numberTotal;
2877  double * lowerChange = saveLower+numberTotal;
2878  double * upperChange = saveUpper+numberTotal;
2879  double * lowerGap = saveLower+4*numberTotal;
2880  double * lowerCoefficient = lowerGap+numberRows_;
2881  double * upperGap = saveUpper+4*numberTotal;
2882  double * upperCoefficient = upperGap+numberRows_;
2883  int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal+2*numberRows_))+2;
2884  int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal+2*numberRows_))+2;
2885  int * lowerActive = lowerList+numberTotal+1;
2886  int * upperActive = upperList+numberTotal+1;
2887  // To mark as odd
2888  char * markDone = reinterpret_cast<char *>(lowerActive+numberTotal);
2889  //memset(markDone,0,numberTotal);
2890  int * backwardBasic = upperActive+numberTotal;
2891  parametricsData paramData;
2892  paramData.lowerChange = lowerChange;
2893  paramData.lowerList=lowerList;
2894  paramData.upperChange = upperChange;
2895  paramData.upperList=upperList;
2896  paramData.markDone=markDone;
2897  paramData.backwardBasic=backwardBasic;
2898  paramData.lowerActive = lowerActive;
2899  paramData.lowerGap = lowerGap;
2900  paramData.lowerCoefficient = lowerCoefficient;
2901  paramData.upperActive = upperActive;
2902  paramData.upperGap = upperGap;
2903  paramData.upperCoefficient = upperCoefficient;
2904  // Find theta when bounds will cross over and create arrays
2905  memset(lowerChange, 0, numberTotal * sizeof(double));
2906  memset(upperChange, 0, numberTotal * sizeof(double));
2907  if (lowerChangeBound)
2908    memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
2909  if (upperChangeBound)
2910    memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
2911  if (lowerChangeRhs)
2912    memcpy(lowerChange+numberColumns_,
2913           lowerChangeRhs,numberRows_*sizeof(double));
2914  if (upperChangeRhs)
2915    memcpy(upperChange+numberColumns_,
2916           upperChangeRhs,numberRows_*sizeof(double));
2917  int nLowerChange=0;
2918  int nUpperChange=0;
2919  for (int i=0;i<numberColumns_;i++) {
2920    if (lowerChange[i]) { 
2921      lowerList[nLowerChange++]=i;
2922    }
2923    if (upperChange[i]) { 
2924      upperList[nUpperChange++]=i;
2925    }
2926  }
2927  lowerList[-2]=nLowerChange;
2928  upperList[-2]=nUpperChange;
2929  for (int i=numberColumns_;i<numberTotal;i++) {
2930    if (lowerChange[i]) { 
2931      lowerList[nLowerChange++]=i;
2932    }
2933    if (upperChange[i]) { 
2934      upperList[nUpperChange++]=i;
2935    }
2936  }
2937  lowerList[-1]=nLowerChange;
2938  upperList[-1]=nUpperChange;
2939  memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
2940  memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
2941  memcpy(lowerCopy+numberColumns_,
2942         rowLower_,numberRows_*sizeof(double));
2943  memcpy(upperCopy+numberColumns_,
2944         rowUpper_,numberRows_*sizeof(double));
2945  {
2946    //  extra for unscaled
2947    double * unscaledCopy;
2948    unscaledCopy = lowerCopy + numberTotal; 
2949    memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
2950    memcpy(unscaledCopy+numberColumns_,
2951           rowLower_,numberRows_*sizeof(double));
2952    unscaledCopy = upperCopy + numberTotal; 
2953    memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
2954    memcpy(unscaledCopy+numberColumns_,
2955           rowUpper_,numberRows_*sizeof(double));
2956  }
2957  double maxTheta = 1.0e50;
2958  double largestChange=0.0;
2959  for (int iRow = 0; iRow < numberRows_; iRow++) {
2960    double lower = rowLower_[iRow];
2961    double upper = rowUpper_[iRow];
2962    if (lower<-1.0e30)
2963      lowerChange[numberColumns_+iRow]=0.0;
2964    double chgLower = lowerChange[numberColumns_+iRow];
2965    largestChange=CoinMax(largestChange,fabs(chgLower));
2966    if (upper>1.0e30)
2967      upperChange[numberColumns_+iRow]=0.0;
2968    double chgUpper = upperChange[numberColumns_+iRow];
2969    largestChange=CoinMax(largestChange,fabs(chgUpper));
2970    if (lower > -1.0e30 && upper < 1.0e30) {
2971      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2972        maxTheta = (upper - lower) / (chgLower - chgUpper);
2973      }
2974    }
2975    lower+=startingTheta*chgLower;
2976    upper+=startingTheta*chgUpper;
2977#ifndef CLP_USER_DRIVEN
2978    if (lower > upper) {
2979      maxTheta = -1.0;
2980      break;
2981    }
2982#endif
2983    rowLower_[iRow]=lower;
2984    rowUpper_[iRow]=upper;
2985  }
2986  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2987    double lower = columnLower_[iColumn];
2988    double upper = columnUpper_[iColumn];
2989    if (lower<-1.0e30)
2990      lowerChange[iColumn]=0.0;
2991    double chgLower = lowerChange[iColumn];
2992    largestChange=CoinMax(largestChange,fabs(chgLower));
2993    if (upper>1.0e30)
2994      upperChange[iColumn]=0.0;
2995    double chgUpper = upperChange[iColumn];
2996    largestChange=CoinMax(largestChange,fabs(chgUpper));
2997    if (lower > -1.0e30 && upper < 1.0e30) {
2998      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2999        maxTheta = (upper - lower) / (chgLower - chgUpper);
3000      }
3001    }
3002    lower+=startingTheta*chgLower;
3003    upper+=startingTheta*chgUpper;
3004#ifndef CLP_USER_DRIVEN
3005    if (lower > upper) {
3006      maxTheta = -1.0;
3007      break;
3008    }
3009#endif
3010    columnLower_[iColumn]=lower;
3011    columnUpper_[iColumn]=upper;
3012  }
3013  if (maxTheta == 1.0e50)
3014    maxTheta = COIN_DBL_MAX;
3015  int returnCode=0;
3016  /*
3017    If user is in charge then they are sophisticated enough
3018    to know what they are doing.  They must be doing something clever
3019    with event handlers!
3020   */
3021#ifndef CLP_USER_DRIVEN
3022  if (maxTheta < 0.0) {
3023    // bad ranges or initial
3024    returnCode = -1;
3025  }
3026  if (maxTheta < endingTheta) {
3027    char line[100];
3028    sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
3029            endingTheta,maxTheta);
3030    handler_->message(CLP_GENERAL,messages_)
3031      << line << CoinMessageEol;
3032    endingTheta = maxTheta;
3033  }
3034  if (endingTheta < startingTheta) {
3035    // bad initial
3036    returnCode = -2;
3037  }
3038  paramData.maxTheta=CoinMin(maxTheta,endingTheta);
3039#else
3040  // accept user's version
3041  paramData.maxTheta=endingTheta;
3042#endif
3043  /* given largest change element choose acceptable end
3044     be safe and make sure difference < 0.1*tolerance */
3045  double acceptableDifference=0.1*primalTolerance_/
3046    CoinMax(largestChange,1.0);
3047  paramData.acceptableMaxTheta=maxTheta-acceptableDifference;
3048  bool swapped=false;
3049  // Dantzig
3050#define ALL_DANTZIG
3051#ifdef ALL_DANTZIG
3052  ClpDualRowPivot * savePivot = dualRowPivot_;
3053  dualRowPivot_ = new ClpDualRowDantzig();
3054  dualRowPivot_->setModel(this);
3055#else
3056  ClpDualRowPivot * savePivot = NULL;
3057#endif
3058  if (!returnCode) {
3059    assert (objective_->type()==1);
3060    objective_->setType(2); // in case matrix empty
3061    returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
3062    objective_->setType(1);
3063    if (!returnCode) {
3064      double saveDualBound=dualBound_;
3065      dualBound_=CoinMax(dualBound_,1.0e15);
3066      swapped=true;
3067      double * temp;
3068      memcpy(saveLower,lower_,numberTotal*sizeof(double));
3069      temp=saveLower;
3070      saveLower=lower_;
3071      lower_=temp;
3072      //columnLowerWork_ = lower_;
3073      //rowLowerWork_ = lower_ + numberColumns_;
3074      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
3075      temp=saveUpper;
3076      saveUpper=upper_;
3077      upper_=temp;
3078      //columnUpperWork_ = upper_;
3079      //rowUpperWork_ = upper_ + numberColumns_;
3080      if (rowScale_) {
3081        // scale saved and change arrays
3082        double * lowerChange = lower_+numberTotal;
3083        double * upperChange = upper_+numberTotal;
3084        double * lowerSave = lowerChange+numberTotal;
3085        double * upperSave = upperChange+numberTotal;
3086        for (int i=0;i<numberColumns_;i++) {
3087          double multiplier = inverseColumnScale_[i];
3088          if (lowerSave[i]>-1.0e20)
3089            lowerSave[i] *= multiplier;
3090          if (upperSave[i]<1.0e20)
3091            upperSave[i] *= multiplier;
3092          lowerChange[i] *= multiplier;
3093          upperChange[i] *= multiplier;
3094        }
3095        lowerChange += numberColumns_;
3096        upperChange += numberColumns_;
3097        lowerSave += numberColumns_;
3098        upperSave += numberColumns_;
3099        for (int i=0;i<numberRows_;i++) {
3100          double multiplier = rowScale_[i];
3101          if (lowerSave[i]>-1.0e20)
3102            lowerSave[i] *= multiplier;
3103          if (upperSave[i]<1.0e20)
3104            upperSave[i] *= multiplier;
3105          lowerChange[i] *= multiplier;
3106          upperChange[i] *= multiplier;
3107        }
3108      }
3109      //double saveEndingTheta = endingTheta;
3110      double * saveDuals = NULL;
3111      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3112      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3113        // probably can get rid of this if we adjust every change in theta
3114        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3115        //   sumPrimalInfeasibilities_);
3116        int pass=100;
3117        while(sumPrimalInfeasibilities_) {
3118          pass--;
3119          if (!pass)
3120            break;
3121          problemStatus_=-1;
3122          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3123            double value=solution_[iSequence];
3124            // remember scaling
3125            if (value<lower_[iSequence]-1.0e-9) {
3126              lower_[iSequence]=value;
3127              lowerCopy[iSequence]=value;
3128            } else if (value>upper_[iSequence]+1.0e-9) {
3129              upper_[iSequence]=value;
3130              upperCopy[iSequence]=value;
3131            }
3132          }
3133          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3134        }
3135      }
3136      if (!problemStatus_) {
3137        if (nLowerChange||nUpperChange) {
3138#ifndef ALL_DANTZIG
3139          // Dantzig
3140          savePivot = dualRowPivot_;
3141          dualRowPivot_ = new ClpDualRowDantzig();
3142          dualRowPivot_->setModel(this);
3143#endif
3144          //for (int i=0;i<numberRows_+numberColumns_;i++)
3145          //setFakeBound(i, noFake);
3146          // Now do parametrics
3147          handler_->message(CLP_PARAMETRICS_STATS, messages_)
3148            << startingTheta << objectiveValue() << CoinMessageEol;
3149          bool canSkipFactorization=true;
3150          while (!returnCode) {
3151            paramData.startingTheta=startingTheta;
3152            paramData.endingTheta=endingTheta;
3153            returnCode = parametricsLoop(paramData,
3154                                         data,canSkipFactorization);
3155            startingTheta=paramData.startingTheta;
3156            endingTheta=paramData.endingTheta;
3157            canSkipFactorization=false;
3158            if (!returnCode) {
3159              //startingTheta = endingTheta;
3160              //endingTheta = saveEndingTheta;
3161              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3162                << startingTheta << objectiveValue() << CoinMessageEol;
3163              if (startingTheta >= endingTheta-primalTolerance_
3164                  ||problemStatus_==2)
3165                break;
3166            } else if (returnCode == -1) {
3167              // trouble - do external solve
3168              abort(); //needToDoSomething = true;
3169            } else if (problemStatus_==1) {
3170              // can't move any further
3171              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3172                << endingTheta << objectiveValue() << CoinMessageEol;
3173              problemStatus_=0;
3174            }
3175          }
3176        }
3177        dualBound_ = saveDualBound;
3178        //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3179      } else {
3180        // check if empty
3181        //if (!numberRows_||!matrix_->getNumElements()) {
3182        // success
3183#ifdef CLP_USER_DRIVEN
3184        //theta_ = endingTheta;
3185        //eventHandler_->event(ClpEventHandler::theta);
3186#endif
3187        //}
3188      }
3189    }
3190    if (problemStatus_==2) {
3191      delete [] ray_;
3192      ray_ = new double [numberColumns_];
3193    }
3194    if (swapped&&lower_) {
3195      double * temp=saveLower;
3196      saveLower=lower_;
3197      lower_=temp;
3198      temp=saveUpper;
3199      saveUpper=upper_;
3200      upper_=temp;
3201    }
3202    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
3203  }   
3204  if (!scalingFlag_) {
3205    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3206    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3207    memcpy(rowLower_,lowerCopy+numberColumns_,
3208           numberRows_*sizeof(double));
3209    memcpy(rowUpper_,upperCopy+numberColumns_,
3210           numberRows_*sizeof(double));
3211  } else {
3212    //  extra for unscaled
3213    double * unscaledCopy;
3214    unscaledCopy = lowerCopy + numberTotal; 
3215    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3216    memcpy(rowLower_,unscaledCopy+numberColumns_,
3217           numberRows_*sizeof(double));
3218    unscaledCopy = upperCopy + numberTotal; 
3219    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3220    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3221           numberRows_*sizeof(double));
3222  }
3223  delete [] saveLower;
3224  delete [] saveUpper;
3225#ifdef ALL_DANTZIG
3226  if (savePivot) {
3227#endif
3228    delete dualRowPivot_;
3229    dualRowPivot_ = savePivot;
3230#ifdef ALL_DANTZIG
3231  }
3232#endif
3233  // Restore any saved stuff
3234  restoreData(data);
3235  perturbation_ = savePerturbation;
3236  delete rowArray_[4];
3237  rowArray_[4]=NULL;
3238  delete rowArray_[5];
3239  rowArray_[5]=NULL;
3240  char line[100];
3241  sprintf(line,"Ending theta %g\n", endingTheta);
3242  handler_->message(CLP_GENERAL,messages_)
3243    << line << CoinMessageEol;
3244  return problemStatus_;
3245}
3246int
3247ClpSimplexOther::parametricsLoop(parametricsData & paramData,
3248                                 ClpDataSave & data,bool canSkipFactorization)
3249{
3250  double & startingTheta = paramData.startingTheta;
3251  double & endingTheta = paramData.endingTheta;
3252  int numberTotal = numberRows_+numberColumns_;
3253  // stuff is already at starting
3254  const int * lowerList = paramData.lowerList;
3255  const int * upperList = paramData.upperList;
3256  problemStatus_ = -1;
3257  //double saveEndingTheta=endingTheta;
3258
3259  // This says whether to restore things etc
3260  // startup will have factorized so can skip
3261  int factorType = 0;
3262  // Start check for cycles
3263  progress_.startCheck();
3264  // Say change made on first iteration
3265  changeMade_ = 1;
3266  /*
3267    Status of problem:
3268    0 - optimal
3269    1 - infeasible
3270    2 - unbounded
3271    -1 - iterating
3272    -2 - factorization wanted
3273    -3 - redo checking without factorization
3274    -4 - looks infeasible
3275  */
3276  while (problemStatus_ < 0) {
3277    int iRow, iColumn;
3278    // clear
3279    for (iRow = 0; iRow < 6; iRow++) {
3280      rowArray_[iRow]->clear();
3281    }
3282   
3283    for (iColumn = 0; iColumn < 2; iColumn++) {
3284      columnArray_[iColumn]->clear();
3285    }
3286   
3287    // give matrix (and model costs and bounds a chance to be
3288    // refreshed (normally null)
3289    matrix_->refresh(this);
3290    // may factorize, checks if problem finished
3291    if (!canSkipFactorization)
3292      statusOfProblemInParametrics(factorType, data);
3293    canSkipFactorization=false;
3294    if (numberPrimalInfeasibilities_) {
3295      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3296        // treat as success
3297        problemStatus_=0;
3298        endingTheta=startingTheta;
3299        break;
3300      }
3301      // probably can get rid of this if we adjust every change in theta
3302      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3303      //     sumPrimalInfeasibilities_);
3304      const double * lowerChange = lower_+numberTotal;
3305      const double * upperChange = upper_+numberTotal;
3306      const double * startLower = lowerChange+numberTotal;
3307      const double * startUpper = upperChange+numberTotal;
3308      //startingTheta -= 1.0e-7;
3309      int nLowerChange = lowerList[-1];
3310      for (int i = 0; i < nLowerChange; i++) {
3311        int iSequence = lowerList[i];
3312        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3313      }
3314      int nUpperChange = upperList[-1];
3315      for (int i = 0; i < nUpperChange; i++) {
3316        int iSequence = upperList[i];
3317        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3318      }
3319      // adjust rhs in case dual uses
3320      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3321      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3322      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3323      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3324      if (rowScale_) {
3325        for (int i=0;i<numberColumns_;i++) {
3326          double multiplier = columnScale_[i];
3327          if (columnLower_[i]>-1.0e20)
3328            columnLower_[i] *= multiplier;
3329          if (columnUpper_[i]<1.0e20)
3330            columnUpper_[i] *= multiplier;
3331        }
3332        for (int i=0;i<numberRows_;i++) {
3333          double multiplier = inverseRowScale_[i];
3334          if (rowLower_[i]>-1.0e20)
3335            rowLower_[i] *= multiplier;
3336          if (rowUpper_[i]<1.0e20)
3337            rowUpper_[i] *= multiplier;
3338        }
3339      }
3340      double * saveDuals = NULL;
3341      problemStatus_=-1;
3342      ClpObjective * saveObjective = objective_;
3343      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3344      if (saveObjective!=objective_) {
3345        delete objective_;
3346        objective_=saveObjective;
3347      }
3348      int pass=100;
3349      double moved=0.0;
3350      while(sumPrimalInfeasibilities_) {
3351        //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_,
3352        //     sumPrimalInfeasibilities_);
3353        pass--;
3354        if (!pass)
3355          break;
3356        problemStatus_=-1;
3357        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3358          double value=solution_[iSequence];
3359          if (value<lower_[iSequence]-1.0e-9) {
3360            moved += lower_[iSequence]-value;
3361            lower_[iSequence]=value;
3362          } else if (value>upper_[iSequence]+1.0e-9) {
3363            moved += upper_[iSequence]-value;
3364            upper_[iSequence]=value;
3365          }
3366        }
3367        if (!moved) {
3368          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3369            double value=solution_[iSequence];
3370            if (value<lower_[iSequence]-1.0e-9) {
3371              moved += lower_[iSequence]-value;
3372              lower_[iSequence]=value;
3373            } else if (value>upper_[iSequence]+1.0e-9) {
3374              moved += upper_[iSequence]-value;
3375              upper_[iSequence]=value;
3376            }
3377          }
3378        }
3379        assert (moved);
3380        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3381      }
3382      // adjust
3383      //printf("Should adjust - moved %g\n",moved);
3384    }
3385    // Say good factorization
3386    factorType = 1;
3387    if (data.sparseThreshold_) {
3388      // use default at present
3389      factorization_->sparseThreshold(0);
3390      factorization_->goSparse();
3391    }
3392   
3393    // exit if victory declared
3394    if ((problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7) ||
3395        (problemStatus_ == 0 && startingTheta>=CoinMin(endingTheta-1.0e-7,1.0e49) )) {
3396      startingTheta=endingTheta;
3397      break;
3398    }
3399   
3400    // test for maximum iterations
3401    if (hitMaximumIterations()) {
3402      problemStatus_ = 3;
3403      break;
3404    }
3405#ifdef CLP_USER_DRIVEN
3406    // Check event
3407    {
3408      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3409      if (status >= 0) {
3410        problemStatus_ = 5;
3411        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3412        break;
3413      }
3414    }
3415#endif
3416    // Do iterations
3417    problemStatus_=-1;
3418    whileIterating(paramData, 0.0,
3419                   NULL);
3420    //startingTheta = endingTheta;
3421    //endingTheta = saveEndingTheta;
3422  }
3423  if (!problemStatus_/*||problemStatus_==2*/) {
3424    theta_ = endingTheta;
3425#ifdef CLP_USER_DRIVEN
3426    {
3427      double saveTheta=theta_;
3428      theta_ = endingTheta;
3429      int status=eventHandler_->event(ClpEventHandler::theta);
3430      if (status>=0&&status<10) {
3431        endingTheta=theta_;
3432        theta_=saveTheta;
3433        problemStatus_=-1;
3434      } else {
3435        if (status>=10) {
3436          problemStatus_=status-10;
3437          startingTheta=endingTheta;
3438        }
3439        theta_=saveTheta;
3440      }
3441    }
3442#endif
3443    return 0;
3444  } else if (problemStatus_ == 10) {
3445    return -1;
3446  } else {
3447    return problemStatus_;
3448  }
3449}
3450/* Checks if finished.  Updates status */
3451void
3452ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3453{
3454     if (type == 2) {
3455          // trouble - go to recovery
3456          problemStatus_ = 10;
3457          return;
3458     }
3459     if (problemStatus_ > -3 || factorization_->pivots()) {
3460          // factorize
3461          // later on we will need to recover from singularities
3462          // also we could skip if first time
3463          if (type) {
3464               // is factorization okay?
3465               if (internalFactorize(1)) {
3466                    // trouble - go to recovery
3467                    problemStatus_ = 10;
3468                    return;
3469               }
3470          }
3471          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3472               problemStatus_ = -3;
3473     }
3474     // at this stage status is -3 or -4 if looks infeasible
3475     // get primal and dual solutions
3476     gutsOfSolution(NULL, NULL);
3477     double realDualInfeasibilities = sumDualInfeasibilities_;
3478     // If bad accuracy treat as singular
3479     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3480          // trouble - go to recovery
3481          problemStatus_ = 10;
3482          return;
3483     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3484          // Can reduce tolerance
3485          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3486          factorization_->pivotTolerance(newTolerance);
3487     }
3488     // Check if looping
3489     int loop;
3490     if (type != 2)
3491          loop = progress_.looping();
3492     else
3493          loop = -1;
3494     if (loop >= 0) {
3495          problemStatus_ = loop; //exit if in loop
3496          if (!problemStatus_) {
3497               // declaring victory
3498               numberPrimalInfeasibilities_ = 0;
3499               sumPrimalInfeasibilities_ = 0.0;
3500          } else {
3501               problemStatus_ = 10; // instead - try other algorithm
3502          }
3503          return;
3504     } else if (loop < -1) {
3505          // something may have changed
3506          gutsOfSolution(NULL, NULL);
3507     }
3508     progressFlag_ = 0; //reset progress flag
3509     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3510          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3511                    << numberIterations_ << objectiveValue();
3512          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3513                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3514          handler_->printing(sumDualInfeasibilities_ > 0.0)
3515                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3516          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3517                             < numberDualInfeasibilities_)
3518                    << numberDualInfeasibilitiesWithoutFree_;
3519          handler_->message() << CoinMessageEol;
3520     }
3521#ifdef CLP_USER_DRIVEN
3522     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3523       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3524       if (status>=0) {
3525         // fix up
3526         for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3527           double value=solution_[iSequence];
3528           if (value<=lower_[iSequence]-primalTolerance_) {
3529             lower_[iSequence]=value;
3530           } else if (value>=upper_[iSequence]+primalTolerance_) {
3531             upper_[iSequence]=value;
3532           }
3533         }
3534         numberPrimalInfeasibilities_ = 0;
3535         sumPrimalInfeasibilities_ = 0.0;
3536       }
3537     }
3538#endif
3539     /* If we are primal feasible and any dual infeasibilities are on
3540        free variables then it is better to go to primal */
3541     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3542               numberDualInfeasibilities_) {
3543          problemStatus_ = 10;
3544          return;
3545     }
3546
3547     // check optimal
3548     // give code benefit of doubt
3549     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3550               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3551          // say optimal (with these bounds etc)
3552          numberDualInfeasibilities_ = 0;
3553          sumDualInfeasibilities_ = 0.0;
3554          numberPrimalInfeasibilities_ = 0;
3555          sumPrimalInfeasibilities_ = 0.0;
3556     }
3557     if (dualFeasible() || problemStatus_ == -4) {
3558          progress_.modifyObjective(objectiveValue_
3559                                    - sumDualInfeasibilities_ * dualBound_);
3560     }
3561     if (numberPrimalInfeasibilities_) {
3562          if (problemStatus_ == -4 || problemStatus_ == -5) {
3563               problemStatus_ = 1; // infeasible
3564          }
3565     } else if (numberDualInfeasibilities_) {
3566          // clean up
3567          problemStatus_ = 10;
3568     } else {
3569          problemStatus_ = 0;
3570     }
3571     lastGoodIteration_ = numberIterations_;
3572     if (problemStatus_ < 0) {
3573          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3574          if (sumDualInfeasibilities_)
3575               numberDualInfeasibilities_ = 1;
3576     }
3577     // Allow matrices to be sorted etc
3578     int fake = -999; // signal sort
3579     matrix_->correctSequence(this, fake, fake);
3580}
3581//static double lastThetaX=0.0;
3582/* This has the flow between re-factorizations
3583   Reasons to come out:
3584   -1 iterations etc
3585   -2 inaccuracy
3586   -3 slight inaccuracy (and done iterations)
3587   +0 looks optimal (might be unbounded - but we will investigate)
3588   +1 looks infeasible
3589   +3 max iterations
3590   +4 accuracy problems
3591*/
3592int
3593ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
3594                                const double * /*changeObjective*/)
3595{
3596  double & startingTheta = paramData.startingTheta;
3597  double & endingTheta = paramData.endingTheta;
3598  const double * lowerChange = paramData.lowerChange;
3599  const double * upperChange = paramData.upperChange;
3600  int numberTotal = numberColumns_ + numberRows_;
3601  const int * lowerList = paramData.lowerList;
3602  const int * upperList = paramData.upperList;
3603  //#define CLP_PARAMETRIC_DENSE_ARRAYS 2
3604#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3605  double * lowerGap = paramData.lowerGap;
3606  double * upperGap = paramData.upperGap;
3607  double * lowerCoefficient = paramData.lowerCoefficient;
3608  double * upperCoefficient = paramData.upperCoefficient;
3609#endif
3610  // do basic pointers
3611  int * backwardBasic = paramData.backwardBasic;
3612  for (int i=0;i<numberTotal;i++)
3613    backwardBasic[i]=-1;
3614  for (int i=0;i<numberRows_;i++) {
3615    int iPivot=pivotVariable_[i];
3616    backwardBasic[iPivot]=i;
3617  }
3618     {
3619          int i;
3620          for (i = 0; i < 4; i++) {
3621               rowArray_[i]->clear();
3622          }
3623          for (i = 0; i < 2; i++) {
3624               columnArray_[i]->clear();
3625          }
3626     }
3627     // if can't trust much and long way from optimal then relax
3628     if (largestPrimalError_ > 10.0)
3629          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3630     else
3631          factorization_->relaxAccuracyCheck(1.0);
3632     // status stays at -1 while iterating, >=0 finished, -2 to invert
3633     // status -3 to go to top without an invert
3634     int returnCode = -1;
3635     double lastTheta = startingTheta;
3636     double useTheta = startingTheta;
3637     while (problemStatus_ == -1) {
3638          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3639          // Get theta for bounds - we know can't crossover
3640          int pivotType = nextTheta(1, increaseTheta, paramData,
3641                                     NULL);
3642          useTheta += theta_;
3643          double change = useTheta - lastTheta;
3644          if (change>1.0e-14) {
3645            int n;
3646            n=lowerList[-1];
3647            for (int i=0;i<n;i++) {
3648              int iSequence = lowerList[i];
3649              double thisChange = change * lowerChange[iSequence];
3650              double newValue = lower_[iSequence] + thisChange;
3651              lower_[iSequence] = newValue;
3652#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3653              if (getStatus(iSequence)==basic) {
3654                int iRow=backwardBasic[iSequence];
3655                lowerGap[iRow] -= thisChange;
3656              } else if(getStatus(iSequence)==atLowerBound) {
3657                solution_[iSequence] = newValue;
3658              }
3659#else
3660              if(getStatus(iSequence)==atLowerBound) {
3661                solution_[iSequence] = newValue;
3662              }
3663#endif
3664#if 0
3665              // may have to adjust other bound
3666              double otherValue = upper_[iSequence];
3667              if (otherValue-newValue<dualBound_) {
3668                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3669                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3670                //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5);
3671              }
3672#endif
3673            }
3674            n=upperList[-1];
3675            for (int i=0;i<n;i++) {
3676              int iSequence = upperList[i];
3677              double thisChange = change * upperChange[iSequence];
3678              double newValue = upper_[iSequence] + thisChange;
3679              upper_[iSequence] = newValue;
3680              if(getStatus(iSequence)==atUpperBound||
3681                 getStatus(iSequence)==isFixed) {
3682                solution_[iSequence] = newValue;
3683#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3684              } else if (getStatus(iSequence)==basic) {
3685                int iRow=backwardBasic[iSequence];
3686                upperGap[iRow] += thisChange;
3687#endif
3688              }
3689              // may have to adjust other bound
3690              double otherValue = lower_[iSequence];
3691              if (newValue-otherValue<dualBound_) {
3692                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3693                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3694                //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5);
3695              }
3696            }
3697          }
3698          sequenceIn_=-1;
3699          if (pivotType) {
3700            if (useTheta>lastTheta+1.0e-9) {
3701              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3702                << useTheta << objectiveValue() << CoinMessageEol;
3703              lastTheta = useTheta;
3704            }
3705            problemStatus_ = -2;
3706            if (!factorization_->pivots()&&pivotRow_<0)
3707              problemStatus_=2;
3708#ifdef CLP_USER_DRIVEN
3709            {
3710              double saveTheta=theta_;
3711              theta_ = endingTheta;
3712              if (problemStatus_==2&&theta_>paramData.acceptableMaxTheta)
3713                theta_=COIN_DBL_MAX; // we have finished
3714              int status=eventHandler_->event(ClpEventHandler::theta);
3715              if (status>=0&&status<10) {
3716                endingTheta=theta_;
3717                problemStatus_=-1;
3718                continue;
3719              } else {
3720                if (status>=10)
3721                  problemStatus_=status-10;
3722                if (status<0)
3723                  startingTheta = useTheta;
3724              }
3725              theta_=saveTheta;
3726            }
3727#else
3728            startingTheta = useTheta;
3729#endif
3730            return 4;
3731          }
3732          // choose row to go out
3733          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3734          if (pivotRow_ >= 0) {
3735               // we found a pivot row
3736               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3737                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3738                              << pivotRow_
3739                              << CoinMessageEol;
3740               }
3741               // check accuracy of weights
3742               dualRowPivot_->checkAccuracy();
3743               // do ratio test for normal iteration
3744               double bestPossiblePivot = bestPivot();
3745               if (sequenceIn_ >= 0) {
3746                    // normal iteration
3747                    // update the incoming column
3748                    double btranAlpha = -alpha_ * directionOut_; // for check
3749#ifndef COIN_FAC_NEW
3750                    unpackPacked(rowArray_[1]);
3751#else
3752                    unpack(rowArray_[1]);
3753#endif
3754                    // and update dual weights (can do in parallel - with extra array)
3755                    rowArray_[2]->clear();
3756                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3757                                                          rowArray_[2],
3758                                                          rowArray_[3],
3759                                                          rowArray_[1]);
3760                    // see if update stable
3761#ifdef CLP_DEBUG
3762                    if ((handler_->logLevel() & 32))
3763                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3764#endif
3765                    double checkValue = 1.0e-7;
3766                    // if can't trust much and long way from optimal then relax
3767                    if (largestPrimalError_ > 10.0)
3768                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3769                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3770                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3771                         handler_->message(CLP_DUAL_CHECK, messages_)
3772                                   << btranAlpha
3773                                   << alpha_
3774                                   << CoinMessageEol;
3775                         // clear arrays
3776                         rowArray_[4]->clear();
3777                         if (factorization_->pivots()) {
3778                              dualRowPivot_->unrollWeights();
3779                              problemStatus_ = -2; // factorize now
3780                              rowArray_[0]->clear();
3781                              rowArray_[1]->clear();
3782                              columnArray_[0]->clear();
3783                              returnCode = -2;
3784                              break;
3785                         } else {
3786                              // take on more relaxed criterion
3787                              double test;
3788                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3789                                   test = 1.0e-1 * fabs(alpha_);
3790                              else
3791                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
3792                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3793                                        fabs(btranAlpha - alpha_) > test) {
3794                                   dualRowPivot_->unrollWeights();
3795                                   // need to reject something
3796                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
3797                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
3798                                             << x << sequenceWithin(sequenceOut_)
3799                                             << CoinMessageEol;
3800                                   setFlagged(sequenceOut_);
3801                                   progress_.clearBadTimes();
3802                                   lastBadIteration_ = numberIterations_; // say be more cautious
3803                                   rowArray_[0]->clear();
3804                                   rowArray_[1]->clear();
3805                                   columnArray_[0]->clear();
3806                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3807                                        //printf("I think should declare infeasible\n");
3808                                        problemStatus_ = 1;
3809                                        returnCode = 1;
3810                                        break;
3811                                   }
3812                                   continue;
3813                              }
3814                         }
3815                    }
3816                    // update duals BEFORE replaceColumn so can do updateColumn
3817                    double objectiveChange = 0.0;
3818                    // do duals first as variables may flip bounds
3819                    // rowArray_[0] and columnArray_[0] may have flips
3820                    // so use rowArray_[3] for work array from here on
3821                    int nswapped = 0;
3822                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3823                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3824#if CLP_CAN_HAVE_ZERO_OBJ
3825                    if ((specialOptions_&2097152)==0) {
3826#endif
3827                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3828                                                                                               rowArray_[2], theta_,
3829                                                                                               objectiveChange, false);
3830                      assert (!nswapped);
3831#if CLP_CAN_HAVE_ZERO_OBJ
3832                    } else {
3833                      rowArray_[0]->clear();
3834                      rowArray_[2]->clear();
3835                      columnArray_[0]->clear();
3836                    }
3837#endif
3838                    // which will change basic solution
3839                    if (nswapped) {
3840                      abort(); //needs testing
3841                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3842                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
3843                                                             1.0, objectiveChange);
3844                         // recompute dualOut_
3845                         valueOut_ = solution_[sequenceOut_];
3846                         if (directionOut_ < 0) {
3847                              dualOut_ = valueOut_ - upperOut_;
3848                         } else {
3849                              dualOut_ = lowerOut_ - valueOut_;
3850                         }
3851                    }
3852                    // amount primal will move
3853                    double movement = -dualOut_ * directionOut_ / alpha_;
3854                    // so objective should increase by fabs(dj)*movement
3855                    // but we already have objective change - so check will be good
3856                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3857#ifdef CLP_DEBUG
3858                         if (handler_->logLevel() & 32)
3859                              printf("movement %g, swap change %g, rest %g  * %g\n",
3860                                     objectiveChange + fabs(movement * dualIn_),
3861                                     objectiveChange, movement, dualIn_);
3862#endif
3863                         assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3864                         if(factorization_->pivots()) {
3865                              // going backwards - factorize
3866                              dualRowPivot_->unrollWeights();
3867                              problemStatus_ = -2; // factorize now
3868                              returnCode = -2;
3869                              break;
3870                         }
3871                    }
3872                    CoinAssert(fabs(dualOut_) < 1.0e50);
3873                    // if stable replace in basis
3874                    int updateStatus = factorization_->replaceColumn(this,
3875                                       rowArray_[2],
3876                                       rowArray_[1],
3877                                       pivotRow_,
3878                                       alpha_);
3879                    // if no pivots, bad update but reasonable alpha - take and invert
3880                    if (updateStatus == 2 &&
3881                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3882                         updateStatus = 4;
3883                    if (updateStatus == 1 || updateStatus == 4) {
3884                         // slight error
3885                         if (factorization_->pivots() > 5 || updateStatus == 4) {
3886                              problemStatus_ = -2; // factorize now
3887                              returnCode = -3;
3888                         }
3889                    } else if (updateStatus == 2) {
3890                         // major error
3891                         dualRowPivot_->unrollWeights();
3892                         // later we may need to unwind more e.g. fake bounds
3893                         if (factorization_->pivots()) {
3894                              problemStatus_ = -2; // factorize now
3895                              returnCode = -2;
3896                              break;
3897                         } else {
3898                              // need to reject something
3899                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3900                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3901                                        << x << sequenceWithin(sequenceOut_)
3902                                        << CoinMessageEol;
3903                              setFlagged(sequenceOut_);
3904                              progress_.clearBadTimes();
3905                              lastBadIteration_ = numberIterations_; // say be more cautious
3906                              rowArray_[0]->clear();
3907                              rowArray_[1]->clear();
3908                              columnArray_[0]->clear();
3909                              // make sure dual feasible
3910                              // look at all rows and columns
3911                              double objectiveChange = 0.0;
3912                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3913                                        0.0, objectiveChange, true);
3914                              continue;
3915                         }
3916                    } else if (updateStatus == 3) {
3917                         // out of memory
3918                         // increase space if not many iterations
3919                         if (factorization_->pivots() <
3920                                   0.5 * factorization_->maximumPivots() &&
3921                                   factorization_->pivots() < 200)
3922                              factorization_->areaFactor(
3923                                   factorization_->areaFactor() * 1.1);
3924                         problemStatus_ = -2; // factorize now
3925                    } else if (updateStatus == 5) {
3926                         problemStatus_ = -2; // factorize now
3927                    }
3928                    int * lowerActive = paramData.lowerActive;
3929                    int * upperActive = paramData.upperActive;
3930                    // update change vector
3931                    {
3932                      double * work = rowArray_[1]->denseVector();
3933                      int number = rowArray_[1]->getNumElements();
3934                      int * which = rowArray_[1]->getIndices();
3935                      assert (!rowArray_[4]->packedMode());
3936#ifndef COIN_FAC_NEW
3937                      assert (rowArray_[1]->packedMode());
3938#else
3939                      assert (!rowArray_[1]->packedMode());
3940#endif
3941                      double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3942                      double multiplier = -pivotValue/alpha_;
3943                      double * array=rowArray_[4]->denseVector();
3944#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3945                      int lowerN=lowerActive[-1];
3946                      int upperN=upperActive[-1];
3947#endif
3948                      if (multiplier) {
3949                        for (int i = 0; i < number; i++) {
3950                          int iRow = which[i];
3951#ifndef COIN_FAC_NEW
3952                          double alpha=multiplier*work[i];
3953#else
3954                          double alpha=multiplier*work[iRow];
3955#endif
3956#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3957                          double alpha3 = alpha+array[iRow];
3958                          int iSequence = pivotVariable_[iRow];
3959                          double oldLower = lowerCoefficient[iRow];
3960                          double oldUpper = upperCoefficient[iRow];
3961                          if (lower_[iSequence]>-1.0e30) {
3962                            //lowerGap[iRow]=value-lower_[iSequence];
3963                            double alpha2 = alpha3 + lowerChange[iSequence];
3964                            if (alpha2>1.0e-8)  {
3965                              lowerCoefficient[iRow]=alpha2;
3966                              if (!oldLower)
3967                                lowerActive[lowerN++]=iRow;
3968                            } else {
3969                              if (oldLower)
3970                                lowerCoefficient[iRow]=COIN_DBL_MIN;
3971                            }
3972                          } else {
3973                            if (oldLower)
3974                              lowerCoefficient[iRow]=COIN_DBL_MIN;
3975                          }
3976                          if (upper_[iSequence]<1.0e30) {
3977                            //upperGap[iRow]=-(value-upper_[iSequence]);
3978                            double alpha2 = -(alpha3+upperChange[iSequence]);
3979                            if (alpha2>1.0e-8) {
3980                              upperCoefficient[iRow]=alpha2;
3981                              if (!oldUpper)
3982                                upperActive[upperN++]=iRow;
3983                            } else {
3984                              if (oldUpper)
3985                                upperCoefficient[iRow]=COIN_DBL_MIN;
3986                            }
3987                          } else {
3988                            if (oldUpper)
3989                              upperCoefficient[iRow]=COIN_DBL_MIN;
3990                          }
3991#endif
3992                          rowArray_[4]->quickAdd(iRow,alpha);
3993                        }
3994                      }
3995                      pivotValue = array[pivotRow_];
3996                      // we want pivot to be -multiplier
3997                      rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3998#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3999                      assert (lowerN>=0&&lowerN<=numberRows_);
4000                      lowerActive[-1]=lowerN;
4001                      upperActive[-1]=upperN;
4002#endif
4003                    }
4004                    // update primal solution
4005#if CLP_CAN_HAVE_ZERO_OBJ
4006                    if ((specialOptions_&2097152)!=0) 
4007                      theta_=0.0;
4008#endif
4009                    if (theta_ < 0.0) {
4010#ifdef CLP_DEBUG
4011                         if (handler_->logLevel() & 32)
4012                              printf("negative theta %g\n", theta_);
4013#endif
4014                         theta_ = 0.0;
4015                    }
4016                    // do actual flips
4017                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
4018                    //rowArray_[1]->expand();
4019#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4020                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
4021                                                        movement,
4022                                                        objectiveChange);
4023#else
4024                    // do by hand
4025                    {
4026                      double * work = rowArray_[1]->denseVector();
4027                      int number = rowArray_[1]->getNumElements();
4028                      int * which = rowArray_[1]->getIndices();
4029                      int i;
4030                      if (rowArray_[1]->packedMode()) {
4031                        for (i = 0; i < number; i++) {
4032                          int iRow = which[i];
4033                          int iSequence = pivotVariable_[iRow];
4034                          double value = solution_[iSequence];
4035                          double change = movement * work[i];
4036                          value -= change;
4037                          if (lower_[iSequence]>-1.0e30)
4038                            lowerGap[iRow]=value-lower_[iSequence];
4039                          if (upper_[iSequence]<1.0e30)
4040                            upperGap[iRow]=-(value-upper_[iSequence]);
4041                          solution_[iSequence] = value;
4042                          objectiveChange -= change * cost_[iSequence];
4043                          work[i] = 0.0;
4044                        }
4045                      } else {
4046                        for (i = 0; i < number; i++) {
4047                          int iRow = which[i];
4048                          int iSequence = pivotVariable_[iRow];
4049                          double value = solution_[iSequence];
4050                          double change = movement * work[iRow];
4051                          value -= change;
4052                          solution_[iSequence] = value;
4053                          objectiveChange -= change * cost_[iSequence];
4054                          work[iRow] = 0.0;
4055                        }
4056                      }
4057                      rowArray_[1]->setNumElements(0);
4058                    }
4059#endif
4060                    // modify dualout
4061                    dualOut_ /= alpha_;
4062                    dualOut_ *= -directionOut_;
4063                    //setStatus(sequenceIn_,basic);
4064                    dj_[sequenceIn_] = 0.0;
4065                    //double oldValue = valueIn_;
4066                    if (directionIn_ == -1) {
4067                         // as if from upper bound
4068                         valueIn_ = upperIn_ + dualOut_;
4069                    } else {
4070                         // as if from lower bound
4071                         valueIn_ = lowerIn_ + dualOut_;
4072                    }
4073                    objectiveChange = 0.0;
4074#if CLP_CAN_HAVE_ZERO_OBJ
4075                    if ((specialOptions_&2097152)==0) {
4076#endif
4077                      for (int i=0;i<numberTotal;i++)
4078                        objectiveChange += solution_[i]*cost_[i];
4079                      objectiveChange -= objectiveValue_;
4080#if CLP_CAN_HAVE_ZERO_OBJ
4081                    }
4082#endif
4083                    // outgoing
4084                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
4085                    lowerOut_=lower_[sequenceOut_];
4086                    upperOut_=upper_[sequenceOut_];
4087                    // set dj to zero unless values pass
4088                    if (directionOut_ > 0) {
4089                         valueOut_ = lowerOut_;
4090                         dj_[sequenceOut_] = theta_;
4091#if CLP_CAN_HAVE_ZERO_OBJ>1
4092#ifdef COIN_REUSE_RANDOM
4093                         if ((specialOptions_&2097152)!=0) {
4094                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
4095                         }
4096#endif
4097#endif
4098                    } else {
4099                         valueOut_ = upperOut_;
4100                         dj_[sequenceOut_] = -theta_;
4101#if CLP_CAN_HAVE_ZERO_OBJ>1
4102#ifdef COIN_REUSE_RANDOM
4103                         if ((specialOptions_&2097152)!=0) {
4104                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
4105                         }
4106#endif
4107#endif
4108                    }
4109                    solution_[sequenceOut_] = valueOut_;
4110                    int whatNext = housekeeping(objectiveChange);
4111                    reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_);
4112                    assert (backwardBasic[sequenceOut_]==pivotRow_);
4113                    backwardBasic[sequenceOut_]=-1;
4114                    backwardBasic[sequenceIn_]=pivotRow_;
4115#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4116                    double value = solution_[sequenceIn_];
4117                    double alpha = rowArray_[4]->denseVector()[pivotRow_];
4118                    double oldLower = lowerCoefficient[pivotRow_];
4119                    double oldUpper = upperCoefficient[pivotRow_];
4120                    if (lower_[sequenceIn_]>-1.0e30) {
4121                      lowerGap[pivotRow_]=value-lower_[sequenceIn_];
4122                      double alpha2 = alpha + lowerChange[sequenceIn_];
4123                      if (alpha2>1.0e-8)  {
4124                        lowerCoefficient[pivotRow_]=alpha2;
4125                        if (!oldLower) {
4126                          int lowerN=lowerActive[-1];
4127                          assert (lowerN>=0&&lowerN<numberRows_);
4128                          lowerActive[lowerN]=pivotRow_;
4129                          lowerActive[-1]=lowerN+1;
4130                        }
4131                      } else {
4132                        if (oldLower)
4133                          lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4134                      }
4135                    } else {
4136                      if (oldLower)
4137                        lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4138                    }
4139                    if (upper_[sequenceIn_]<1.0e30) {
4140                      upperGap[pivotRow_]=-(value-upper_[sequenceIn_]);
4141                      double alpha2 = -(alpha+upperChange[sequenceIn_]);
4142                      if (alpha2>1.0e-8) {
4143                        upperCoefficient[pivotRow_]=alpha2;
4144                        if (!oldUpper) {
4145                          int upperN=upperActive[-1];
4146                          assert (upperN>=0&&upperN<numberRows_);
4147                          upperActive[upperN]=pivotRow_;
4148                          upperActive[-1]=upperN+1;
4149                        }
4150                      } else {
4151                        if (oldUpper)
4152                          upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4153                      }
4154                    } else {
4155                      if (oldUpper)
4156                        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4157                    }
4158#endif
4159                    {
4160                      char in[200],out[200];
4161                      int iSequence=sequenceIn_;
4162                      if (iSequence<numberColumns_) {
4163                        if (lengthNames_) 
4164                          strcpy(in,columnNames_[iSequence].c_str());
4165                         else 
4166                          sprintf(in,"C%7.7d",iSequence);
4167                      } else {
4168                        iSequence -= numberColumns_;
4169                        if (lengthNames_) 
4170                          strcpy(in,rowNames_[iSequence].c_str());
4171                         else 
4172                          sprintf(in,"R%7.7d",iSequence);
4173                      }
4174                      iSequence=sequenceOut_;
4175                      if (iSequence<numberColumns_) {
4176                        if (lengthNames_) 
4177                          strcpy(out,columnNames_[iSequence].c_str());
4178                         else 
4179                          sprintf(out,"C%7.7d",iSequence);
4180                      } else {
4181                        iSequence -= numberColumns_;
4182                        if (lengthNames_) 
4183                          strcpy(out,rowNames_[iSequence].c_str());
4184                         else 
4185                          sprintf(out,"R%7.7d",iSequence);
4186                      }
4187                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
4188                        << useTheta << objectiveValue() 
4189                        << in << out << CoinMessageEol;
4190                    }
4191                    if (useTheta>lastTheta+1.0e-9) {
4192                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
4193                        << useTheta << objectiveValue() << CoinMessageEol;
4194                      lastTheta = useTheta;
4195                    }
4196                    // and set bounds correctly
4197                    originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
4198                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
4199                    if (whatNext == 1) {
4200                         problemStatus_ = -2; // refactorize
4201                    } else if (whatNext == 2) {
4202                         // maximum iterations or equivalent
4203                         problemStatus_ = 3;
4204                         returnCode = 3;
4205                         break;
4206                    }
4207#ifdef CLP_USER_DRIVEN
4208                    // Check event
4209                    {
4210                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
4211                         if (status >= 0) {
4212                              problemStatus_ = 5;
4213                              secondaryStatus_ = ClpEventHandler::endOfIteration;
4214                              returnCode = 4;
4215                              break;
4216                         }
4217                    }
4218#endif
4219               } else {
4220                    // no incoming column is valid
4221#ifdef CLP_USER_DRIVEN
4222                 rowArray_[0]->clear();
4223                 columnArray_[0]->clear();
4224                 theta_ = useTheta;
4225                 lastTheta = useTheta;
4226                 int action = eventHandler_->event(ClpEventHandler::noTheta);
4227                 if (action>=0) {
4228                   endingTheta=theta_;
4229                   theta_ = 0.0;
4230                   //adjust [4] from handler - but
4231                   //rowArray_[4]->clear(); // temp
4232                   if (action>=0&&action<10)
4233                     problemStatus_=-1; // carry on
4234                   else if (action==15)
4235                     problemStatus_ =5; // say stopped
4236                   returnCode = 1;
4237                   if (action==0||action>=10) 
4238                     break;
4239                   else
4240                     continue;
4241                 } else {
4242                 theta_ = 0.0;
4243                 }
4244#endif
4245                    pivotRow_ = -1;
4246#ifdef CLP_DEBUG
4247                    if (handler_->logLevel() & 32)
4248                         printf("** no column pivot\n");
4249#endif
4250                    if (factorization_->pivots() < 10) { 
4251                         // If we have just factorized and infeasibility reasonable say infeas
4252                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
4253                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
4254                                        || (specialOptions_ & 64) == 0) {
4255                                   // say infeasible
4256                                   problemStatus_ = 1;
4257                                   // unless primal feasible!!!!
4258                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
4259                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
4260                                   if (numberDualInfeasibilities_)
4261                                        problemStatus_ = 10;
4262                                   rowArray_[0]->clear();
4263                                   columnArray_[0]->clear();
4264                              }
4265                         }
4266                         // If special option set - put off as long as possible
4267                         if ((specialOptions_ & 64) == 0) {
4268                              problemStatus_ = -4; //say looks infeasible
4269                         } else {
4270                              // flag
4271                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
4272                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
4273                                        << x << sequenceWithin(sequenceOut_)
4274                                        << CoinMessageEol;
4275                              setFlagged(sequenceOut_);
4276                              if (!factorization_->pivots()) {
4277                                   rowArray_[0]->clear();
4278                                   columnArray_[0]->clear();
4279                                   continue;
4280                              }
4281                         }
4282                    }
4283                    rowArray_[0]->clear();
4284                    columnArray_[0]->clear();
4285                    returnCode = 1;
4286                    break;
4287               }
4288          } else {
4289               // no pivot row
4290#ifdef CLP_USER_DRIVEN
4291            {
4292              double saveTheta=theta_;
4293              theta_ = endingTheta;
4294              int status=eventHandler_->event(ClpEventHandler::theta);
4295              if (status>=0&&status<10) {
4296                endingTheta=theta_;
4297                theta_=saveTheta;
4298                continue;
4299              } else {
4300                theta_=saveTheta;
4301              }
4302            }
4303#endif
4304#ifdef CLP_DEBUG
4305               if (handler_->logLevel() & 32)
4306                    printf("** no row pivot\n");
4307#endif
4308               int numberPivots = factorization_->pivots();
4309               bool specialCase;
4310               int useNumberFake;
4311               returnCode = 0;
4312               if (numberPivots < 20 &&
4313                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
4314                         && dualBound_ > 1.0e8) {
4315                    specialCase = true;
4316                    // as dual bound high - should be okay
4317                    useNumberFake = 0;
4318               } else {
4319                    specialCase = false;
4320                    useNumberFake = numberFake_;
4321               }
4322               if (!numberPivots || specialCase) {
4323                    // may have crept through - so may be optimal
4324                    // check any flagged variables
4325                    int iRow;
4326                    for (iRow = 0; iRow < numberRows_; iRow++) {
4327                         int iPivot = pivotVariable_[iRow];
4328                         if (flagged(iPivot))
4329                              break;
4330                    }
4331                    if (iRow < numberRows_ && numberPivots) {
4332                         // try factorization
4333                         returnCode = -2;
4334                    }
4335
4336                    if (useNumberFake || numberDualInfeasibilities_) {
4337                         // may be dual infeasible
4338                         problemStatus_ = -5;
4339                    } else {
4340                         if (iRow < numberRows_) {
4341                              problemStatus_ = -5;
4342                         } else {
4343                              if (numberPivots) {
4344                                   // objective may be wrong
4345                                   objectiveValue_ = innerProduct(cost_,
4346                                                                  numberColumns_ + numberRows_,
4347                                                                  solution_);
4348                                   objectiveValue_ += objective_->nonlinearOffset();
4349                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
4350                                   if ((specialOptions_ & 16384) == 0) {
4351                                        // and dual_ may be wrong (i.e. for fixed or basic)
4352                                        CoinIndexedVector * arrayVector = rowArray_[1];
4353                                        arrayVector->clear();
4354                                        int iRow;
4355                                        double * array = arrayVector->denseVector();
4356                                        /* Use dual_ instead of array
4357                                           Even though dual_ is only numberRows_ long this is
4358                                           okay as gets permuted to longer rowArray_[2]
4359                                        */
4360                                        arrayVector->setDenseVector(dual_);
4361                                        int * index = arrayVector->getIndices();
4362                                        int number = 0;
4363                                        for (iRow = 0; iRow < numberRows_; iRow++) {
4364                                             int iPivot = pivotVariable_[iRow];
4365                                             double value = cost_[iPivot];
4366                                             dual_[iRow] = value;
4367                                             if (value) {
4368                                                  index[number++] = iRow;
4369                                             }
4370                                        }
4371                                        arrayVector->setNumElements(number);
4372                                        // Extended duals before "updateTranspose"
4373                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
4374                                        // Btran basic costs
4375                                        rowArray_[2]->clear();
4376                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
4377                                        // and return vector
4378                                        arrayVector->setDenseVector(array);
4379                                   }
4380                              }
4381                              problemStatus_ = 0;
4382                              sumPrimalInfeasibilities_ = 0.0;
4383                              if ((specialOptions_&(1024 + 16384)) != 0) {
4384                                   CoinIndexedVector * arrayVector = rowArray_[1];
4385                                   arrayVector->clear();
4386                                   double * rhs = arrayVector->denseVector();
4387                                   times(1.0, solution_, rhs);
4388                                   bool bad2 = false;
4389                                   int i;
4390                                   for ( i = 0; i < numberRows_; i++) {
4391                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
4392                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
4393                                             bad2 = true;
4394                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
4395                                        }
4396                                        rhs[i] = 0.0;
4397                                   }
4398                                   for ( i = 0; i < numberColumns_; i++) {
4399                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
4400                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4401                                             bad2 = true;
4402                                        }
4403                                   }
4404                                   if (bad2) {
4405                                        problemStatus_ = -3;
4406                                        returnCode = -2;
4407                                        // Force to re-factorize early next time
4408                                        int numberPivots = factorization_->pivots();
4409                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4410                                   }
4411                              }
4412                         }
4413                    }
4414               } else {
4415                    problemStatus_ = -3;
4416                    returnCode = -2;
4417                    // Force to re-factorize early next time
4418                    int numberPivots = factorization_->pivots();
4419                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4420               }
4421               break;
4422          }
4423     }
4424     startingTheta = lastTheta+theta_;
4425     return returnCode;
4426}
4427#if 0
4428static int zzzzzz=0;
4429int zzzzzzOther=0;
4430#endif
4431// Finds best possible pivot
4432double 
4433ClpSimplexOther::bestPivot(bool justColumns)
4434{
4435  // Get good size for pivot
4436  // Allow first few iterations to take tiny
4437  double acceptablePivot = 1.0e-9;
4438  if (numberIterations_ > 100)
4439    acceptablePivot = 1.0e-8;
4440  if (factorization_->pivots() > 10 ||
4441      (factorization_->pivots() && sumDualInfeasibilities_))
4442    acceptablePivot = 1.0e-5; // if we have iterated be more strict
4443  else if (factorization_->pivots() > 5)
4444    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4445  else if (factorization_->pivots())
4446    acceptablePivot = 1.0e-8; // relax
4447  double bestPossiblePivot = 1.0;
4448  // get sign for finding row of tableau
4449  // normal iteration
4450  // create as packed
4451  double direction = directionOut_;
4452#ifndef COIN_FAC_NEW
4453  rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4454#else
4455  rowArray_[0]->createOneUnpackedElement(pivotRow_, direction);
4456#endif
4457  factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4458  // put row of tableau in rowArray[0] and columnArray[0]
4459  matrix_->transposeTimes(this, -1.0,
4460                          rowArray_[0], rowArray_[3], columnArray_[0]);
4461  sequenceIn_=-1;
4462  if (justColumns)
4463    rowArray_[0]->clear();
4464  // do ratio test for normal iteration
4465  bestPossiblePivot = 
4466    reinterpret_cast<ClpSimplexDual *> 
4467    ( this)->dualColumn(rowArray_[0],
4468                        columnArray_[0], columnArray_[1],
4469                        rowArray_[3], acceptablePivot, NULL);
4470  return bestPossiblePivot;
4471}
4472// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4473int
4474ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData,
4475                           const double * /*changeObjective*/)
4476{
4477  const double * lowerChange = paramData.lowerChange;
4478  const double * upperChange = paramData.upperChange;
4479  const int * lowerList = paramData.lowerList;
4480  const int * upperList = paramData.upperList;
4481  int iSequence;
4482  bool toLower = false;
4483  //assert (type==1);
4484  // may need to decide based on model?
4485  bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4486  double * array = rowArray_[4]->denseVector();
4487  //rowArray_[4]->checkClean();
4488  const int * row = matrix_->getIndices();
4489  const int * columnLength = matrix_->getVectorLengths();
4490  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4491  const double * elementByColumn = matrix_->getElements();
4492#if 0
4493  double tempArray[5000];
4494  bool checkIt=false;
4495  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
4496    memcpy(tempArray,array,numberRows_*sizeof(double));
4497    checkIt=true;
4498    needFullUpdate=true;
4499  }
4500#endif
4501#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4502  double * lowerGap = paramData.lowerGap;
4503  double * upperGap = paramData.upperGap;
4504  double * lowerCoefficient = paramData.lowerCoefficient;
4505  double * upperCoefficient = paramData.upperCoefficient;
4506  int * lowerActive=paramData.lowerActive;
4507  int * upperActive=paramData.upperActive;
4508#endif
4509  if (!factorization_->pivots()||needFullUpdate) {
4510    //zzzzzz=0;
4511    rowArray_[4]->clear();
4512    // get change
4513    if (!rowScale_) {
4514      int n;
4515      n=lowerList[-2];
4516      int i;
4517      for (i=0;i<n;i++) {
4518        int iSequence = lowerList[i];
4519        assert (iSequence<numberColumns_);
4520        if (getColumnStatus(iSequence)==atLowerBound) {
4521          double value=lowerChange[iSequence];
4522          for (CoinBigIndex j = columnStart[iSequence];
4523               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4524            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4525          }
4526        }
4527      }
4528      n=lowerList[-1];
4529      const double * change = lowerChange+numberColumns_;
4530      for (;i<n;i++) {
4531        int iSequence = lowerList[i]-numberColumns_;
4532        assert (iSequence>=0);
4533        if (getRowStatus(iSequence)==atLowerBound) {
4534          double value=change[iSequence];
4535          rowArray_[4]->quickAdd(iSequence, -value);
4536        }
4537      }
4538      n=upperList[-2];
4539      for (i=0;i<n;i++) {
4540        int iSequence = upperList[i];
4541        assert (iSequence<numberColumns_);
4542        if (getColumnStatus(iSequence)==atUpperBound) {
4543          double value=upperChange[iSequence];
4544          for (CoinBigIndex j = columnStart[iSequence];
4545               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4546            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4547          }
4548        }
4549      }
4550      n=upperList[-1];
4551      change = upperChange+numberColumns_;
4552      for (;i<n;i++) {
4553        int iSequence = upperList[i]-numberColumns_;
4554        assert (iSequence>=0);
4555        if (getRowStatus(iSequence)==atUpperBound) {
4556          double value=change[iSequence];
4557          rowArray_[4]->quickAdd(iSequence, -value);
4558        }
4559      }
4560    } else {
4561      int n;
4562      n=lowerList[-2];
4563      int i;
4564      for (i=0;i<n;i++) {
4565        int iSequence = lowerList[i];
4566        assert (iSequence<numberColumns_);
4567        if (getColumnStatus(iSequence)==atLowerBound) {
4568          double value=lowerChange[iSequence];
4569          // apply scaling
4570          double scale = columnScale_[iSequence];
4571          for (CoinBigIndex j = columnStart[iSequence];
4572               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4573            int iRow = row[j];
4574            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4575          }
4576        }
4577      }
4578      n=lowerList[-1];
4579      const double * change = lowerChange+numberColumns_;
4580      for (;i<n;i++) {
4581        int iSequence = lowerList[i]-numberColumns_;
4582        assert (iSequence>=0);
4583        if (getRowStatus(iSequence)==atLowerBound) {
4584          double value=change[iSequence];
4585          rowArray_[4]->quickAdd(iSequence, -value);
4586        }
4587      }
4588      n=upperList[-2];
4589      for (i=0;i<n;i++) {
4590        int iSequence = upperList[i];
4591        assert (iSequence<numberColumns_);
4592        if (getColumnStatus(iSequence)==atUpperBound) {
4593          double value=upperChange[iSequence];
4594          // apply scaling
4595          double scale = columnScale_[iSequence];
4596          for (CoinBigIndex j = columnStart[iSequence];
4597               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4598            int iRow = row[j];
4599            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4600          }
4601        }
4602      }
4603      n=upperList[-1];
4604      change = upperChange+numberColumns_;
4605      for (;i<n;i++) {
4606        int iSequence = upperList[i]-numberColumns_;
4607        assert (iSequence>=0);
4608        if (getRowStatus(iSequence)==atUpperBound) {
4609          double value=change[iSequence];
4610          rowArray_[4]->quickAdd(iSequence, -value);
4611        }
4612      }
4613    }
4614    // ftran it
4615    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4616#if 0
4617    if (checkIt) {
4618      for (int i=0;i<numberRows_;i++) {
4619        assert (fabs(tempArray[i]-array[i])<1.0e-8);
4620      }
4621    }
4622#endif
4623#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4624    /* later for sparse - keep like CoinIndexedvector
4625       and just redo here */
4626    int lowerN=0;
4627    int upperN=0;
4628    memset(lowerCoefficient,0,numberRows_*sizeof(double));
4629    memset(upperCoefficient,0,numberRows_*sizeof(double));
4630    for (int iRow=0;iRow<numberRows_;iRow++) {
4631      iSequence = pivotVariable_[iRow];
4632      double currentSolution = solution_[iSequence];
4633      double alpha = array[iRow];
4634      double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4635      double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4636      if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4637        double currentLower = lower_[iSequence];
4638        ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4639        double gap=currentSolution-currentLower;
4640        lowerGap[iRow]=gap;
4641        lowerCoefficient[iRow]=thetaCoefficientLower;
4642        lowerActive[lowerN++]=iRow;
4643        //} else {
4644        //lowerCoefficient[iRow]=0.0;
4645      }
4646      if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4647        double currentUpper = upper_[iSequence];
4648        ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4649        double gap2=-(currentSolution-currentUpper); //positive
4650        upperGap[iRow]=gap2;
4651        upperCoefficient[iRow]=-thetaCoefficientUpper;
4652        upperActive[upperN++]=iRow;
4653      }
4654    }
4655    assert (lowerN>=0&&lowerN<=numberRows_);
4656    lowerActive[-1]=lowerN;
4657    upperActive[-1]=upperN;
4658#endif
4659  } else if (sequenceIn_>=0) {
4660    //assert (sequenceIn_>=0);
4661    assert (sequenceOut_>=0);
4662    assert (sequenceIn_!=sequenceOut_);
4663    double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4664    int needed=0;
4665    assert (!rowArray_[5]->getNumElements());
4666    if (change) {
4667      if (sequenceIn_<numberColumns_) {
4668        if (!rowScale_) {
4669          for (CoinBigIndex i = columnStart[sequenceIn_];
4670               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4671            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4672          }
4673        } else {
4674          // apply scaling
4675          double scale = columnScale_[sequenceIn_];
4676          for (CoinBigIndex i = columnStart[sequenceIn_];
4677               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4678            int iRow = row[i];
4679            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4680          }
4681        }
4682      } else {
4683        rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4684      }
4685      needed++;
4686    }
4687    if (getStatus(sequenceOut_)==atLowerBound)
4688      change=lowerChange[sequenceOut_];
4689    else
4690      change=upperChange[sequenceOut_];
4691    if (change) {
4692      if (sequenceOut_<numberColumns_) {
4693        if (!rowScale_) {
4694          for (CoinBigIndex i = columnStart[sequenceOut_];
4695               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4696            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4697          }
4698        } else {
4699          // apply scaling
4700          double scale = columnScale_[sequenceOut_];
4701          for (CoinBigIndex i = columnStart[sequenceOut_];
4702               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4703            int iRow = row[i];
4704            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4705          }
4706        }
4707      } else {
4708        rowArray_[5]->quickAdd(sequenceOut_-numberColumns_,-change);
4709      }
4710      needed++;
4711    }
4712    //printf("seqin %d seqout %d needed %d\n",
4713    //     sequenceIn_,sequenceOut_,needed);
4714    if (needed) {
4715      // ftran it
4716      factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4717      // add
4718      double * array5 = rowArray_[5]->denseVector();
4719      int * index5 = rowArray_[5]->getIndices();
4720      int number5 = rowArray_[5]->getNumElements();
4721#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4722      int lowerN=lowerActive[-1];
4723      int upperN=upperActive[-1];
4724      int nIn4=rowArray_[4]->getNumElements();
4725      int * index4 = rowArray_[4]->getIndices();
4726#endif
4727      for (int i = 0; i < number5; i++) {
4728        int iPivot = index5[i];
4729#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4730        rowArray_[4]->quickAdd(iPivot,array5[iPivot]);
4731#else
4732        /* later for sparse - modify here */
4733        int iSequence = pivotVariable_[iPivot];
4734        double currentSolution = solution_[iSequence];
4735        double currentAlpha = array[iPivot];
4736        double alpha5 = array5[iPivot];
4737        double alpha = currentAlpha+alpha5;
4738        if (currentAlpha) {
4739          if (alpha) {
4740            array[iPivot] = alpha;
4741          } else {
4742            array[iPivot] = COIN_DBL_MIN;
4743          }
4744        } else {
4745          index4[nIn4++] = iPivot;
4746          array[iPivot] = alpha;
4747        }
4748        double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4749        double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4750        double oldLower = lowerCoefficient[iPivot];
4751        double oldUpper = upperCoefficient[iPivot];
4752        if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4753          double currentLower = lower_[iSequence];
4754          ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4755          double gap=currentSolution-currentLower;
4756          lowerGap[iPivot]=gap;
4757          lowerCoefficient[iPivot]=thetaCoefficientLower;
4758          if (!oldLower)
4759            lowerActive[lowerN++]=iPivot;
4760        } else {
4761          if (oldLower)
4762            lowerCoefficient[iPivot]=COIN_DBL_MIN;
4763        }
4764        if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4765          double currentUpper = upper_[iSequence];
4766          ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4767          double gap2=-(currentSolution-currentUpper); //positive
4768          upperGap[iPivot]=gap2;
4769          upperCoefficient[iPivot]=-thetaCoefficientUpper;
4770          if (!oldUpper)
4771            upperActive[upperN++]=iPivot;
4772        } else {
4773          if (oldUpper)
4774            upperCoefficient[iPivot]=COIN_DBL_MIN;
4775        }
4776#endif
4777        array5[iPivot]=0.0;
4778      }
4779      rowArray_[5]->setNumElements(0);
4780#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4781      rowArray_[4]->setNumElements(nIn4);
4782      assert (lowerN>=0&&lowerN<=numberRows_);
4783      lowerActive[-1]=lowerN;
4784      upperActive[-1]=upperN;
4785#endif
4786    }
4787  }
4788  const int * index = rowArray_[4]->getIndices();
4789  int number = rowArray_[4]->getNumElements();
4790#define TESTXX 0
4791#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4792  int * markDone = reinterpret_cast<int *>(paramData.markDone);
4793  int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT;
4794  memset(markDone,0,nToZero*sizeof(int));
4795  const int * backwardBasic = paramData.backwardBasic;
4796#endif
4797  // first ones with alpha
4798  double theta1=maxTheta;
4799  int pivotRow1=-1;
4800#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4801  int pivotRow2=-1;
4802  double theta2=maxTheta;
4803#endif
4804#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4805  for (int i=0;i<number;i++) {
4806    int iPivot=index[i];
4807    iSequence = pivotVariable_[iPivot];
4808    //assert(!markDone[iSequence]);
4809    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4810    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4811    markDone[word] |= ( 1 << bit );
4812    // solution value will be sol - theta*alpha
4813    // bounds will be bounds + change *theta
4814    double currentSolution = solution_[iSequence];
4815    double alpha = array[iPivot];
4816    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4817    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4818    if (thetaCoefficientLower > 1.0e-8) {
4819      double currentLower = lower_[iSequence];
4820      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4821      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4822      double gap=currentSolution-currentLower;
4823      if (thetaCoefficientLower*theta1>gap) {
4824        theta1 = gap/thetaCoefficientLower;
4825        //toLower=true;
4826        pivotRow1=iPivot;
4827      }
4828    }
4829    if (thetaCoefficientUpper < -1.0e-8) {
4830      double currentUpper = upper_[iSequence];
4831      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4832      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4833      double gap2=currentSolution-currentUpper; //negative
4834      if (thetaCoefficientUpper*theta2<gap2) {
4835        theta2 = gap2/thetaCoefficientUpper;
4836        //toLower=false;
4837        pivotRow2=iPivot;
4838      }
4839    }
4840  }
4841  // now others
4842  int nLook=lowerList[-1];
4843  for (int i=0;i<nLook;i++) {
4844    int iSequence = lowerList[i];
4845    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4846    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4847    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4848      double currentSolution = solution_[iSequence];
4849      double currentLower = lower_[iSequence];
4850      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4851      double thetaCoefficient = lowerChange[iSequence];
4852      if (thetaCoefficient > 0.0) {
4853        double gap=currentSolution-currentLower;
4854        if (thetaCoefficient*theta1>gap) {
4855          theta1 = gap/thetaCoefficient;
4856          //toLower=true;
4857          pivotRow1 = backwardBasic[iSequence];
4858        }
4859      }
4860    }
4861  }
4862  nLook=upperList[-1];
4863  for (int i=0;i<nLook;i++) {
4864    int iSequence = upperList[i];
4865    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4866    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4867    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4868      double currentSolution = solution_[iSequence];
4869      double currentUpper = upper_[iSequence];
4870      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4871      double thetaCoefficient = upperChange[iSequence];
4872      if (thetaCoefficient < 0) {
4873        double gap=currentSolution-currentUpper; //negative
4874        if (thetaCoefficient*theta2<gap) {
4875          theta2 = gap/thetaCoefficient;
4876          //toLower=false;
4877          pivotRow2 = backwardBasic[iSequence];
4878        }
4879      }
4880    }
4881  }
4882  if (theta2<theta1) {
4883    theta_=theta2;
4884    toLower=false;
4885    pivotRow_=pivotRow2;
4886  } else {
4887    theta_=theta1;
4888    toLower=true;
4889    pivotRow_=pivotRow1;
4890  }
4891#if 0 //TESTXX
4892#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4893  {
4894    double * checkArray = new double[numberRows_];
4895    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4896    int lowerN=lowerActive[-1];
4897    for (int i=0;i<lowerN;i++) {
4898      int iRow=lowerActive[i];
4899      int iSequence = pivotVariable_[iRow];
4900      double alpha = array[iRow];
4901      double thetaCoefficient = lowerChange[iSequence] + alpha;
4902      if (thetaCoefficient > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4903        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4904        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4905          abort();
4906        }
4907      } else {
4908        assert (fabs(checkArray[iRow])<1.0e-12);
4909        if (fabs(checkArray[iRow])>1.0e-12) {
4910          abort();
4911        }
4912      }
4913      checkArray[iRow]=0.0;
4914    }
4915    for (int i=0;i<numberRows_;i++) {
4916      assert (!checkArray[i]);
4917      if (checkArray[i])
4918        abort();
4919    }
4920    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
4921    int upperN=upperActive[-1];
4922    for (int i=0;i<upperN;i++) {
4923      int iRow=upperActive[i];
4924      int iSequence = pivotVariable_[iRow];
4925      double alpha = array[iRow];
4926      double thetaCoefficient = -(upperChange[iSequence] + alpha);
4927      if (thetaCoefficient > 1.0e-8&&upper_[iSequence]<1.0e30) {
4928        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4929        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4930          abort();
4931        }
4932      } else {
4933        assert (fabs(checkArray[iRow])<1.0e-12);
4934        if (fabs(checkArray[iRow])>1.0e-12) {
4935          abort();
4936        }
4937      }
4938      checkArray[iRow]=0.0;
4939    }
4940    for (int i=0;i<numberRows_;i++) {
4941      assert (!checkArray[i]);
4942      if (checkArray[i])
4943        abort();
4944    }
4945    delete [] checkArray;
4946  }
4947  double theta3=maxTheta;
4948  int pivotRow3=-1;
4949  int lowerN=lowerActive[-1];
4950  for (int i=0;i<lowerN;i++) {
4951    int iRow=lowerActive[i];
4952    double lowerC = lowerCoefficient[iRow];
4953    double gap=lowerGap[iRow];
4954    if (toLower&&iRow==pivotRow_) {
4955      assert (lowerC*theta3>gap-1.0e-8);
4956      if (lowerC*theta3<gap-1.0e-8)
4957        abort();
4958    }
4959    if (lowerC*theta3>gap&&lowerC!=COIN_DBL_MIN) {
4960      theta3 = gap/lowerC;
4961      pivotRow3=iRow;
4962    }
4963  }
4964  int pivotRow4=pivotRow3;
4965  double theta4=theta3;
4966  int upperN=upperActive[-1];
4967  for (int i=0;i<upperN;i++) {
4968    int iRow=upperActive[i];
4969    double upperC = upperCoefficient[iRow];
4970    double gap=upperGap[iRow];
4971    if (!toLower&&iRow==pivotRow_) {
4972      assert (upperC*theta3>gap-1.0e-8);
4973      if (upperC*theta3<gap-1.0e-8)
4974        abort();
4975    }
4976    if (upperC*theta4>gap&&upperC!=COIN_DBL_MIN) {
4977      theta4 = gap/upperC;
4978      pivotRow4=iRow;
4979    }
4980  }
4981  bool toLower3;
4982  if (theta4<theta3) {
4983    theta3=theta4;
4984    toLower3=false;
4985    pivotRow3=pivotRow4;
4986  } else {
4987    toLower3=true;
4988  }
4989  if (fabs(theta3-theta_)>1.0e-8)
4990    abort();
4991  if (toLower!=toLower3||pivotRow_!=pivotRow3) {
4992    printf("bad piv - good %d %g %s, bad %d %g %s\n",pivotRow_,theta_,toLower ? "toLower" : "toUpper",
4993           pivotRow3,theta3,toLower3 ? "toLower" : "toUpper");
4994    //zzzzzz++;
4995    if (true/*zzzzzz>zzzzzzOther*/) {
4996      printf("Swapping\n");
4997      pivotRow_=pivotRow3;
4998      theta_=theta3;
4999      toLower=toLower3;
5000    }
5001  }
5002#endif
5003#endif
5004#else
5005#if 0 //CLP_PARAMETRIC_DENSE_ARRAYS==2
5006  {
5007    double * checkArray = new double[numberRows_];
5008    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
5009    int lowerN=lowerActive[-1];
5010    for (int i=0;i<lowerN;i++) {
5011      int iRow=lowerActive[i];
5012      checkArray[iRow]=0.0;
5013    }
5014    for (int i=0;i<numberRows_;i++) {
5015      assert (!checkArray[i]);
5016      if (checkArray[i])
5017        abort();
5018    }
5019    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
5020    int upperN=upperActive[-1];
5021    for (int i=0;i<upperN;i++) {
5022      int iRow=upperActive[i];
5023      checkArray[iRow]=0.0;
5024    }
5025    for (int i=0;i<numberRows_;i++) {
5026      assert (!checkArray[i]);
5027      if (checkArray[i])
5028        abort();
5029    }
5030    delete [] checkArray;
5031  }
5032#endif
5033  int lowerN=lowerActive[-1];
5034  for (int i=0;i<lowerN;i++) {
5035    int iRow=lowerActive[i];
5036    double lowerC = lowerCoefficient[iRow];
5037    double gap=lowerGap[iRow];
5038    if (lowerC*theta1>gap&&lowerC!=COIN_DBL_MIN) {
5039      theta1 = gap/lowerC;
5040      pivotRow1=iRow;
5041    }
5042  }
5043  pivotRow_=pivotRow1;
5044  theta_=theta1;
5045  int upperN=upperActive[-1];
5046  for (int i=0;i<upperN;i++) {
5047    int iRow=upperActive[i];
5048    double upperC = upperCoefficient[iRow];
5049    double gap=upperGap[iRow];
5050    if (upperC*theta1>gap&&upperC!=COIN_DBL_MIN) {
5051      theta1 = gap/upperC;
5052      pivotRow1=iRow;
5053    }
5054  }
5055  if (theta1<theta_) {
5056    theta_=theta1;
5057    toLower=false;
5058    pivotRow_=pivotRow1;
5059  } else {
5060    toLower=true;
5061  }
5062#endif
5063  theta_ = CoinMax(theta_,0.0);
5064  if (theta_>1.0e-15) {
5065    // update solution
5066    for (int iRow = 0; iRow < number; iRow++) {
5067      int iPivot = index[iRow];
5068      iSequence = pivotVariable_[iPivot];
5069      // solution value will be sol - theta*alpha
5070      double alpha = array[iPivot];
5071      double currentSolution = solution_[iSequence] - theta_ * alpha;
5072      solution_[iSequence] =currentSolution;
5073#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5074      if (lower_[iSequence]>-1.0e30)
5075        lowerGap[iPivot]=currentSolution-lower_[iSequence];
5076      if (upper_[iSequence]<1.0e30)
5077        upperGap[iPivot]=-(currentSolution-upper_[iSequence]);
5078#endif
5079    }
5080  }
5081#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5082  if (pivotRow_>=0&&false) {
5083    double oldValue = upperCoefficient[pivotRow_];
5084    double value = array[pivotRow_];
5085    if (value) {
5086      if (!oldValue) {
5087        int upperN=upperActive[-1];
5088        assert (upperN>=0&&upperN<numberRows_);
5089        upperActive[upperN]=pivotRow_;
5090        upperActive[-1]=upperN+1;
5091      }
5092    } else {
5093      if (oldValue)
5094        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
5095    }
5096  }
5097#endif
5098#if 0
5099  for (int i=0;i<numberTotal;i++)
5100    assert(!markDone[i]);
5101#endif
5102  if (pivotRow_ >= 0) {
5103    sequenceOut_ = pivotVariable_[pivotRow_];
5104    valueOut_ = solution_[sequenceOut_];
5105    lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
5106    upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
5107    if (!toLower) {
5108      directionOut_ = -1;
5109      dualOut_ = valueOut_ - upperOut_;
5110    } else {
5111      directionOut_ = 1;
5112      dualOut_ = lowerOut_ - valueOut_;
5113    }
5114    return 0;
5115  } else {
5116    //theta_=0.0;
5117    return -1;
5118  }
5119}
5120// Restores bound to original bound
5121void 
5122ClpSimplexOther::originalBound(int iSequence, double theta, 
5123                               const double * lowerChange,
5124                               const double * upperChange)
5125{
5126     if (getFakeBound(iSequence) != noFake) {
5127          numberFake_--;
5128          setFakeBound(iSequence, noFake);
5129          if (iSequence >= numberColumns_) {
5130               // rows
5131               int iRow = iSequence - numberColumns_;
5132               rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
5133               rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
5134               if (rowScale_) {
5135                    if (rowLowerWork_[iRow] > -1.0e50)
5136                         rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5137                    if (rowUpperWork_[iRow] < 1.0e50)
5138                         rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5139               } else if (rhsScale_ != 1.0) {
5140                    if (rowLowerWork_[iRow] > -1.0e50)
5141                         rowLowerWork_[iRow] *= rhsScale_;
5142                    if (rowUpperWork_[iRow] < 1.0e50)
5143                         rowUpperWork_[iRow] *= rhsScale_;
5144               }
5145          } else {
5146               // columns
5147               columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
5148               columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
5149               if (rowScale_) {
5150                    double multiplier = 1.0 * inverseColumnScale_[iSequence];
5151                    if (columnLowerWork_[iSequence] > -1.0e50)
5152                         columnLowerWork_[iSequence] *= multiplier * rhsScale_;
5153                    if (columnUpperWork_[iSequence] < 1.0e50)
5154                         columnUpperWork_[iSequence] *= multiplier * rhsScale_;
5155               } else if (rhsScale_ != 1.0) {
5156                    if (columnLowerWork_[iSequence] > -1.0e50)
5157                         columnLowerWork_[iSequence] *= rhsScale_;
5158                    if (columnUpperWork_[iSequence] < 1.0e50)
5159                         columnUpperWork_[iSequence] *= rhsScale_;
5160               }
5161          }
5162     }
5163}
5164/* Expands out all possible combinations for a knapsack
5165   If buildObj NULL then just computes space needed - returns number elements
5166   On entry numberOutput is maximum allowed, on exit it is number needed or
5167   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
5168   least space to return values which reconstruct input.
5169   Rows returned will be original rows but no entries will be returned for
5170   any rows all of whose entries are in knapsack.  So up to user to allow for this.
5171   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
5172   in expanded knapsack.  Values in buildRow and buildElement;
5173*/
5174int
5175ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
5176                                double * buildObj, CoinBigIndex * buildStart,
5177                                int * buildRow, double * buildElement, int reConstruct) const
5178{
5179     int iRow;
5180     int iColumn;
5181     // Get column copy
5182     CoinPackedMatrix * columnCopy = matrix();
5183     // Get a row copy in standard format
5184     CoinPackedMatrix matrixByRow;
5185     matrixByRow.reverseOrderedCopyOf(*columnCopy);
5186     const double * elementByRow = matrixByRow.getElements();
5187     const int * column = matrixByRow.getIndices();
5188     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5189     const int * rowLength = matrixByRow.getVectorLengths();
5190     CoinBigIndex j;
5191     int * whichColumn = new int [numberColumns_];
5192     int * whichRow = new int [numberRows_];
5193     int numJ = 0;
5194     // Get what other columns can compensate for
5195     double * lo = new double [numberRows_];
5196     double * high = new double [numberRows_];
5197     {
5198          // Use to get tight column bounds
5199          ClpSimplex tempModel(*this);
5200          tempModel.tightenPrimalBounds(0.0, 0, true);
5201          // Now another model without knapsacks
5202          int nCol = 0;
5203          for (iRow = 0; iRow < numberRows_; iRow++) {
5204               whichRow[iRow] = iRow;
5205          }
5206          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
5207               whichColumn[iColumn] = -1;
5208          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
5209               int iColumn = column[j];
5210               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
5211                    whichColumn[iColumn] = 0;
5212               } else {
5213                    assert (!columnLower_[iColumn]); // fix later
5214               }
5215          }
5216          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5217               if (whichColumn[iColumn] < 0)
5218                    whichColumn[nCol++] = iColumn;
5219          }
5220          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
5221          // Row copy
5222          CoinPackedMatrix matrixByRow;
5223          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
5224          const double * elementByRow = matrixByRow.getElements();
5225          const int * column = matrixByRow.getIndices();
5226          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5227          const int * rowLength = matrixByRow.getVectorLengths();
5228          const double * columnLower = tempModel2.getColLower();
5229          const double * columnUpper = tempModel2.getColUpper();
5230          for (iRow = 0; iRow < numberRows_; iRow++) {
5231               lo[iRow] = -COIN_DBL_MAX;
5232               high[iRow] = COIN_DBL_MAX;
5233               if (rowLower_[iRow] >