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

Last change on this file since 2329 was 2329, checked in by forrest, 9 months ago

new not malloc

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