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

Last change on this file since 1929 was 1929, checked in by stefan, 7 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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