source: stable/1.16/Clp/src/ClpSimplexOther.cpp @ 2382

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

should not have had a printf

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