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

Last change on this file since 2226 was 2226, checked in by forrest, 3 years ago

slight mistake

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