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

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

minor changes to implement Aboca

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