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

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

minor update to improve accuarcy in parametrics

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