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

Last change on this file since 2380 was 2380, checked in by forrest, 4 months ago

copy handler in dualize

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