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

Last change on this file since 2030 was 2028, checked in by forrest, 5 years ago

contributed fix to parametrics

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