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

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

out abort

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