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

Last change on this file since 1861 was 1861, checked in by forrest, 8 years ago

allow user to override crossover considerations in parametrics

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 355.0 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1861 2012-06-08 10:16:41Z 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    if (lower > upper) {
2978      maxTheta = -1.0;
2979      break;
2980    }
2981    rowLower_[iRow]=lower;
2982    rowUpper_[iRow]=upper;
2983  }
2984  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2985    double lower = columnLower_[iColumn];
2986    double upper = columnUpper_[iColumn];
2987    if (lower<-1.0e30)
2988      lowerChange[iColumn]=0.0;
2989    double chgLower = lowerChange[iColumn];
2990    largestChange=CoinMax(largestChange,fabs(chgLower));
2991    if (upper>1.0e30)
2992      upperChange[iColumn]=0.0;
2993    double chgUpper = upperChange[iColumn];
2994    largestChange=CoinMax(largestChange,fabs(chgUpper));
2995    if (lower > -1.0e30 && upper < 1.0e30) {
2996      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2997        maxTheta = (upper - lower) / (chgLower - chgUpper);
2998      }
2999    }
3000    lower+=startingTheta*chgLower;
3001    upper+=startingTheta*chgUpper;
3002    if (lower > upper) {
3003      maxTheta = -1.0;
3004      break;
3005    }
3006    columnLower_[iColumn]=lower;
3007    columnUpper_[iColumn]=upper;
3008  }
3009  if (maxTheta == 1.0e50)
3010    maxTheta = COIN_DBL_MAX;
3011  int returnCode=0;
3012  if (maxTheta < 0.0) {
3013    // bad ranges or initial
3014    returnCode = -1;
3015  }
3016  if (maxTheta < endingTheta) {
3017    char line[100];
3018    sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
3019            endingTheta,maxTheta);
3020    handler_->message(CLP_GENERAL,messages_)
3021      << line << CoinMessageEol;
3022    endingTheta = maxTheta;
3023  }
3024  if (endingTheta < startingTheta) {
3025    // bad initial
3026    returnCode = -2;
3027  }
3028#ifdef CLP_USER_DRIVEN
3029  // accept user's version
3030  paramData.maxTheta=endingTheta;
3031#else
3032  paramData.maxTheta=CoinMin(maxTheta,endingTheta);
3033#endif
3034  /* given largest change element choose acceptable end
3035     be safe and make sure difference < 0.1*tolerance */
3036  double acceptableDifference=0.1*primalTolerance_/
3037    CoinMax(largestChange,1.0);
3038  paramData.acceptableMaxTheta=maxTheta-acceptableDifference;
3039  bool swapped=false;
3040  // Dantzig
3041#define ALL_DANTZIG
3042#ifdef ALL_DANTZIG
3043  ClpDualRowPivot * savePivot = dualRowPivot_;
3044  dualRowPivot_ = new ClpDualRowDantzig();
3045  dualRowPivot_->setModel(this);
3046#else
3047  ClpDualRowPivot * savePivot = NULL;
3048#endif
3049  if (!returnCode) {
3050    assert (objective_->type()==1);
3051    objective_->setType(2); // in case matrix empty
3052    returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
3053    objective_->setType(1);
3054    if (!returnCode) {
3055      double saveDualBound=dualBound_;
3056      dualBound_=CoinMax(dualBound_,1.0e15);
3057      swapped=true;
3058      double * temp;
3059      memcpy(saveLower,lower_,numberTotal*sizeof(double));
3060      temp=saveLower;
3061      saveLower=lower_;
3062      lower_=temp;
3063      //columnLowerWork_ = lower_;
3064      //rowLowerWork_ = lower_ + numberColumns_;
3065      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
3066      temp=saveUpper;
3067      saveUpper=upper_;
3068      upper_=temp;
3069      //columnUpperWork_ = upper_;
3070      //rowUpperWork_ = upper_ + numberColumns_;
3071      if (rowScale_) {
3072        // scale saved and change arrays
3073        double * lowerChange = lower_+numberTotal;
3074        double * upperChange = upper_+numberTotal;
3075        double * lowerSave = lowerChange+numberTotal;
3076        double * upperSave = upperChange+numberTotal;
3077        for (int i=0;i<numberColumns_;i++) {
3078          double multiplier = inverseColumnScale_[i];
3079          if (lowerSave[i]>-1.0e20)
3080            lowerSave[i] *= multiplier;
3081          if (upperSave[i]<1.0e20)
3082            upperSave[i] *= multiplier;
3083          lowerChange[i] *= multiplier;
3084          upperChange[i] *= multiplier;
3085        }
3086        lowerChange += numberColumns_;
3087        upperChange += numberColumns_;
3088        lowerSave += numberColumns_;
3089        upperSave += numberColumns_;
3090        for (int i=0;i<numberRows_;i++) {
3091          double multiplier = rowScale_[i];
3092          if (lowerSave[i]>-1.0e20)
3093            lowerSave[i] *= multiplier;
3094          if (upperSave[i]<1.0e20)
3095            upperSave[i] *= multiplier;
3096          lowerChange[i] *= multiplier;
3097          upperChange[i] *= multiplier;
3098        }
3099      }
3100      //double saveEndingTheta = endingTheta;
3101      double * saveDuals = NULL;
3102      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3103      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3104        // probably can get rid of this if we adjust every change in theta
3105        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3106        //   sumPrimalInfeasibilities_);
3107        int pass=100;
3108        while(sumPrimalInfeasibilities_) {
3109          pass--;
3110          if (!pass)
3111            break;
3112          problemStatus_=-1;
3113          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3114            double value=solution_[iSequence];
3115            // remember scaling
3116            if (value<lower_[iSequence]-1.0e-9) {
3117              lower_[iSequence]=value;
3118              lowerCopy[iSequence]=value;
3119            } else if (value>upper_[iSequence]+1.0e-9) {
3120              upper_[iSequence]=value;
3121              upperCopy[iSequence]=value;
3122            }
3123          }
3124          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3125        }
3126      }
3127      if (!problemStatus_) {
3128        if (nLowerChange||nUpperChange) {
3129#ifndef ALL_DANTZIG
3130          // Dantzig
3131          savePivot = dualRowPivot_;
3132          dualRowPivot_ = new ClpDualRowDantzig();
3133          dualRowPivot_->setModel(this);
3134#endif
3135          //for (int i=0;i<numberRows_+numberColumns_;i++)
3136          //setFakeBound(i, noFake);
3137          // Now do parametrics
3138          handler_->message(CLP_PARAMETRICS_STATS, messages_)
3139            << startingTheta << objectiveValue() << CoinMessageEol;
3140          bool canSkipFactorization=true;
3141          while (!returnCode) {
3142            paramData.startingTheta=startingTheta;
3143            paramData.endingTheta=endingTheta;
3144            returnCode = parametricsLoop(paramData,
3145                                         data,canSkipFactorization);
3146            startingTheta=paramData.startingTheta;
3147            endingTheta=paramData.endingTheta;
3148            canSkipFactorization=false;
3149            if (!returnCode) {
3150              //startingTheta = endingTheta;
3151              //endingTheta = saveEndingTheta;
3152              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3153                << startingTheta << objectiveValue() << CoinMessageEol;
3154              if (startingTheta >= endingTheta-primalTolerance_
3155                  ||problemStatus_==2)
3156                break;
3157            } else if (returnCode == -1) {
3158              // trouble - do external solve
3159              abort(); //needToDoSomething = true;
3160            } else if (problemStatus_==1) {
3161              // can't move any further
3162              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3163                << endingTheta << objectiveValue() << CoinMessageEol;
3164              problemStatus_=0;
3165            }
3166          }
3167        }
3168        dualBound_ = saveDualBound;
3169        //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3170      } else {
3171        // check if empty
3172        //if (!numberRows_||!matrix_->getNumElements()) {
3173        // success
3174#ifdef CLP_USER_DRIVEN
3175        //theta_ = endingTheta;
3176        //eventHandler_->event(ClpEventHandler::theta);
3177#endif
3178        //}
3179      }
3180    }
3181    if (problemStatus_==2) {
3182      delete [] ray_;
3183      ray_ = new double [numberColumns_];
3184    }
3185    if (swapped&&lower_) {
3186      double * temp=saveLower;
3187      saveLower=lower_;
3188      lower_=temp;
3189      temp=saveUpper;
3190      saveUpper=upper_;
3191      upper_=temp;
3192    }
3193    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
3194  }   
3195  if (!scalingFlag_) {
3196    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3197    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3198    memcpy(rowLower_,lowerCopy+numberColumns_,
3199           numberRows_*sizeof(double));
3200    memcpy(rowUpper_,upperCopy+numberColumns_,
3201           numberRows_*sizeof(double));
3202  } else {
3203    //  extra for unscaled
3204    double * unscaledCopy;
3205    unscaledCopy = lowerCopy + numberTotal; 
3206    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3207    memcpy(rowLower_,unscaledCopy+numberColumns_,
3208           numberRows_*sizeof(double));
3209    unscaledCopy = upperCopy + numberTotal; 
3210    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3211    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3212           numberRows_*sizeof(double));
3213  }
3214  delete [] saveLower;
3215  delete [] saveUpper;
3216#ifdef ALL_DANTZIG
3217  if (savePivot) {
3218#endif
3219    delete dualRowPivot_;
3220    dualRowPivot_ = savePivot;
3221#ifdef ALL_DANTZIG
3222  }
3223#endif
3224  // Restore any saved stuff
3225  restoreData(data);
3226  perturbation_ = savePerturbation;
3227  delete rowArray_[4];
3228  rowArray_[4]=NULL;
3229  delete rowArray_[5];
3230  rowArray_[5]=NULL;
3231  char line[100];
3232  sprintf(line,"Ending theta %g\n", endingTheta);
3233  handler_->message(CLP_GENERAL,messages_)
3234    << line << CoinMessageEol;
3235  return problemStatus_;
3236}
3237int
3238ClpSimplexOther::parametricsLoop(parametricsData & paramData,
3239                                 ClpDataSave & data,bool canSkipFactorization)
3240{
3241  double & startingTheta = paramData.startingTheta;
3242  double & endingTheta = paramData.endingTheta;
3243  int numberTotal = numberRows_+numberColumns_;
3244  // stuff is already at starting
3245  const int * lowerList = paramData.lowerList;
3246  const int * upperList = paramData.upperList;
3247  problemStatus_ = -1;
3248  //double saveEndingTheta=endingTheta;
3249
3250  // This says whether to restore things etc
3251  // startup will have factorized so can skip
3252  int factorType = 0;
3253  // Start check for cycles
3254  progress_.startCheck();
3255  // Say change made on first iteration
3256  changeMade_ = 1;
3257  /*
3258    Status of problem:
3259    0 - optimal
3260    1 - infeasible
3261    2 - unbounded
3262    -1 - iterating
3263    -2 - factorization wanted
3264    -3 - redo checking without factorization
3265    -4 - looks infeasible
3266  */
3267  while (problemStatus_ < 0) {
3268    int iRow, iColumn;
3269    // clear
3270    for (iRow = 0; iRow < 6; iRow++) {
3271      rowArray_[iRow]->clear();
3272    }
3273   
3274    for (iColumn = 0; iColumn < 2; iColumn++) {
3275      columnArray_[iColumn]->clear();
3276    }
3277   
3278    // give matrix (and model costs and bounds a chance to be
3279    // refreshed (normally null)
3280    matrix_->refresh(this);
3281    // may factorize, checks if problem finished
3282    if (!canSkipFactorization)
3283      statusOfProblemInParametrics(factorType, data);
3284    canSkipFactorization=false;
3285    if (numberPrimalInfeasibilities_) {
3286      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3287        // treat as success
3288        problemStatus_=0;
3289        endingTheta=startingTheta;
3290        break;
3291      }
3292      // probably can get rid of this if we adjust every change in theta
3293      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3294      //     sumPrimalInfeasibilities_);
3295      const double * lowerChange = lower_+numberTotal;
3296      const double * upperChange = upper_+numberTotal;
3297      const double * startLower = lowerChange+numberTotal;
3298      const double * startUpper = upperChange+numberTotal;
3299      //startingTheta -= 1.0e-7;
3300      int nLowerChange = lowerList[-1];
3301      for (int i = 0; i < nLowerChange; i++) {
3302        int iSequence = lowerList[i];
3303        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3304      }
3305      int nUpperChange = upperList[-1];
3306      for (int i = 0; i < nUpperChange; i++) {
3307        int iSequence = upperList[i];
3308        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3309      }
3310      // adjust rhs in case dual uses
3311      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3312      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3313      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3314      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3315      if (rowScale_) {
3316        for (int i=0;i<numberColumns_;i++) {
3317          double multiplier = columnScale_[i];
3318          if (columnLower_[i]>-1.0e20)
3319            columnLower_[i] *= multiplier;
3320          if (columnUpper_[i]<1.0e20)
3321            columnUpper_[i] *= multiplier;
3322        }
3323        for (int i=0;i<numberRows_;i++) {
3324          double multiplier = inverseRowScale_[i];
3325          if (rowLower_[i]>-1.0e20)
3326            rowLower_[i] *= multiplier;
3327          if (rowUpper_[i]<1.0e20)
3328            rowUpper_[i] *= multiplier;
3329        }
3330      }
3331      double * saveDuals = NULL;
3332      problemStatus_=-1;
3333      ClpObjective * saveObjective = objective_;
3334      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3335      if (saveObjective!=objective_) {
3336        delete objective_;
3337        objective_=saveObjective;
3338      }
3339      int pass=100;
3340      double moved=0.0;
3341      while(sumPrimalInfeasibilities_) {
3342        //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_,
3343        //     sumPrimalInfeasibilities_);
3344        pass--;
3345        if (!pass)
3346          break;
3347        problemStatus_=-1;
3348        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3349          double value=solution_[iSequence];
3350          if (value<lower_[iSequence]-1.0e-9) {
3351            moved += lower_[iSequence]-value;
3352            lower_[iSequence]=value;
3353          } else if (value>upper_[iSequence]+1.0e-9) {
3354            moved += upper_[iSequence]-value;
3355            upper_[iSequence]=value;
3356          }
3357        }
3358        if (!moved) {
3359          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3360            double value=solution_[iSequence];
3361            if (value<lower_[iSequence]-1.0e-9) {
3362              moved += lower_[iSequence]-value;
3363              lower_[iSequence]=value;
3364            } else if (value>upper_[iSequence]+1.0e-9) {
3365              moved += upper_[iSequence]-value;
3366              upper_[iSequence]=value;
3367            }
3368          }
3369        }
3370        assert (moved);
3371        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3372      }
3373      // adjust
3374      //printf("Should adjust - moved %g\n",moved);
3375    }
3376    // Say good factorization
3377    factorType = 1;
3378    if (data.sparseThreshold_) {
3379      // use default at present
3380      factorization_->sparseThreshold(0);
3381      factorization_->goSparse();
3382    }
3383   
3384    // exit if victory declared
3385    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3386      break;
3387   
3388    // test for maximum iterations
3389    if (hitMaximumIterations()) {
3390      problemStatus_ = 3;
3391      break;
3392    }
3393#ifdef CLP_USER_DRIVEN
3394    // Check event
3395    {
3396      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3397      if (status >= 0) {
3398        problemStatus_ = 5;
3399        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3400        break;
3401      }
3402    }
3403#endif
3404    // Do iterations
3405    problemStatus_=-1;
3406    whileIterating(paramData, 0.0,
3407                   NULL);
3408    //startingTheta = endingTheta;
3409    //endingTheta = saveEndingTheta;
3410  }
3411  if (!problemStatus_/*||problemStatus_==2*/) {
3412    theta_ = endingTheta;
3413#ifdef CLP_USER_DRIVEN
3414    {
3415      double saveTheta=theta_;
3416      theta_ = endingTheta;
3417      int status=eventHandler_->event(ClpEventHandler::theta);
3418      if (status>=0&&status<10) {
3419        endingTheta=theta_;
3420        theta_=saveTheta;
3421        problemStatus_=-1;
3422      } else {
3423        if (status>=10) {
3424          problemStatus_=status-10;
3425          startingTheta=endingTheta;
3426        }
3427        theta_=saveTheta;
3428      }
3429    }
3430#endif
3431    return 0;
3432  } else if (problemStatus_ == 10) {
3433    return -1;
3434  } else {
3435    return problemStatus_;
3436  }
3437}
3438/* Checks if finished.  Updates status */
3439void
3440ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3441{
3442     if (type == 2) {
3443          // trouble - go to recovery
3444          problemStatus_ = 10;
3445          return;
3446     }
3447     if (problemStatus_ > -3 || factorization_->pivots()) {
3448          // factorize
3449          // later on we will need to recover from singularities
3450          // also we could skip if first time
3451          if (type) {
3452               // is factorization okay?
3453               if (internalFactorize(1)) {
3454                    // trouble - go to recovery
3455                    problemStatus_ = 10;
3456                    return;
3457               }
3458          }
3459          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3460               problemStatus_ = -3;
3461     }
3462     // at this stage status is -3 or -4 if looks infeasible
3463     // get primal and dual solutions
3464     gutsOfSolution(NULL, NULL);
3465     double realDualInfeasibilities = sumDualInfeasibilities_;
3466     // If bad accuracy treat as singular
3467     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3468          // trouble - go to recovery
3469          problemStatus_ = 10;
3470          return;
3471     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3472          // Can reduce tolerance
3473          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3474          factorization_->pivotTolerance(newTolerance);
3475     }
3476     // Check if looping
3477     int loop;
3478     if (type != 2)
3479          loop = progress_.looping();
3480     else
3481          loop = -1;
3482     if (loop >= 0) {
3483          problemStatus_ = loop; //exit if in loop
3484          if (!problemStatus_) {
3485               // declaring victory
3486               numberPrimalInfeasibilities_ = 0;
3487               sumPrimalInfeasibilities_ = 0.0;
3488          } else {
3489               problemStatus_ = 10; // instead - try other algorithm
3490          }
3491          return;
3492     } else if (loop < -1) {
3493          // something may have changed
3494          gutsOfSolution(NULL, NULL);
3495     }
3496     progressFlag_ = 0; //reset progress flag
3497     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3498          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3499                    << numberIterations_ << objectiveValue();
3500          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3501                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3502          handler_->printing(sumDualInfeasibilities_ > 0.0)
3503                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3504          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3505                             < numberDualInfeasibilities_)
3506                    << numberDualInfeasibilitiesWithoutFree_;
3507          handler_->message() << CoinMessageEol;
3508     }
3509#ifdef CLP_USER_DRIVEN
3510     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3511       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3512       if (status>=0) {
3513         // fix up
3514         for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3515           double value=solution_[iSequence];
3516           if (value<=lower_[iSequence]-primalTolerance_) {
3517             lower_[iSequence]=value;
3518           } else if (value>=upper_[iSequence]+primalTolerance_) {
3519             upper_[iSequence]=value;
3520           }
3521         }
3522         numberPrimalInfeasibilities_ = 0;
3523         sumPrimalInfeasibilities_ = 0.0;
3524       }
3525     }
3526#endif
3527     /* If we are primal feasible and any dual infeasibilities are on
3528        free variables then it is better to go to primal */
3529     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3530               numberDualInfeasibilities_) {
3531          problemStatus_ = 10;
3532          return;
3533     }
3534
3535     // check optimal
3536     // give code benefit of doubt
3537     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3538               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3539          // say optimal (with these bounds etc)
3540          numberDualInfeasibilities_ = 0;
3541          sumDualInfeasibilities_ = 0.0;
3542          numberPrimalInfeasibilities_ = 0;
3543          sumPrimalInfeasibilities_ = 0.0;
3544     }
3545     if (dualFeasible() || problemStatus_ == -4) {
3546          progress_.modifyObjective(objectiveValue_
3547                                    - sumDualInfeasibilities_ * dualBound_);
3548     }
3549     if (numberPrimalInfeasibilities_) {
3550          if (problemStatus_ == -4 || problemStatus_ == -5) {
3551               problemStatus_ = 1; // infeasible
3552          }
3553     } else if (numberDualInfeasibilities_) {
3554          // clean up
3555          problemStatus_ = 10;
3556     } else {
3557          problemStatus_ = 0;
3558     }
3559     lastGoodIteration_ = numberIterations_;
3560     if (problemStatus_ < 0) {
3561          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3562          if (sumDualInfeasibilities_)
3563               numberDualInfeasibilities_ = 1;
3564     }
3565     // Allow matrices to be sorted etc
3566     int fake = -999; // signal sort
3567     matrix_->correctSequence(this, fake, fake);
3568}
3569//static double lastThetaX=0.0;
3570/* This has the flow between re-factorizations
3571   Reasons to come out:
3572   -1 iterations etc
3573   -2 inaccuracy
3574   -3 slight inaccuracy (and done iterations)
3575   +0 looks optimal (might be unbounded - but we will investigate)
3576   +1 looks infeasible
3577   +3 max iterations
3578   +4 accuracy problems
3579*/
3580int
3581ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
3582                                const double * /*changeObjective*/)
3583{
3584  double & startingTheta = paramData.startingTheta;
3585  double & endingTheta = paramData.endingTheta;
3586  const double * lowerChange = paramData.lowerChange;
3587  const double * upperChange = paramData.upperChange;
3588  int numberTotal = numberColumns_ + numberRows_;
3589  const int * lowerList = paramData.lowerList;
3590  const int * upperList = paramData.upperList;
3591  //#define CLP_PARAMETRIC_DENSE_ARRAYS 2
3592#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3593  double * lowerGap = paramData.lowerGap;
3594  double * upperGap = paramData.upperGap;
3595  double * lowerCoefficient = paramData.lowerCoefficient;
3596  double * upperCoefficient = paramData.upperCoefficient;
3597#endif
3598  // do basic pointers
3599  int * backwardBasic = paramData.backwardBasic;
3600  for (int i=0;i<numberTotal;i++)
3601    backwardBasic[i]=-1;
3602  for (int i=0;i<numberRows_;i++) {
3603    int iPivot=pivotVariable_[i];
3604    backwardBasic[iPivot]=i;
3605  }
3606     {
3607          int i;
3608          for (i = 0; i < 4; i++) {
3609               rowArray_[i]->clear();
3610          }
3611          for (i = 0; i < 2; i++) {
3612               columnArray_[i]->clear();
3613          }
3614     }
3615     // if can't trust much and long way from optimal then relax
3616     if (largestPrimalError_ > 10.0)
3617          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3618     else
3619          factorization_->relaxAccuracyCheck(1.0);
3620     // status stays at -1 while iterating, >=0 finished, -2 to invert
3621     // status -3 to go to top without an invert
3622     int returnCode = -1;
3623     double lastTheta = startingTheta;
3624     double useTheta = startingTheta;
3625     while (problemStatus_ == -1) {
3626          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3627          // Get theta for bounds - we know can't crossover
3628          int pivotType = nextTheta(1, increaseTheta, paramData,
3629                                     NULL);
3630          useTheta += theta_;
3631          double change = useTheta - lastTheta;
3632          if (change>1.0e-14) {
3633            int n;
3634            n=lowerList[-1];
3635            for (int i=0;i<n;i++) {
3636              int iSequence = lowerList[i];
3637              double thisChange = change * lowerChange[iSequence];
3638              double newValue = lower_[iSequence] + thisChange;
3639              lower_[iSequence] = newValue;
3640#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3641              if (getStatus(iSequence)==basic) {
3642                int iRow=backwardBasic[iSequence];
3643                lowerGap[iRow] -= thisChange;
3644              } else if(getStatus(iSequence)==atLowerBound) {
3645                solution_[iSequence] = newValue;
3646              }
3647#else
3648              if(getStatus(iSequence)==atLowerBound) {
3649                solution_[iSequence] = newValue;
3650              }
3651#endif
3652#if 0
3653              // may have to adjust other bound
3654              double otherValue = upper_[iSequence];
3655              if (otherValue-newValue<dualBound_) {
3656                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3657                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3658                //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5);
3659              }
3660#endif
3661            }
3662            n=upperList[-1];
3663            for (int i=0;i<n;i++) {
3664              int iSequence = upperList[i];
3665              double thisChange = change * upperChange[iSequence];
3666              double newValue = upper_[iSequence] + thisChange;
3667              upper_[iSequence] = newValue;
3668              if(getStatus(iSequence)==atUpperBound||
3669                 getStatus(iSequence)==isFixed) {
3670                solution_[iSequence] = newValue;
3671#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3672              } else if (getStatus(iSequence)==basic) {
3673                int iRow=backwardBasic[iSequence];
3674                upperGap[iRow] += thisChange;
3675#endif
3676              }
3677              // may have to adjust other bound
3678              double otherValue = lower_[iSequence];
3679              if (newValue-otherValue<dualBound_) {
3680                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3681                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3682                //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5);
3683              }
3684            }
3685          }
3686          sequenceIn_=-1;
3687          if (pivotType) {
3688            if (useTheta>lastTheta+1.0e-9) {
3689              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3690                << useTheta << objectiveValue() << CoinMessageEol;
3691              lastTheta = useTheta;
3692            }
3693            problemStatus_ = -2;
3694            if (!factorization_->pivots()&&pivotRow_<0)
3695              problemStatus_=2;
3696#ifdef CLP_USER_DRIVEN
3697            {
3698              double saveTheta=theta_;
3699              theta_ = endingTheta;
3700              if (problemStatus_==2&&theta_>paramData.acceptableMaxTheta)
3701                theta_=COIN_DBL_MAX; // we have finished
3702              int status=eventHandler_->event(ClpEventHandler::theta);
3703              if (status>=0&&status<10) {
3704                endingTheta=theta_;
3705                problemStatus_=-1;
3706                continue;
3707              } else {
3708                if (status>=10)
3709                  problemStatus_=status-10;
3710                if (status<0)
3711                  startingTheta = useTheta;
3712              }
3713              theta_=saveTheta;
3714            }
3715#else
3716            startingTheta = useTheta;
3717#endif
3718            return 4;
3719          }
3720          // choose row to go out
3721          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3722          if (pivotRow_ >= 0) {
3723               // we found a pivot row
3724               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3725                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3726                              << pivotRow_
3727                              << CoinMessageEol;
3728               }
3729               // check accuracy of weights
3730               dualRowPivot_->checkAccuracy();
3731               // do ratio test for normal iteration
3732               double bestPossiblePivot = bestPivot();
3733               if (sequenceIn_ >= 0) {
3734                    // normal iteration
3735                    // update the incoming column
3736                    double btranAlpha = -alpha_ * directionOut_; // for check
3737#ifndef COIN_FAC_NEW
3738                    unpackPacked(rowArray_[1]);
3739#else
3740                    unpack(rowArray_[1]);
3741#endif
3742                    // and update dual weights (can do in parallel - with extra array)
3743                    rowArray_[2]->clear();
3744                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3745                                                          rowArray_[2],
3746                                                          rowArray_[3],
3747                                                          rowArray_[1]);
3748                    // see if update stable
3749#ifdef CLP_DEBUG
3750                    if ((handler_->logLevel() & 32))
3751                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3752#endif
3753                    double checkValue = 1.0e-7;
3754                    // if can't trust much and long way from optimal then relax
3755                    if (largestPrimalError_ > 10.0)
3756                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3757                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3758                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3759                         handler_->message(CLP_DUAL_CHECK, messages_)
3760                                   << btranAlpha
3761                                   << alpha_
3762                                   << CoinMessageEol;
3763                         // clear arrays
3764                         rowArray_[4]->clear();
3765                         if (factorization_->pivots()) {
3766                              dualRowPivot_->unrollWeights();
3767                              problemStatus_ = -2; // factorize now
3768                              rowArray_[0]->clear();
3769                              rowArray_[1]->clear();
3770                              columnArray_[0]->clear();
3771                              returnCode = -2;
3772                              break;
3773                         } else {
3774                              // take on more relaxed criterion
3775                              double test;
3776                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3777                                   test = 1.0e-1 * fabs(alpha_);
3778                              else
3779                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
3780                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3781                                        fabs(btranAlpha - alpha_) > test) {
3782                                   dualRowPivot_->unrollWeights();
3783                                   // need to reject something
3784                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
3785                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
3786                                             << x << sequenceWithin(sequenceOut_)
3787                                             << CoinMessageEol;
3788                                   setFlagged(sequenceOut_);
3789                                   progress_.clearBadTimes();
3790                                   lastBadIteration_ = numberIterations_; // say be more cautious
3791                                   rowArray_[0]->clear();
3792                                   rowArray_[1]->clear();
3793                                   columnArray_[0]->clear();
3794                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3795                                        //printf("I think should declare infeasible\n");
3796                                        problemStatus_ = 1;
3797                                        returnCode = 1;
3798                                        break;
3799                                   }
3800                                   continue;
3801                              }
3802                         }
3803                    }
3804                    // update duals BEFORE replaceColumn so can do updateColumn
3805                    double objectiveChange = 0.0;
3806                    // do duals first as variables may flip bounds
3807                    // rowArray_[0] and columnArray_[0] may have flips
3808                    // so use rowArray_[3] for work array from here on
3809                    int nswapped = 0;
3810                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3811                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3812#if CLP_CAN_HAVE_ZERO_OBJ
3813                    if ((specialOptions_&2097152)==0) {
3814#endif
3815                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3816                                                                                               rowArray_[2], theta_,
3817                                                                                               objectiveChange, false);
3818                      assert (!nswapped);
3819#if CLP_CAN_HAVE_ZERO_OBJ
3820                    } else {
3821                      rowArray_[0]->clear();
3822                      rowArray_[2]->clear();
3823                      columnArray_[0]->clear();
3824                    }
3825#endif
3826                    // which will change basic solution
3827                    if (nswapped) {
3828                      abort(); //needs testing
3829                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3830                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
3831                                                             1.0, objectiveChange);
3832                         // recompute dualOut_
3833                         valueOut_ = solution_[sequenceOut_];
3834                         if (directionOut_ < 0) {
3835                              dualOut_ = valueOut_ - upperOut_;
3836                         } else {
3837                              dualOut_ = lowerOut_ - valueOut_;
3838                         }
3839                    }
3840                    // amount primal will move
3841                    double movement = -dualOut_ * directionOut_ / alpha_;
3842                    // so objective should increase by fabs(dj)*movement
3843                    // but we already have objective change - so check will be good
3844                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3845#ifdef CLP_DEBUG
3846                         if (handler_->logLevel() & 32)
3847                              printf("movement %g, swap change %g, rest %g  * %g\n",
3848                                     objectiveChange + fabs(movement * dualIn_),
3849                                     objectiveChange, movement, dualIn_);
3850#endif
3851                         assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3852                         if(factorization_->pivots()) {
3853                              // going backwards - factorize
3854                              dualRowPivot_->unrollWeights();
3855                              problemStatus_ = -2; // factorize now
3856                              returnCode = -2;
3857                              break;
3858                         }
3859                    }
3860                    CoinAssert(fabs(dualOut_) < 1.0e50);
3861                    // if stable replace in basis
3862                    int updateStatus = factorization_->replaceColumn(this,
3863                                       rowArray_[2],
3864                                       rowArray_[1],
3865                                       pivotRow_,
3866                                       alpha_);
3867                    // if no pivots, bad update but reasonable alpha - take and invert
3868                    if (updateStatus == 2 &&
3869                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3870                         updateStatus = 4;
3871                    if (updateStatus == 1 || updateStatus == 4) {
3872                         // slight error
3873                         if (factorization_->pivots() > 5 || updateStatus == 4) {
3874                              problemStatus_ = -2; // factorize now
3875                              returnCode = -3;
3876                         }
3877                    } else if (updateStatus == 2) {
3878                         // major error
3879                         dualRowPivot_->unrollWeights();
3880                         // later we may need to unwind more e.g. fake bounds
3881                         if (factorization_->pivots()) {
3882                              problemStatus_ = -2; // factorize now
3883                              returnCode = -2;
3884                              break;
3885                         } else {
3886                              // need to reject something
3887                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3888                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3889                                        << x << sequenceWithin(sequenceOut_)
3890                                        << CoinMessageEol;
3891                              setFlagged(sequenceOut_);
3892                              progress_.clearBadTimes();
3893                              lastBadIteration_ = numberIterations_; // say be more cautious
3894                              rowArray_[0]->clear();
3895                              rowArray_[1]->clear();
3896                              columnArray_[0]->clear();
3897                              // make sure dual feasible
3898                              // look at all rows and columns
3899                              double objectiveChange = 0.0;
3900                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3901                                        0.0, objectiveChange, true);
3902                              continue;
3903                         }
3904                    } else if (updateStatus == 3) {
3905                         // out of memory
3906                         // increase space if not many iterations
3907                         if (factorization_->pivots() <
3908                                   0.5 * factorization_->maximumPivots() &&
3909                                   factorization_->pivots() < 200)
3910                              factorization_->areaFactor(
3911                                   factorization_->areaFactor() * 1.1);
3912                         problemStatus_ = -2; // factorize now
3913                    } else if (updateStatus == 5) {
3914                         problemStatus_ = -2; // factorize now
3915                    }
3916                    int * lowerActive = paramData.lowerActive;
3917                    int * upperActive = paramData.upperActive;
3918                    // update change vector
3919                    {
3920                      double * work = rowArray_[1]->denseVector();
3921                      int number = rowArray_[1]->getNumElements();
3922                      int * which = rowArray_[1]->getIndices();
3923                      assert (!rowArray_[4]->packedMode());
3924#ifndef COIN_FAC_NEW
3925                      assert (rowArray_[1]->packedMode());
3926#else
3927                      assert (!rowArray_[1]->packedMode());
3928#endif
3929                      double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3930                      double multiplier = -pivotValue/alpha_;
3931                      double * array=rowArray_[4]->denseVector();
3932#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3933                      int lowerN=lowerActive[-1];
3934                      int upperN=upperActive[-1];
3935#endif
3936                      if (multiplier) {
3937                        for (int i = 0; i < number; i++) {
3938                          int iRow = which[i];
3939#ifndef COIN_FAC_NEW
3940                          double alpha=multiplier*work[i];
3941#else
3942                          double alpha=multiplier*work[iRow];
3943#endif
3944#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3945                          double alpha3 = alpha+array[iRow];
3946                          int iSequence = pivotVariable_[iRow];
3947                          double oldLower = lowerCoefficient[iRow];
3948                          double oldUpper = upperCoefficient[iRow];
3949                          if (lower_[iSequence]>-1.0e30) {
3950                            //lowerGap[iRow]=value-lower_[iSequence];
3951                            double alpha2 = alpha3 + lowerChange[iSequence];
3952                            if (alpha2>1.0e-8)  {
3953                              lowerCoefficient[iRow]=alpha2;
3954                              if (!oldLower)
3955                                lowerActive[lowerN++]=iRow;
3956                            } else {
3957                              if (oldLower)
3958                                lowerCoefficient[iRow]=COIN_DBL_MIN;
3959                            }
3960                          } else {
3961                            if (oldLower)
3962                              lowerCoefficient[iRow]=COIN_DBL_MIN;
3963                          }
3964                          if (upper_[iSequence]<1.0e30) {
3965                            //upperGap[iRow]=-(value-upper_[iSequence]);
3966                            double alpha2 = -(alpha3+upperChange[iSequence]);
3967                            if (alpha2>1.0e-8) {
3968                              upperCoefficient[iRow]=alpha2;
3969                              if (!oldUpper)
3970                                upperActive[upperN++]=iRow;
3971                            } else {
3972                              if (oldUpper)
3973                                upperCoefficient[iRow]=COIN_DBL_MIN;
3974                            }
3975                          } else {
3976                            if (oldUpper)
3977                              upperCoefficient[iRow]=COIN_DBL_MIN;
3978                          }
3979#endif
3980                          rowArray_[4]->quickAdd(iRow,alpha);
3981                        }
3982                      }
3983                      pivotValue = array[pivotRow_];
3984                      // we want pivot to be -multiplier
3985                      rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3986#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3987                      assert (lowerN>=0&&lowerN<=numberRows_);
3988                      lowerActive[-1]=lowerN;
3989                      upperActive[-1]=upperN;
3990#endif
3991                    }
3992                    // update primal solution
3993#if CLP_CAN_HAVE_ZERO_OBJ
3994                    if ((specialOptions_&2097152)!=0) 
3995                      theta_=0.0;
3996#endif
3997                    if (theta_ < 0.0) {
3998#ifdef CLP_DEBUG
3999                         if (handler_->logLevel() & 32)
4000                              printf("negative theta %g\n", theta_);
4001#endif
4002                         theta_ = 0.0;
4003                    }
4004                    // do actual flips
4005                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
4006                    //rowArray_[1]->expand();
4007#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4008                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
4009                                                        movement,
4010                                                        objectiveChange);
4011#else
4012                    // do by hand
4013                    {
4014                      double * work = rowArray_[1]->denseVector();
4015                      int number = rowArray_[1]->getNumElements();
4016                      int * which = rowArray_[1]->getIndices();
4017                      int i;
4018                      if (rowArray_[1]->packedMode()) {
4019                        for (i = 0; i < number; i++) {
4020                          int iRow = which[i];
4021                          int iSequence = pivotVariable_[iRow];
4022                          double value = solution_[iSequence];
4023                          double change = movement * work[i];
4024                          value -= change;
4025                          if (lower_[iSequence]>-1.0e30)
4026                            lowerGap[iRow]=value-lower_[iSequence];
4027                          if (upper_[iSequence]<1.0e30)
4028                            upperGap[iRow]=-(value-upper_[iSequence]);
4029                          solution_[iSequence] = value;
4030                          objectiveChange -= change * cost_[iSequence];
4031                          work[i] = 0.0;
4032                        }
4033                      } else {
4034                        for (i = 0; i < number; i++) {
4035                          int iRow = which[i];
4036                          int iSequence = pivotVariable_[iRow];
4037                          double value = solution_[iSequence];
4038                          double change = movement * work[iRow];
4039                          value -= change;
4040                          solution_[iSequence] = value;
4041                          objectiveChange -= change * cost_[iSequence];
4042                          work[iRow] = 0.0;
4043                        }
4044                      }
4045                      rowArray_[1]->setNumElements(0);
4046                    }
4047#endif
4048                    // modify dualout
4049                    dualOut_ /= alpha_;
4050                    dualOut_ *= -directionOut_;
4051                    //setStatus(sequenceIn_,basic);
4052                    dj_[sequenceIn_] = 0.0;
4053                    //double oldValue = valueIn_;
4054                    if (directionIn_ == -1) {
4055                         // as if from upper bound
4056                         valueIn_ = upperIn_ + dualOut_;
4057                    } else {
4058                         // as if from lower bound
4059                         valueIn_ = lowerIn_ + dualOut_;
4060                    }
4061                    objectiveChange = 0.0;
4062#if CLP_CAN_HAVE_ZERO_OBJ
4063                    if ((specialOptions_&2097152)==0) {
4064#endif
4065                      for (int i=0;i<numberTotal;i++)
4066                        objectiveChange += solution_[i]*cost_[i];
4067                      objectiveChange -= objectiveValue_;
4068#if CLP_CAN_HAVE_ZERO_OBJ
4069                    }
4070#endif
4071                    // outgoing
4072                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
4073                    lowerOut_=lower_[sequenceOut_];
4074                    upperOut_=upper_[sequenceOut_];
4075                    // set dj to zero unless values pass
4076                    if (directionOut_ > 0) {
4077                         valueOut_ = lowerOut_;
4078                         dj_[sequenceOut_] = theta_;
4079#if CLP_CAN_HAVE_ZERO_OBJ>1
4080#ifdef COIN_REUSE_RANDOM
4081                         if ((specialOptions_&2097152)!=0) {
4082                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
4083                         }
4084#endif
4085#endif
4086                    } else {
4087                         valueOut_ = upperOut_;
4088                         dj_[sequenceOut_] = -theta_;
4089#if CLP_CAN_HAVE_ZERO_OBJ>1
4090#ifdef COIN_REUSE_RANDOM
4091                         if ((specialOptions_&2097152)!=0) {
4092                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
4093                         }
4094#endif
4095#endif
4096                    }
4097                    solution_[sequenceOut_] = valueOut_;
4098                    int whatNext = housekeeping(objectiveChange);
4099                    reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_);
4100                    assert (backwardBasic[sequenceOut_]==pivotRow_);
4101                    backwardBasic[sequenceOut_]=-1;
4102                    backwardBasic[sequenceIn_]=pivotRow_;
4103#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4104                    double value = solution_[sequenceIn_];
4105                    double alpha = rowArray_[4]->denseVector()[pivotRow_];
4106                    double oldLower = lowerCoefficient[pivotRow_];
4107                    double oldUpper = upperCoefficient[pivotRow_];
4108                    if (lower_[sequenceIn_]>-1.0e30) {
4109                      lowerGap[pivotRow_]=value-lower_[sequenceIn_];
4110                      double alpha2 = alpha + lowerChange[sequenceIn_];
4111                      if (alpha2>1.0e-8)  {
4112                        lowerCoefficient[pivotRow_]=alpha2;
4113                        if (!oldLower) {
4114                          int lowerN=lowerActive[-1];
4115                          assert (lowerN>=0&&lowerN<numberRows_);
4116                          lowerActive[lowerN]=pivotRow_;
4117                          lowerActive[-1]=lowerN+1;
4118                        }
4119                      } else {
4120                        if (oldLower)
4121                          lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4122                      }
4123                    } else {
4124                      if (oldLower)
4125                        lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4126                    }
4127                    if (upper_[sequenceIn_]<1.0e30) {
4128                      upperGap[pivotRow_]=-(value-upper_[sequenceIn_]);
4129                      double alpha2 = -(alpha+upperChange[sequenceIn_]);
4130                      if (alpha2>1.0e-8) {
4131                        upperCoefficient[pivotRow_]=alpha2;
4132                        if (!oldUpper) {
4133                          int upperN=upperActive[-1];
4134                          assert (upperN>=0&&upperN<numberRows_);
4135                          upperActive[upperN]=pivotRow_;
4136                          upperActive[-1]=upperN+1;
4137                        }
4138                      } else {
4139                        if (oldUpper)
4140                          upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4141                      }
4142                    } else {
4143                      if (oldUpper)
4144                        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4145                    }
4146#endif
4147                    {
4148                      char in[200],out[200];
4149                      int iSequence=sequenceIn_;
4150                      if (iSequence<numberColumns_) {
4151                        if (lengthNames_) 
4152                          strcpy(in,columnNames_[iSequence].c_str());
4153                         else 
4154                          sprintf(in,"C%7.7d",iSequence);
4155                      } else {
4156                        iSequence -= numberColumns_;
4157                        if (lengthNames_) 
4158                          strcpy(in,rowNames_[iSequence].c_str());
4159                         else 
4160                          sprintf(in,"R%7.7d",iSequence);
4161                      }
4162                      iSequence=sequenceOut_;
4163                      if (iSequence<numberColumns_) {
4164                        if (lengthNames_) 
4165                          strcpy(out,columnNames_[iSequence].c_str());
4166                         else 
4167                          sprintf(out,"C%7.7d",iSequence);
4168                      } else {
4169                        iSequence -= numberColumns_;
4170                        if (lengthNames_) 
4171                          strcpy(out,rowNames_[iSequence].c_str());
4172                         else 
4173                          sprintf(out,"R%7.7d",iSequence);
4174                      }
4175                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
4176                        << useTheta << objectiveValue() 
4177                        << in << out << CoinMessageEol;
4178                    }
4179                    if (useTheta>lastTheta+1.0e-9) {
4180                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
4181                        << useTheta << objectiveValue() << CoinMessageEol;
4182                      lastTheta = useTheta;
4183                    }
4184                    // and set bounds correctly
4185                    originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
4186                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
4187                    if (whatNext == 1) {
4188                         problemStatus_ = -2; // refactorize
4189                    } else if (whatNext == 2) {
4190                         // maximum iterations or equivalent
4191                         problemStatus_ = 3;
4192                         returnCode = 3;
4193                         break;
4194                    }
4195#ifdef CLP_USER_DRIVEN
4196                    // Check event
4197                    {
4198                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
4199                         if (status >= 0) {
4200                              problemStatus_ = 5;
4201                              secondaryStatus_ = ClpEventHandler::endOfIteration;
4202                              returnCode = 4;
4203                              break;
4204                         }
4205                    }
4206#endif
4207               } else {
4208                    // no incoming column is valid
4209#ifdef CLP_USER_DRIVEN
4210                 rowArray_[0]->clear();
4211                 columnArray_[0]->clear();
4212                 theta_ = useTheta;
4213                 lastTheta = useTheta;
4214                 int action = eventHandler_->event(ClpEventHandler::noTheta);
4215                 if (action>=0) {
4216                   endingTheta=theta_;
4217                   theta_ = 0.0;
4218                   //adjust [4] from handler - but
4219                   //rowArray_[4]->clear(); // temp
4220                   if (action>=0&&action<10)
4221                     problemStatus_=-1; // carry on
4222                   else if (action==15)
4223                     problemStatus_ =5; // say stopped
4224                   returnCode = 1;
4225                   if (action==0||action>=10) 
4226                     break;
4227                   else
4228                     continue;
4229                 } else {
4230                 theta_ = 0.0;
4231                 }
4232#endif
4233                    pivotRow_ = -1;
4234#ifdef CLP_DEBUG
4235                    if (handler_->logLevel() & 32)
4236                         printf("** no column pivot\n");
4237#endif
4238                    if (factorization_->pivots() < 10) { 
4239                         // If we have just factorized and infeasibility reasonable say infeas
4240                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
4241                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
4242                                        || (specialOptions_ & 64) == 0) {
4243                                   // say infeasible
4244                                   problemStatus_ = 1;
4245                                   // unless primal feasible!!!!
4246                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
4247                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
4248                                   if (numberDualInfeasibilities_)
4249                                        problemStatus_ = 10;
4250                                   rowArray_[0]->clear();
4251                                   columnArray_[0]->clear();
4252                              }
4253                         }
4254                         // If special option set - put off as long as possible
4255                         if ((specialOptions_ & 64) == 0) {
4256                              problemStatus_ = -4; //say looks infeasible
4257                         } else {
4258                              // flag
4259                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
4260                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
4261                                        << x << sequenceWithin(sequenceOut_)
4262                                        << CoinMessageEol;
4263                              setFlagged(sequenceOut_);
4264                              if (!factorization_->pivots()) {
4265                                   rowArray_[0]->clear();
4266                                   columnArray_[0]->clear();
4267                                   continue;
4268                              }
4269                         }
4270                    }
4271                    rowArray_[0]->clear();
4272                    columnArray_[0]->clear();
4273                    returnCode = 1;
4274                    break;
4275               }
4276          } else {
4277               // no pivot row
4278#ifdef CLP_USER_DRIVEN
4279            {
4280              double saveTheta=theta_;
4281              theta_ = endingTheta;
4282              int status=eventHandler_->event(ClpEventHandler::theta);
4283              if (status>=0&&status<10) {
4284                endingTheta=theta_;
4285                theta_=saveTheta;
4286                continue;
4287              } else {
4288                theta_=saveTheta;
4289              }
4290            }
4291#endif
4292#ifdef CLP_DEBUG
4293               if (handler_->logLevel() & 32)
4294                    printf("** no row pivot\n");
4295#endif
4296               int numberPivots = factorization_->pivots();
4297               bool specialCase;
4298               int useNumberFake;
4299               returnCode = 0;
4300               if (numberPivots < 20 &&
4301                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
4302                         && dualBound_ > 1.0e8) {
4303                    specialCase = true;
4304                    // as dual bound high - should be okay
4305                    useNumberFake = 0;
4306               } else {
4307                    specialCase = false;
4308                    useNumberFake = numberFake_;
4309               }
4310               if (!numberPivots || specialCase) {
4311                    // may have crept through - so may be optimal
4312                    // check any flagged variables
4313                    int iRow;
4314                    for (iRow = 0; iRow < numberRows_; iRow++) {
4315                         int iPivot = pivotVariable_[iRow];
4316                         if (flagged(iPivot))
4317                              break;
4318                    }
4319                    if (iRow < numberRows_ && numberPivots) {
4320                         // try factorization
4321                         returnCode = -2;
4322                    }
4323
4324                    if (useNumberFake || numberDualInfeasibilities_) {
4325                         // may be dual infeasible
4326                         problemStatus_ = -5;
4327                    } else {
4328                         if (iRow < numberRows_) {
4329                              problemStatus_ = -5;
4330                         } else {
4331                              if (numberPivots) {
4332                                   // objective may be wrong
4333                                   objectiveValue_ = innerProduct(cost_,
4334                                                                  numberColumns_ + numberRows_,
4335                                                                  solution_);
4336                                   objectiveValue_ += objective_->nonlinearOffset();
4337                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
4338                                   if ((specialOptions_ & 16384) == 0) {
4339                                        // and dual_ may be wrong (i.e. for fixed or basic)
4340                                        CoinIndexedVector * arrayVector = rowArray_[1];
4341                                        arrayVector->clear();
4342                                        int iRow;
4343                                        double * array = arrayVector->denseVector();
4344                                        /* Use dual_ instead of array
4345                                           Even though dual_ is only numberRows_ long this is
4346                                           okay as gets permuted to longer rowArray_[2]
4347                                        */
4348                                        arrayVector->setDenseVector(dual_);
4349                                        int * index = arrayVector->getIndices();
4350                                        int number = 0;
4351                                        for (iRow = 0; iRow < numberRows_; iRow++) {
4352                                             int iPivot = pivotVariable_[iRow];
4353                                             double value = cost_[iPivot];
4354                                             dual_[iRow] = value;
4355                                             if (value) {
4356                                                  index[number++] = iRow;
4357                                             }
4358                                        }
4359                                        arrayVector->setNumElements(number);
4360                                        // Extended duals before "updateTranspose"
4361                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
4362                                        // Btran basic costs
4363                                        rowArray_[2]->clear();
4364                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
4365                                        // and return vector
4366                                        arrayVector->setDenseVector(array);
4367                                   }
4368                              }
4369                              problemStatus_ = 0;
4370                              sumPrimalInfeasibilities_ = 0.0;
4371                              if ((specialOptions_&(1024 + 16384)) != 0) {
4372                                   CoinIndexedVector * arrayVector = rowArray_[1];
4373                                   arrayVector->clear();
4374                                   double * rhs = arrayVector->denseVector();
4375                                   times(1.0, solution_, rhs);
4376                                   bool bad2 = false;
4377                                   int i;
4378                                   for ( i = 0; i < numberRows_; i++) {
4379                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
4380                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
4381                                             bad2 = true;
4382                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
4383                                        }
4384                                        rhs[i] = 0.0;
4385                                   }
4386                                   for ( i = 0; i < numberColumns_; i++) {
4387                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
4388                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4389                                             bad2 = true;
4390                                        }
4391                                   }
4392                                   if (bad2) {
4393                                        problemStatus_ = -3;
4394                                        returnCode = -2;
4395                                        // Force to re-factorize early next time
4396                                        int numberPivots = factorization_->pivots();
4397                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4398                                   }
4399                              }
4400                         }
4401                    }
4402               } else {
4403                    problemStatus_ = -3;
4404                    returnCode = -2;
4405                    // Force to re-factorize early next time
4406                    int numberPivots = factorization_->pivots();
4407                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4408               }
4409               break;
4410          }
4411     }
4412     startingTheta = lastTheta+theta_;
4413     return returnCode;
4414}
4415#if 0
4416static int zzzzzz=0;
4417int zzzzzzOther=0;
4418#endif
4419// Finds best possible pivot
4420double 
4421ClpSimplexOther::bestPivot(bool justColumns)
4422{
4423  // Get good size for pivot
4424  // Allow first few iterations to take tiny
4425  double acceptablePivot = 1.0e-9;
4426  if (numberIterations_ > 100)
4427    acceptablePivot = 1.0e-8;
4428  if (factorization_->pivots() > 10 ||
4429      (factorization_->pivots() && sumDualInfeasibilities_))
4430    acceptablePivot = 1.0e-5; // if we have iterated be more strict
4431  else if (factorization_->pivots() > 5)
4432    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4433  else if (factorization_->pivots())
4434    acceptablePivot = 1.0e-8; // relax
4435  double bestPossiblePivot = 1.0;
4436  // get sign for finding row of tableau
4437  // normal iteration
4438  // create as packed
4439  double direction = directionOut_;
4440#ifndef COIN_FAC_NEW
4441  rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4442#else
4443  rowArray_[0]->createOneUnpackedElement(pivotRow_, direction);
4444#endif
4445  factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4446  // put row of tableau in rowArray[0] and columnArray[0]
4447  matrix_->transposeTimes(this, -1.0,
4448                          rowArray_[0], rowArray_[3], columnArray_[0]);
4449  sequenceIn_=-1;
4450  if (justColumns)
4451    rowArray_[0]->clear();
4452  // do ratio test for normal iteration
4453  bestPossiblePivot = 
4454    reinterpret_cast<ClpSimplexDual *> 
4455    ( this)->dualColumn(rowArray_[0],
4456                        columnArray_[0], columnArray_[1],
4457                        rowArray_[3], acceptablePivot, NULL);
4458  return bestPossiblePivot;
4459}
4460// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4461int
4462ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData,
4463                           const double * /*changeObjective*/)
4464{
4465  const double * lowerChange = paramData.lowerChange;
4466  const double * upperChange = paramData.upperChange;
4467  const int * lowerList = paramData.lowerList;
4468  const int * upperList = paramData.upperList;
4469  int iSequence;
4470  bool toLower = false;
4471  //assert (type==1);
4472  // may need to decide based on model?
4473  bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4474  double * array = rowArray_[4]->denseVector();
4475  //rowArray_[4]->checkClean();
4476  const int * row = matrix_->getIndices();
4477  const int * columnLength = matrix_->getVectorLengths();
4478  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4479  const double * elementByColumn = matrix_->getElements();
4480#if 0
4481  double tempArray[5000];
4482  bool checkIt=false;
4483  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
4484    memcpy(tempArray,array,numberRows_*sizeof(double));
4485    checkIt=true;
4486    needFullUpdate=true;
4487  }
4488#endif
4489#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4490  double * lowerGap = paramData.lowerGap;
4491  double * upperGap = paramData.upperGap;
4492  double * lowerCoefficient = paramData.lowerCoefficient;
4493  double * upperCoefficient = paramData.upperCoefficient;
4494  int * lowerActive=paramData.lowerActive;
4495  int * upperActive=paramData.upperActive;
4496#endif
4497  if (!factorization_->pivots()||needFullUpdate) {
4498    //zzzzzz=0;
4499    rowArray_[4]->clear();
4500    // get change
4501    if (!rowScale_) {
4502      int n;
4503      n=lowerList[-2];
4504      int i;
4505      for (i=0;i<n;i++) {
4506        int iSequence = lowerList[i];
4507        assert (iSequence<numberColumns_);
4508        if (getColumnStatus(iSequence)==atLowerBound) {
4509          double value=lowerChange[iSequence];
4510          for (CoinBigIndex j = columnStart[iSequence];
4511               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4512            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4513          }
4514        }
4515      }
4516      n=lowerList[-1];
4517      const double * change = lowerChange+numberColumns_;
4518      for (;i<n;i++) {
4519        int iSequence = lowerList[i]-numberColumns_;
4520        assert (iSequence>=0);
4521        if (getRowStatus(iSequence)==atLowerBound) {
4522          double value=change[iSequence];
4523          rowArray_[4]->quickAdd(iSequence, -value);
4524        }
4525      }
4526      n=upperList[-2];
4527      for (i=0;i<n;i++) {
4528        int iSequence = upperList[i];
4529        assert (iSequence<numberColumns_);
4530        if (getColumnStatus(iSequence)==atUpperBound) {
4531          double value=upperChange[iSequence];
4532          for (CoinBigIndex j = columnStart[iSequence];
4533               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4534            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4535          }
4536        }
4537      }
4538      n=upperList[-1];
4539      change = upperChange+numberColumns_;
4540      for (;i<n;i++) {
4541        int iSequence = upperList[i]-numberColumns_;
4542        assert (iSequence>=0);
4543        if (getRowStatus(iSequence)==atUpperBound) {
4544          double value=change[iSequence];
4545          rowArray_[4]->quickAdd(iSequence, -value);
4546        }
4547      }
4548    } else {
4549      int n;
4550      n=lowerList[-2];
4551      int i;
4552      for (i=0;i<n;i++) {
4553        int iSequence = lowerList[i];
4554        assert (iSequence<numberColumns_);
4555        if (getColumnStatus(iSequence)==atLowerBound) {
4556          double value=lowerChange[iSequence];
4557          // apply scaling
4558          double scale = columnScale_[iSequence];
4559          for (CoinBigIndex j = columnStart[iSequence];
4560               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4561            int iRow = row[j];
4562            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4563          }
4564        }
4565      }
4566      n=lowerList[-1];
4567      const double * change = lowerChange+numberColumns_;
4568      for (;i<n;i++) {
4569        int iSequence = lowerList[i]-numberColumns_;
4570        assert (iSequence>=0);
4571        if (getRowStatus(iSequence)==atLowerBound) {
4572          double value=change[iSequence];
4573          rowArray_[4]->quickAdd(iSequence, -value);
4574        }
4575      }
4576      n=upperList[-2];
4577      for (i=0;i<n;i++) {
4578        int iSequence = upperList[i];
4579        assert (iSequence<numberColumns_);
4580        if (getColumnStatus(iSequence)==atUpperBound) {
4581          double value=upperChange[iSequence];
4582          // apply scaling
4583          double scale = columnScale_[iSequence];
4584          for (CoinBigIndex j = columnStart[iSequence];
4585               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4586            int iRow = row[j];
4587            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4588          }
4589        }
4590      }
4591      n=upperList[-1];
4592      change = upperChange+numberColumns_;
4593      for (;i<n;i++) {
4594        int iSequence = upperList[i]-numberColumns_;
4595        assert (iSequence>=0);
4596        if (getRowStatus(iSequence)==atUpperBound) {
4597          double value=change[iSequence];
4598          rowArray_[4]->quickAdd(iSequence, -value);
4599        }
4600      }
4601    }
4602    // ftran it
4603    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4604#if 0
4605    if (checkIt) {
4606      for (int i=0;i<numberRows_;i++) {
4607        assert (fabs(tempArray[i]-array[i])<1.0e-8);
4608      }
4609    }
4610#endif
4611#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4612    /* later for sparse - keep like CoinIndexedvector
4613       and just redo here */
4614    int lowerN=0;
4615    int upperN=0;
4616    memset(lowerCoefficient,0,numberRows_*sizeof(double));
4617    memset(upperCoefficient,0,numberRows_*sizeof(double));
4618    for (int iRow=0;iRow<numberRows_;iRow++) {
4619      iSequence = pivotVariable_[iRow];
4620      double currentSolution = solution_[iSequence];
4621      double alpha = array[iRow];
4622      double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4623      double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4624      if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4625        double currentLower = lower_[iSequence];
4626        ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4627        double gap=currentSolution-currentLower;
4628        lowerGap[iRow]=gap;
4629        lowerCoefficient[iRow]=thetaCoefficientLower;
4630        lowerActive[lowerN++]=iRow;
4631        //} else {
4632        //lowerCoefficient[iRow]=0.0;
4633      }
4634      if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4635        double currentUpper = upper_[iSequence];
4636        ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4637        double gap2=-(currentSolution-currentUpper); //positive
4638        upperGap[iRow]=gap2;
4639        upperCoefficient[iRow]=-thetaCoefficientUpper;
4640        upperActive[upperN++]=iRow;
4641      }
4642    }
4643    assert (lowerN>=0&&lowerN<=numberRows_);
4644    lowerActive[-1]=lowerN;
4645    upperActive[-1]=upperN;
4646#endif
4647  } else if (sequenceIn_>=0) {
4648    //assert (sequenceIn_>=0);
4649    assert (sequenceOut_>=0);
4650    assert (sequenceIn_!=sequenceOut_);
4651    double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4652    int needed=0;
4653    assert (!rowArray_[5]->getNumElements());
4654    if (change) {
4655      if (sequenceIn_<numberColumns_) {
4656        if (!rowScale_) {
4657          for (CoinBigIndex i = columnStart[sequenceIn_];
4658               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4659            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4660          }
4661        } else {
4662          // apply scaling
4663          double scale = columnScale_[sequenceIn_];
4664          for (CoinBigIndex i = columnStart[sequenceIn_];
4665               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4666            int iRow = row[i];
4667            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4668          }
4669        }
4670      } else {
4671        rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4672      }
4673      needed++;
4674    }
4675    if (getStatus(sequenceOut_)==atLowerBound)
4676      change=lowerChange[sequenceOut_];
4677    else
4678      change=upperChange[sequenceOut_];
4679    if (change) {
4680      if (sequenceOut_<numberColumns_) {
4681        if (!rowScale_) {
4682          for (CoinBigIndex i = columnStart[sequenceOut_];
4683               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4684            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4685          }
4686        } else {
4687          // apply scaling
4688          double scale = columnScale_[sequenceOut_];
4689          for (CoinBigIndex i = columnStart[sequenceOut_];
4690               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4691            int iRow = row[i];
4692            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4693          }
4694        }
4695      } else {
4696        rowArray_[5]->quickAdd(sequenceOut_-numberColumns_,-change);
4697      }
4698      needed++;
4699    }
4700    //printf("seqin %d seqout %d needed %d\n",
4701    //     sequenceIn_,sequenceOut_,needed);
4702    if (needed) {
4703      // ftran it
4704      factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4705      // add
4706      double * array5 = rowArray_[5]->denseVector();
4707      int * index5 = rowArray_[5]->getIndices();
4708      int number5 = rowArray_[5]->getNumElements();
4709#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4710      int lowerN=lowerActive[-1];
4711      int upperN=upperActive[-1];
4712      int nIn4=rowArray_[4]->getNumElements();
4713      int * index4 = rowArray_[4]->getIndices();
4714#endif
4715      for (int i = 0; i < number5; i++) {
4716        int iPivot = index5[i];
4717#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4718        rowArray_[4]->quickAdd(iPivot,array5[iPivot]);
4719#else
4720        /* later for sparse - modify here */
4721        int iSequence = pivotVariable_[iPivot];
4722        double currentSolution = solution_[iSequence];
4723        double currentAlpha = array[iPivot];
4724        double alpha5 = array5[iPivot];
4725        double alpha = currentAlpha+alpha5;
4726        if (currentAlpha) {
4727          if (alpha) {
4728            array[iPivot] = alpha;
4729          } else {
4730            array[iPivot] = COIN_DBL_MIN;
4731          }
4732        } else {
4733          index4[nIn4++] = iPivot;
4734          array[iPivot] = alpha;
4735        }
4736        double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4737        double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4738        double oldLower = lowerCoefficient[iPivot];
4739        double oldUpper = upperCoefficient[iPivot];
4740        if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4741          double currentLower = lower_[iSequence];
4742          ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4743          double gap=currentSolution-currentLower;
4744          lowerGap[iPivot]=gap;
4745          lowerCoefficient[iPivot]=thetaCoefficientLower;
4746          if (!oldLower)
4747            lowerActive[lowerN++]=iPivot;
4748        } else {
4749          if (oldLower)
4750            lowerCoefficient[iPivot]=COIN_DBL_MIN;
4751        }
4752        if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4753          double currentUpper = upper_[iSequence];
4754          ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4755          double gap2=-(currentSolution-currentUpper); //positive
4756          upperGap[iPivot]=gap2;
4757          upperCoefficient[iPivot]=-thetaCoefficientUpper;
4758          if (!oldUpper)
4759            upperActive[upperN++]=iPivot;
4760        } else {
4761          if (oldUpper)
4762            upperCoefficient[iPivot]=COIN_DBL_MIN;
4763        }
4764#endif
4765        array5[iPivot]=0.0;
4766      }
4767      rowArray_[5]->setNumElements(0);
4768#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4769      rowArray_[4]->setNumElements(nIn4);
4770      assert (lowerN>=0&&lowerN<=numberRows_);
4771      lowerActive[-1]=lowerN;
4772      upperActive[-1]=upperN;
4773#endif
4774    }
4775  }
4776  const int * index = rowArray_[4]->getIndices();
4777  int number = rowArray_[4]->getNumElements();
4778#define TESTXX 0
4779#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4780  int * markDone = reinterpret_cast<int *>(paramData.markDone);
4781  int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT;
4782  memset(markDone,0,nToZero*sizeof(int));
4783  const int * backwardBasic = paramData.backwardBasic;
4784#endif
4785  // first ones with alpha
4786  double theta1=maxTheta;
4787  int pivotRow1=-1;
4788#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4789  int pivotRow2=-1;
4790  double theta2=maxTheta;
4791#endif
4792#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4793  for (int i=0;i<number;i++) {
4794    int iPivot=index[i];
4795    iSequence = pivotVariable_[iPivot];
4796    //assert(!markDone[iSequence]);
4797    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4798    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4799    markDone[word] |= ( 1 << bit );
4800    // solution value will be sol - theta*alpha
4801    // bounds will be bounds + change *theta
4802    double currentSolution = solution_[iSequence];
4803    double alpha = array[iPivot];
4804    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4805    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4806    if (thetaCoefficientLower > 1.0e-8) {
4807      double currentLower = lower_[iSequence];
4808      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4809      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4810      double gap=currentSolution-currentLower;
4811      if (thetaCoefficientLower*theta1>gap) {
4812        theta1 = gap/thetaCoefficientLower;
4813        //toLower=true;
4814        pivotRow1=iPivot;
4815      }
4816    }
4817    if (thetaCoefficientUpper < -1.0e-8) {
4818      double currentUpper = upper_[iSequence];
4819      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4820      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4821      double gap2=currentSolution-currentUpper; //negative
4822      if (thetaCoefficientUpper*theta2<gap2) {
4823        theta2 = gap2/thetaCoefficientUpper;
4824        //toLower=false;
4825        pivotRow2=iPivot;
4826      }
4827    }
4828  }
4829  // now others
4830  int nLook=lowerList[-1];
4831  for (int i=0;i<nLook;i++) {
4832    int iSequence = lowerList[i];
4833    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4834    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4835    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4836      double currentSolution = solution_[iSequence];
4837      double currentLower = lower_[iSequence];
4838      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4839      double thetaCoefficient = lowerChange[iSequence];
4840      if (thetaCoefficient > 0.0) {
4841        double gap=currentSolution-currentLower;
4842        if (thetaCoefficient*theta1>gap) {
4843          theta1 = gap/thetaCoefficient;
4844          //toLower=true;
4845          pivotRow1 = backwardBasic[iSequence];
4846        }
4847      }
4848    }
4849  }
4850  nLook=upperList[-1];
4851  for (int i=0;i<nLook;i++) {
4852    int iSequence = upperList[i];
4853    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4854    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4855    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4856      double currentSolution = solution_[iSequence];
4857      double currentUpper = upper_[iSequence];
4858      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4859      double thetaCoefficient = upperChange[iSequence];
4860      if (thetaCoefficient < 0) {
4861        double gap=currentSolution-currentUpper; //negative
4862        if (thetaCoefficient*theta2<gap) {
4863          theta2 = gap/thetaCoefficient;
4864          //toLower=false;
4865          pivotRow2 = backwardBasic[iSequence];
4866        }
4867      }
4868    }
4869  }
4870  if (theta2<theta1) {
4871    theta_=theta2;
4872    toLower=false;
4873    pivotRow_=pivotRow2;
4874  } else {
4875    theta_=theta1;
4876    toLower=true;
4877    pivotRow_=pivotRow1;
4878  }
4879#if 0 //TESTXX
4880#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4881  {
4882    double * checkArray = new double[numberRows_];
4883    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4884    int lowerN=lowerActive[-1];
4885    for (int i=0;i<lowerN;i++) {
4886      int iRow=lowerActive[i];
4887      int iSequence = pivotVariable_[iRow];
4888      double alpha = array[iRow];
4889      double thetaCoefficient = lowerChange[iSequence] + alpha;
4890      if (thetaCoefficient > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4891        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4892        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4893          abort();
4894        }
4895      } else {
4896        assert (fabs(checkArray[iRow])<1.0e-12);
4897        if (fabs(checkArray[iRow])>1.0e-12) {
4898          abort();
4899        }
4900      }
4901      checkArray[iRow]=0.0;
4902    }
4903    for (int i=0;i<numberRows_;i++) {
4904      assert (!checkArray[i]);
4905      if (checkArray[i])
4906        abort();
4907    }
4908    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
4909    int upperN=upperActive[-1];
4910    for (int i=0;i<upperN;i++) {
4911      int iRow=upperActive[i];
4912      int iSequence = pivotVariable_[iRow];
4913      double alpha = array[iRow];
4914      double thetaCoefficient = -(upperChange[iSequence] + alpha);
4915      if (thetaCoefficient > 1.0e-8&&upper_[iSequence]<1.0e30) {
4916        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4917        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4918          abort();
4919        }
4920      } else {
4921        assert (fabs(checkArray[iRow])<1.0e-12);
4922        if (fabs(checkArray[iRow])>1.0e-12) {
4923          abort();
4924        }
4925      }
4926      checkArray[iRow]=0.0;
4927    }
4928    for (int i=0;i<numberRows_;i++) {
4929      assert (!checkArray[i]);
4930      if (checkArray[i])
4931        abort();
4932    }
4933    delete [] checkArray;
4934  }
4935  double theta3=maxTheta;
4936  int pivotRow3=-1;
4937  int lowerN=lowerActive[-1];
4938  for (int i=0;i<lowerN;i++) {
4939    int iRow=lowerActive[i];
4940    double lowerC = lowerCoefficient[iRow];
4941    double gap=lowerGap[iRow];
4942    if (toLower&&iRow==pivotRow_) {
4943      assert (lowerC*theta3>gap-1.0e-8);
4944      if (lowerC*theta3<gap-1.0e-8)
4945        abort();
4946    }
4947    if (lowerC*theta3>gap&&lowerC!=COIN_DBL_MIN) {
4948      theta3 = gap/lowerC;
4949      pivotRow3=iRow;
4950    }
4951  }
4952  int pivotRow4=pivotRow3;
4953  double theta4=theta3;
4954  int upperN=upperActive[-1];
4955  for (int i=0;i<upperN;i++) {
4956    int iRow=upperActive[i];
4957    double upperC = upperCoefficient[iRow];
4958    double gap=upperGap[iRow];
4959    if (!toLower&&iRow==pivotRow_) {
4960      assert (upperC*theta3>gap-1.0e-8);
4961      if (upperC*theta3<gap-1.0e-8)
4962        abort();
4963    }
4964    if (upperC*theta4>gap&&upperC!=COIN_DBL_MIN) {
4965      theta4 = gap/upperC;
4966      pivotRow4=iRow;
4967    }
4968  }
4969  bool toLower3;
4970  if (theta4<theta3) {
4971    theta3=theta4;
4972    toLower3=false;
4973    pivotRow3=pivotRow4;
4974  } else {
4975    toLower3=true;
4976  }
4977  if (fabs(theta3-theta_)>1.0e-8)
4978    abort();
4979  if (toLower!=toLower3||pivotRow_!=pivotRow3) {
4980    printf("bad piv - good %d %g %s, bad %d %g %s\n",pivotRow_,theta_,toLower ? "toLower" : "toUpper",
4981           pivotRow3,theta3,toLower3 ? "toLower" : "toUpper");
4982    //zzzzzz++;
4983    if (true/*zzzzzz>zzzzzzOther*/) {
4984      printf("Swapping\n");
4985      pivotRow_=pivotRow3;
4986      theta_=theta3;
4987      toLower=toLower3;
4988    }
4989  }
4990#endif
4991#endif
4992#else
4993#if 0 //CLP_PARAMETRIC_DENSE_ARRAYS==2
4994  {
4995    double * checkArray = new double[numberRows_];
4996    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4997    int lowerN=lowerActive[-1];
4998    for (int i=0;i<lowerN;i++) {
4999      int iRow=lowerActive[i];
5000      checkArray[iRow]=0.0;
5001    }
5002    for (int i=0;i<numberRows_;i++) {
5003      assert (!checkArray[i]);
5004      if (checkArray[i])
5005        abort();
5006    }
5007    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
5008    int upperN=upperActive[-1];
5009    for (int i=0;i<upperN;i++) {
5010      int iRow=upperActive[i];
5011      checkArray[iRow]=0.0;
5012    }
5013    for (int i=0;i<numberRows_;i++) {
5014      assert (!checkArray[i]);
5015      if (checkArray[i])
5016        abort();
5017    }
5018    delete [] checkArray;
5019  }
5020#endif
5021  int lowerN=lowerActive[-1];
5022  for (int i=0;i<lowerN;i++) {
5023    int iRow=lowerActive[i];
5024    double lowerC = lowerCoefficient[iRow];
5025    double gap=lowerGap[iRow];
5026    if (lowerC*theta1>gap&&lowerC!=COIN_DBL_MIN) {
5027      theta1 = gap/lowerC;
5028      pivotRow1=iRow;
5029    }
5030  }
5031  pivotRow_=pivotRow1;
5032  theta_=theta1;
5033  int upperN=upperActive[-1];
5034  for (int i=0;i<upperN;i++) {
5035    int iRow=upperActive[i];
5036    double upperC = upperCoefficient[iRow];
5037    double gap=upperGap[iRow];
5038    if (upperC*theta1>gap&&upperC!=COIN_DBL_MIN) {
5039      theta1 = gap/upperC;
5040      pivotRow1=iRow;
5041    }
5042  }
5043  if (theta1<theta_) {
5044    theta_=theta1;
5045    toLower=false;
5046    pivotRow_=pivotRow1;
5047  } else {
5048    toLower=true;
5049  }
5050#endif
5051  theta_ = CoinMax(theta_,0.0);
5052  if (theta_>1.0e-15) {
5053    // update solution
5054    for (int iRow = 0; iRow < number; iRow++) {
5055      int iPivot = index[iRow];
5056      iSequence = pivotVariable_[iPivot];
5057      // solution value will be sol - theta*alpha
5058      double alpha = array[iPivot];
5059      double currentSolution = solution_[iSequence] - theta_ * alpha;
5060      solution_[iSequence] =currentSolution;
5061#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5062      if (lower_[iSequence]>-1.0e30)
5063        lowerGap[iPivot]=currentSolution-lower_[iSequence];
5064      if (upper_[iSequence]<1.0e30)
5065        upperGap[iPivot]=-(currentSolution-upper_[iSequence]);
5066#endif
5067    }
5068  }
5069#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5070  if (pivotRow_>=0&&false) {
5071    double oldValue = upperCoefficient[pivotRow_];
5072    double value = array[pivotRow_];
5073    if (value) {
5074      if (!oldValue) {
5075        int upperN=upperActive[-1];
5076        assert (upperN>=0&&upperN<numberRows_);
5077        upperActive[upperN]=pivotRow_;
5078        upperActive[-1]=upperN+1;
5079      }
5080    } else {
5081      if (oldValue)
5082        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
5083    }
5084  }
5085#endif
5086#if 0
5087  for (int i=0;i<numberTotal;i++)
5088    assert(!markDone[i]);
5089#endif
5090  if (pivotRow_ >= 0) {
5091    sequenceOut_ = pivotVariable_[pivotRow_];
5092    valueOut_ = solution_[sequenceOut_];
5093    lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
5094    upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
5095    if (!toLower) {
5096      directionOut_ = -1;
5097      dualOut_ = valueOut_ - upperOut_;
5098    } else {
5099      directionOut_ = 1;
5100      dualOut_ = lowerOut_ - valueOut_;
5101    }
5102    return 0;
5103  } else {
5104    //theta_=0.0;
5105    return -1;
5106  }
5107}
5108// Restores bound to original bound
5109void 
5110ClpSimplexOther::originalBound(int iSequence, double theta, 
5111                               const double * lowerChange,
5112                               const double * upperChange)
5113{
5114     if (getFakeBound(iSequence) != noFake) {
5115          numberFake_--;
5116          setFakeBound(iSequence, noFake);
5117          if (iSequence >= numberColumns_) {
5118               // rows
5119               int iRow = iSequence - numberColumns_;
5120               rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
5121               rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
5122               if (rowScale_) {
5123                    if (rowLowerWork_[iRow] > -1.0e50)
5124                         rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5125                    if (rowUpperWork_[iRow] < 1.0e50)
5126                         rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5127               } else if (rhsScale_ != 1.0) {
5128                    if (rowLowerWork_[iRow] > -1.0e50)
5129                         rowLowerWork_[iRow] *= rhsScale_;
5130                    if (rowUpperWork_[iRow] < 1.0e50)
5131                         rowUpperWork_[iRow] *= rhsScale_;
5132               }
5133          } else {
5134               // columns
5135               columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
5136               columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
5137               if (rowScale_) {
5138                    double multiplier = 1.0 * inverseColumnScale_[iSequence];
5139                    if (columnLowerWork_[iSequence] > -1.0e50)
5140                         columnLowerWork_[iSequence] *= multiplier * rhsScale_;
5141                    if (columnUpperWork_[iSequence] < 1.0e50)
5142                         columnUpperWork_[iSequence] *= multiplier * rhsScale_;
5143               } else if (rhsScale_ != 1.0) {
5144                    if (columnLowerWork_[iSequence] > -1.0e50)
5145                         columnLowerWork_[iSequence] *= rhsScale_;
5146                    if (columnUpperWork_[iSequence] < 1.0e50)
5147                         columnUpperWork_[iSequence] *= rhsScale_;
5148               }
5149          }
5150     }
5151}
5152/* Expands out all possible combinations for a knapsack
5153   If buildObj NULL then just computes space needed - returns number elements
5154   On entry numberOutput is maximum allowed, on exit it is number needed or
5155   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
5156   least space to return values which reconstruct input.
5157   Rows returned will be original rows but no entries will be returned for
5158   any rows all of whose entries are in knapsack.  So up to user to allow for this.
5159   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
5160   in expanded knapsack.  Values in buildRow and buildElement;
5161*/
5162int
5163ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
5164                                double * buildObj, CoinBigIndex * buildStart,
5165                                int * buildRow, double * buildElement, int reConstruct) const
5166{
5167     int iRow;
5168     int iColumn;
5169     // Get column copy
5170     CoinPackedMatrix * columnCopy = matrix();
5171     // Get a row copy in standard format
5172     CoinPackedMatrix matrixByRow;
5173     matrixByRow.reverseOrderedCopyOf(*columnCopy);
5174     const double * elementByRow = matrixByRow.getElements();
5175     const int * column = matrixByRow.getIndices();
5176     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5177     const int * rowLength = matrixByRow.getVectorLengths();
5178     CoinBigIndex j;
5179     int * whichColumn = new int [numberColumns_];
5180     int * whichRow = new int [numberRows_];
5181     int numJ = 0;
5182     // Get what other columns can compensate for
5183     double * lo = new double [numberRows_];
5184     double * high = new double [numberRows_];
5185     {
5186          // Use to get tight column bounds
5187          ClpSimplex tempModel(*this);
5188          tempModel.tightenPrimalBounds(0.0, 0, true);
5189          // Now another model without knapsacks
5190          int nCol = 0;
5191          for (iRow = 0; iRow < numberRows_; iRow++) {
5192               whichRow[iRow] = iRow;
5193          }
5194          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
5195               whichColumn[iColumn] = -1;
5196          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
5197               int iColumn = column[j];
5198               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
5199                    whichColumn[iColumn] = 0;
5200               } else {
5201                    assert (!columnLower_[iColumn]); // fix later
5202               }
5203          }
5204          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5205               if (whichColumn[iColumn] < 0)
5206                    whichColumn[nCol++] = iColumn;
5207          }
5208          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
5209          // Row copy
5210          CoinPackedMatrix matrixByRow;
5211          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
5212          const double * elementByRow = matrixByRow.getElements();
5213          const int * column = matrixByRow.getIndices();
5214          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5215          const int * rowLength = matrixByRow.getVectorLengths();
5216          const double * columnLower = tempModel2.getColLower();
5217          const double * columnUpper = tempModel2.getColUpper();
5218          for (iRow = 0; iRow < numberRows_; iRow++) {
5219               lo[iRow] = -COIN_DBL_MAX;
5220               high[iRow] = COIN_DBL_MAX;
5221               if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
5222
5223                    // possible row
5224                    int infiniteUpper = 0;
5225                    int infiniteLower = 0;
5226                    double maximumUp = 0.0;
5227                    double maximumDown = 0.0;
5228                    CoinBigIndex rStart = rowStart[iRow];
5229                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];