Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/Idiot.cpp

    r2333 r2385  
    4242#define IDIOT_FIX_TOLERANCE 1e-6
    4343#define SMALL_IDIOT_FIX_TOLERANCE 1e-10
    44 int
    45 Idiot::dropping(IdiotResult result,
    46                 double tolerance,
    47                 double small,
    48                 int *nbad)
     44int Idiot::dropping(IdiotResult result,
     45  double tolerance,
     46  double small,
     47  int *nbad)
    4948{
    50      if (result.infeas <= small) {
    51           double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0;
    52           if (result.dropThis > tolerance * value) {
    53                *nbad = 0;
    54                return 1;
    55           } else {
    56                (*nbad)++;
    57                if (*nbad > 4) {
    58                     return 0;
    59                } else {
    60                     return 1;
    61                }
    62           }
    63      } else {
    64           *nbad = 0;
    65           return 1;
    66      }
     49  if (result.infeas <= small) {
     50    double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0;
     51    if (result.dropThis > tolerance * value) {
     52      *nbad = 0;
     53      return 1;
     54    } else {
     55      (*nbad)++;
     56      if (*nbad > 4) {
     57        return 0;
     58      } else {
     59        return 1;
     60      }
     61    }
     62  } else {
     63    *nbad = 0;
     64    return 1;
     65  }
    6766}
    6867// Deals with whenUsed and slacks
    69 int
    70 Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
    71                       double * colsol, const double * lower, const double * upper,
    72                       const double * rowLower, const double * rowUpper,
    73                       const double * cost, const double * element, double fixTolerance,
    74                       double & objValue, double & infValue, double & maxInfeasibility)
     68int Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
     69  double *colsol, const double *lower, const double *upper,
     70  const double *rowLower, const double *rowUpper,
     71  const double *cost, const double *element, double fixTolerance,
     72  double &objValue, double &infValue, double &maxInfeasibility)
    7573{
    76      int n = 0;
    77      if ((strategy_ & 16384) == 0) {
    78           for (int i = ordinaryStart; i < ordinaryEnd; i++) {
    79                if (colsol[i] > lower[i] + fixTolerance) {
    80                     if (colsol[i] < upper[i] - fixTolerance) {
    81                          n++;
    82                     } else {
    83                          colsol[i] = upper[i];
    84                     }
    85                     whenUsed_[i] = iteration;
    86                } else {
    87                     colsol[i] = lower[i];
    88                }
    89           }
    90           return n;
    91      } else {
     74  int n = 0;
     75  if ((strategy_ & 16384) == 0) {
     76    for (int i = ordinaryStart; i < ordinaryEnd; i++) {
     77      if (colsol[i] > lower[i] + fixTolerance) {
     78        if (colsol[i] < upper[i] - fixTolerance) {
     79          n++;
     80        } else {
     81          colsol[i] = upper[i];
     82        }
     83        whenUsed_[i] = iteration;
     84      } else {
     85        colsol[i] = lower[i];
     86      }
     87    }
     88    return n;
     89  } else {
    9290#ifdef COIN_DEVELOP
    93           printf("entering inf %g, obj %g\n", infValue, objValue);
    94 #endif
    95           int nrows = model_->getNumRows();
    96           int ncols = model_->getNumCols();
    97           int * posSlack = whenUsed_ + ncols;
    98           int * negSlack = posSlack + nrows;
    99           int * nextSlack = negSlack + nrows;
    100           double * rowsol = reinterpret_cast<double *> (nextSlack + ncols);
    101           memset(rowsol, 0, nrows * sizeof(double));
     91    printf("entering inf %g, obj %g\n", infValue, objValue);
     92#endif
     93    int nrows = model_->getNumRows();
     94    int ncols = model_->getNumCols();
     95    int *posSlack = whenUsed_ + ncols;
     96    int *negSlack = posSlack + nrows;
     97    int *nextSlack = negSlack + nrows;
     98    double *rowsol = reinterpret_cast< double * >(nextSlack + ncols);
     99    memset(rowsol, 0, nrows * sizeof(double));
    102100#ifdef OSI_IDIOT
    103           const CoinPackedMatrix * matrix = model_->getMatrixByCol();
     101    const CoinPackedMatrix *matrix = model_->getMatrixByCol();
    104102#else
    105           // safer for odd matrices
    106           const CoinPackedMatrix * matrix = model_->matrix();
    107           //ClpMatrixBase * matrix = model_->clpMatrix();
    108 #endif
    109           const int * row = matrix->getIndices();
    110           const CoinBigIndex * columnStart = matrix->getVectorStarts();
    111           const int * columnLength = matrix->getVectorLengths();
    112           //const double * element = matrix->getElements();
    113           int i;
    114           objValue = 0.0;
    115           infValue = 0.0;
    116           maxInfeasibility = 0.0;
    117           for ( i = 0; i < ncols; i++) {
    118                if (nextSlack[i] == -1) {
    119                     // not a slack
    120                     if (colsol[i] > lower[i] + fixTolerance) {
    121                          if (colsol[i] < upper[i] - fixTolerance) {
    122                               n++;
    123                               whenUsed_[i] = iteration;
    124                          } else {
    125                               colsol[i] = upper[i];
    126                          }
    127                          whenUsed_[i] = iteration;
    128                     } else {
    129                          colsol[i] = lower[i];
    130                     }
    131                     double value = colsol[i];
    132                     if (value) {
    133                          objValue += cost[i] * value;
    134                          CoinBigIndex j;
    135                          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    136                               int iRow = row[j];
    137                               rowsol[iRow] += value * element[j];
    138                          }
    139                     }
    140                }
    141           }
    142           // temp fix for infinite lbs - just limit to -1000
    143           for (i = 0; i < nrows; i++) {
    144                double rowSave = rowsol[i];
    145                int iCol;
    146                iCol = posSlack[i];
    147                if (iCol >= 0) {
    148                     // slide all slack down
    149                     double rowValue = rowsol[i];
    150                     CoinBigIndex j = columnStart[iCol];
    151                     double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
    152                     rowSave += (colsol[iCol] - lowerValue) * element[j];
    153                     colsol[iCol] = lowerValue;
    154                     while (nextSlack[iCol] >= 0) {
    155                          iCol = nextSlack[iCol];
    156                          double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
    157                          j = columnStart[iCol];
    158                          rowSave += (colsol[iCol] - lowerValue) * element[j];
    159                          colsol[iCol] = lowerValue;
    160                     }
    161                     iCol = posSlack[i];
    162                     while (rowValue < rowLower[i] && iCol >= 0) {
    163                          // want to increase
    164                          double distance = rowLower[i] - rowValue;
    165                          double value = element[columnStart[iCol]];
    166                          double thisCost = cost[iCol];
    167                          if (distance <= value*(upper[iCol] - colsol[iCol])) {
    168                               // can get there
    169                               double movement = distance / value;
    170                               objValue += movement * thisCost;
    171                               rowValue = rowLower[i];
    172                               colsol[iCol] += movement;
    173                          } else {
    174                               // can't get there
    175                               double movement = upper[iCol] - colsol[iCol];
    176                               objValue += movement * thisCost;
    177                               rowValue += movement * value;
    178                               colsol[iCol] = upper[iCol];
    179                               iCol = nextSlack[iCol];
    180                          }
    181                     }
    182                     if (iCol >= 0) {
    183                          // may want to carry on - because of cost?
    184                          while (iCol >= 0 && cost[iCol] < 0 && rowValue < rowUpper[i]) {
    185                               // want to increase
    186                               double distance = rowUpper[i] - rowValue;
    187                               double value = element[columnStart[iCol]];
    188                               double thisCost = cost[iCol];
    189                               if (distance <= value*(upper[iCol] - colsol[iCol])) {
    190                                    // can get there
    191                                    double movement = distance / value;
    192                                    objValue += movement * thisCost;
    193                                    rowValue = rowUpper[i];
    194                                    colsol[iCol] += movement;
    195                                    iCol = -1;
    196                               } else {
    197                                    // can't get there
    198                                    double movement = upper[iCol] - colsol[iCol];
    199                                    objValue += movement * thisCost;
    200                                    rowValue += movement * value;
    201                                    colsol[iCol] = upper[iCol];
    202                                    iCol = nextSlack[iCol];
    203                               }
    204                          }
    205                          if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance &&
    206                                    colsol[iCol] < upper[iCol] - fixTolerance) {
    207                               whenUsed_[i] = iteration;
    208                               n++;
    209                          }
    210                     }
    211                     rowsol[i] = rowValue;
    212                }
    213                iCol = negSlack[i];
    214                if (iCol >= 0) {
    215                     // slide all slack down
    216                     double rowValue = rowsol[i];
    217                     CoinBigIndex j = columnStart[iCol];
    218                     double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
    219                     rowSave += (colsol[iCol] - lowerValue) * element[j];
    220                     colsol[iCol] = lowerValue;
    221                     while (nextSlack[iCol] >= 0) {
    222                          iCol = nextSlack[iCol];
    223                          j = columnStart[iCol];
    224                          double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
    225                          rowSave += (colsol[iCol] - lowerValue) * element[j];
    226                          colsol[iCol] = lowerValue;
    227                     }
    228                     iCol = negSlack[i];
    229                     while (rowValue > rowUpper[i] && iCol >= 0) {
    230                          // want to increase
    231                          double distance = -(rowUpper[i] - rowValue);
    232                          double value = -element[columnStart[iCol]];
    233                          double thisCost = cost[iCol];
    234                          if (distance <= value*(upper[iCol] - lower[iCol])) {
    235                               // can get there
    236                               double movement = distance / value;
    237                               objValue += movement * thisCost;
    238                               rowValue = rowUpper[i];
    239                               colsol[iCol] += movement;
    240                          } else {
    241                               // can't get there
    242                               double movement = upper[iCol] - lower[iCol];
    243                               objValue += movement * thisCost;
    244                               rowValue -= movement * value;
    245                               colsol[iCol] = upper[iCol];
    246                               iCol = nextSlack[iCol];
    247                          }
    248                     }
    249                     if (iCol >= 0) {
    250                          // may want to carry on - because of cost?
    251                          while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) {
    252                               // want to increase
    253                               double distance = -(rowLower[i] - rowValue);
    254                               double value = -element[columnStart[iCol]];
    255                               double thisCost = cost[iCol];
    256                               if (distance <= value*(upper[iCol] - colsol[iCol])) {
    257                                    // can get there
    258                                    double movement = distance / value;
    259                                    objValue += movement * thisCost;
    260                                    rowValue = rowLower[i];
    261                                    colsol[iCol] += movement;
    262                                    iCol = -1;
    263                               } else {
    264                                    // can't get there
    265                                    double movement = upper[iCol] - colsol[iCol];
    266                                    objValue += movement * thisCost;
    267                                    rowValue -= movement * value;
    268                                    colsol[iCol] = upper[iCol];
    269                                    iCol = nextSlack[iCol];
    270                               }
    271                          }
    272                          if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance &&
    273                                    colsol[iCol] < upper[iCol] - fixTolerance) {
    274                               whenUsed_[i] = iteration;
    275                               n++;
    276                          }
    277                     }
    278                     rowsol[i] = rowValue;
    279                }
    280                double infeasibility = CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
    281                infValue += infeasibility;
    282                maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
    283                // just change
    284                rowsol[i] -= rowSave;
    285           }
    286           return n;
    287      }
     103    // safer for odd matrices
     104    const CoinPackedMatrix *matrix = model_->matrix();
     105    //ClpMatrixBase * matrix = model_->clpMatrix();
     106#endif
     107    const int *row = matrix->getIndices();
     108    const CoinBigIndex *columnStart = matrix->getVectorStarts();
     109    const int *columnLength = matrix->getVectorLengths();
     110    //const double * element = matrix->getElements();
     111    int i;
     112    objValue = 0.0;
     113    infValue = 0.0;
     114    maxInfeasibility = 0.0;
     115    for (i = 0; i < ncols; i++) {
     116      if (nextSlack[i] == -1) {
     117        // not a slack
     118        if (colsol[i] > lower[i] + fixTolerance) {
     119          if (colsol[i] < upper[i] - fixTolerance) {
     120            n++;
     121            whenUsed_[i] = iteration;
     122          } else {
     123            colsol[i] = upper[i];
     124          }
     125          whenUsed_[i] = iteration;
     126        } else {
     127          colsol[i] = lower[i];
     128        }
     129        double value = colsol[i];
     130        if (value) {
     131          objValue += cost[i] * value;
     132          CoinBigIndex j;
     133          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     134            int iRow = row[j];
     135            rowsol[iRow] += value * element[j];
     136          }
     137        }
     138      }
     139    }
     140    // temp fix for infinite lbs - just limit to -1000
     141    for (i = 0; i < nrows; i++) {
     142      double rowSave = rowsol[i];
     143      int iCol;
     144      iCol = posSlack[i];
     145      if (iCol >= 0) {
     146        // slide all slack down
     147        double rowValue = rowsol[i];
     148        CoinBigIndex j = columnStart[iCol];
     149        double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     150        rowSave += (colsol[iCol] - lowerValue) * element[j];
     151        colsol[iCol] = lowerValue;
     152        while (nextSlack[iCol] >= 0) {
     153          iCol = nextSlack[iCol];
     154          double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     155          j = columnStart[iCol];
     156          rowSave += (colsol[iCol] - lowerValue) * element[j];
     157          colsol[iCol] = lowerValue;
     158        }
     159        iCol = posSlack[i];
     160        while (rowValue < rowLower[i] && iCol >= 0) {
     161          // want to increase
     162          double distance = rowLower[i] - rowValue;
     163          double value = element[columnStart[iCol]];
     164          double thisCost = cost[iCol];
     165          if (distance <= value * (upper[iCol] - colsol[iCol])) {
     166            // can get there
     167            double movement = distance / value;
     168            objValue += movement * thisCost;
     169            rowValue = rowLower[i];
     170            colsol[iCol] += movement;
     171          } else {
     172            // can't get there
     173            double movement = upper[iCol] - colsol[iCol];
     174            objValue += movement * thisCost;
     175            rowValue += movement * value;
     176            colsol[iCol] = upper[iCol];
     177            iCol = nextSlack[iCol];
     178          }
     179        }
     180        if (iCol >= 0) {
     181          // may want to carry on - because of cost?
     182          while (iCol >= 0 && cost[iCol] < 0 && rowValue < rowUpper[i]) {
     183            // want to increase
     184            double distance = rowUpper[i] - rowValue;
     185            double value = element[columnStart[iCol]];
     186            double thisCost = cost[iCol];
     187            if (distance <= value * (upper[iCol] - colsol[iCol])) {
     188              // can get there
     189              double movement = distance / value;
     190              objValue += movement * thisCost;
     191              rowValue = rowUpper[i];
     192              colsol[iCol] += movement;
     193              iCol = -1;
     194            } else {
     195              // can't get there
     196              double movement = upper[iCol] - colsol[iCol];
     197              objValue += movement * thisCost;
     198              rowValue += movement * value;
     199              colsol[iCol] = upper[iCol];
     200              iCol = nextSlack[iCol];
     201            }
     202          }
     203          if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && colsol[iCol] < upper[iCol] - fixTolerance) {
     204            whenUsed_[i] = iteration;
     205            n++;
     206          }
     207        }
     208        rowsol[i] = rowValue;
     209      }
     210      iCol = negSlack[i];
     211      if (iCol >= 0) {
     212        // slide all slack down
     213        double rowValue = rowsol[i];
     214        CoinBigIndex j = columnStart[iCol];
     215        double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     216        rowSave += (colsol[iCol] - lowerValue) * element[j];
     217        colsol[iCol] = lowerValue;
     218        while (nextSlack[iCol] >= 0) {
     219          iCol = nextSlack[iCol];
     220          j = columnStart[iCol];
     221          double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     222          rowSave += (colsol[iCol] - lowerValue) * element[j];
     223          colsol[iCol] = lowerValue;
     224        }
     225        iCol = negSlack[i];
     226        while (rowValue > rowUpper[i] && iCol >= 0) {
     227          // want to increase
     228          double distance = -(rowUpper[i] - rowValue);
     229          double value = -element[columnStart[iCol]];
     230          double thisCost = cost[iCol];
     231          if (distance <= value * (upper[iCol] - lower[iCol])) {
     232            // can get there
     233            double movement = distance / value;
     234            objValue += movement * thisCost;
     235            rowValue = rowUpper[i];
     236            colsol[iCol] += movement;
     237          } else {
     238            // can't get there
     239            double movement = upper[iCol] - lower[iCol];
     240            objValue += movement * thisCost;
     241            rowValue -= movement * value;
     242            colsol[iCol] = upper[iCol];
     243            iCol = nextSlack[iCol];
     244          }
     245        }
     246        if (iCol >= 0) {
     247          // may want to carry on - because of cost?
     248          while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) {
     249            // want to increase
     250            double distance = -(rowLower[i] - rowValue);
     251            double value = -element[columnStart[iCol]];
     252            double thisCost = cost[iCol];
     253            if (distance <= value * (upper[iCol] - colsol[iCol])) {
     254              // can get there
     255              double movement = distance / value;
     256              objValue += movement * thisCost;
     257              rowValue = rowLower[i];
     258              colsol[iCol] += movement;
     259              iCol = -1;
     260            } else {
     261              // can't get there
     262              double movement = upper[iCol] - colsol[iCol];
     263              objValue += movement * thisCost;
     264              rowValue -= movement * value;
     265              colsol[iCol] = upper[iCol];
     266              iCol = nextSlack[iCol];
     267            }
     268          }
     269          if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance && colsol[iCol] < upper[iCol] - fixTolerance) {
     270            whenUsed_[i] = iteration;
     271            n++;
     272          }
     273        }
     274        rowsol[i] = rowValue;
     275      }
     276      double infeasibility = CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
     277      infValue += infeasibility;
     278      maxInfeasibility = CoinMax(maxInfeasibility, infeasibility);
     279      // just change
     280      rowsol[i] -= rowSave;
     281    }
     282    return n;
     283  }
    288284}
    289285
    290286/* returns -1 if none or start of costed slacks or -2 if
    291287   there are costed slacks but it is messy */
    292 static int countCostedSlacks(OsiSolverInterface * model)
     288static int countCostedSlacks(OsiSolverInterface *model)
    293289{
    294290#ifdef OSI_IDIOT
    295      const CoinPackedMatrix * matrix = model->getMatrixByCol();
     291  const CoinPackedMatrix *matrix = model->getMatrixByCol();
    296292#else
    297      // safer for odd matrices (note really ClpSimplex not OsiSolverInterface)
    298      const CoinPackedMatrix * matrix = model->matrix();
    299      //ClpMatrixBase * matrix = model_->clpMatrix();
    300 #endif
    301      const int * row = matrix->getIndices();
    302      const CoinBigIndex * columnStart = matrix->getVectorStarts();
    303      const int * columnLength = matrix->getVectorLengths();
    304      const double * element = matrix->getElements();
    305      const double * rowupper = model->getRowUpper();
    306      int nrows = model->getNumRows();
    307      int ncols = model->getNumCols();
    308      int slackStart = ncols - nrows;
    309      int nSlacks = nrows;
    310      int i;
     293  // safer for odd matrices (note really ClpSimplex not OsiSolverInterface)
     294  const CoinPackedMatrix *matrix = model->matrix();
     295  //ClpMatrixBase * matrix = model_->clpMatrix();
     296#endif
     297  const int *row = matrix->getIndices();
     298  const CoinBigIndex *columnStart = matrix->getVectorStarts();
     299  const int *columnLength = matrix->getVectorLengths();
     300  const double *element = matrix->getElements();
     301  const double *rowupper = model->getRowUpper();
     302  int nrows = model->getNumRows();
     303  int ncols = model->getNumCols();
     304  int slackStart = ncols - nrows;
     305  int nSlacks = nrows;
     306  int i;
    311307
    312      if (ncols <= nrows) return -1;
    313      while (1) {
    314           for (i = 0; i < nrows; i++) {
    315                int j = i + slackStart;
    316                CoinBigIndex k = columnStart[j];
    317                if (columnLength[j] == 1) {
    318                     if (row[k] != i || element[k] != 1.0) {
    319                          nSlacks = 0;
    320                          break;
    321                     }
    322                } else {
    323                     nSlacks = 0;
    324                     break;
    325                }
    326                if (rowupper[i] <= 0.0) {
    327                     nSlacks = 0;
    328                     break;
    329                }
    330           }
    331           if (nSlacks || !slackStart) break;
    332           slackStart = 0;
    333      }
    334      if (!nSlacks) slackStart = -1;
    335      return slackStart;
     308  if (ncols <= nrows)
     309    return -1;
     310  while (1) {
     311    for (i = 0; i < nrows; i++) {
     312      int j = i + slackStart;
     313      CoinBigIndex k = columnStart[j];
     314      if (columnLength[j] == 1) {
     315        if (row[k] != i || element[k] != 1.0) {
     316          nSlacks = 0;
     317          break;
     318        }
     319      } else {
     320        nSlacks = 0;
     321        break;
     322      }
     323      if (rowupper[i] <= 0.0) {
     324        nSlacks = 0;
     325        break;
     326      }
     327    }
     328    if (nSlacks || !slackStart)
     329      break;
     330    slackStart = 0;
     331  }
     332  if (!nSlacks)
     333    slackStart = -1;
     334  return slackStart;
    336335}
    337 void
    338 Idiot::crash(int numberPass, CoinMessageHandler * handler,
    339              const CoinMessages *messages, bool doCrossover)
     336void Idiot::crash(int numberPass, CoinMessageHandler *handler,
     337  const CoinMessages *messages, bool doCrossover)
    340338{
    341      // lightweight options
    342      int numberColumns = model_->getNumCols();
    343      const double * objective = model_->getObjCoefficients();
    344      int nnzero = 0;
    345      double sum = 0.0;
    346      int i;
    347      for (i = 0; i < numberColumns; i++) {
    348           if (objective[i]) {
    349                sum += fabs(objective[i]);
    350                nnzero++;
    351           }
    352      }
    353      sum /= static_cast<double> (nnzero + 1);
    354      if (maxIts_ == 5)
    355           maxIts_ = 2;
    356      if (numberPass <= 0)
    357           majorIterations_ = static_cast<int>(2 + log10(static_cast<double>(numberColumns + 1)));
    358      else
    359           majorIterations_ = numberPass;
    360      // If mu not changed then compute
    361      if (mu_ == 1e-4)
    362           mu_ = CoinMax(1.0e-3, sum * 1.0e-5);
    363      if (maxIts2_ == 100) {
    364           if (!lightWeight_) {
    365                maxIts2_ = 105;
    366           } else if (lightWeight_ == 1) {
    367                mu_ *= 1000.0;
    368                maxIts2_ = 23;
    369           } else if (lightWeight_ == 2) {
    370                maxIts2_ = 11;
    371           } else {
    372                maxIts2_ = 23;
    373           }
    374      }
    375      //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
    376      if (numberColumns)
    377        solve2(handler, messages);
     339  // lightweight options
     340  int numberColumns = model_->getNumCols();
     341  const double *objective = model_->getObjCoefficients();
     342  int nnzero = 0;
     343  double sum = 0.0;
     344  int i;
     345  for (i = 0; i < numberColumns; i++) {
     346    if (objective[i]) {
     347      sum += fabs(objective[i]);
     348      nnzero++;
     349    }
     350  }
     351  sum /= static_cast< double >(nnzero + 1);
     352  if (maxIts_ == 5)
     353    maxIts_ = 2;
     354  if (numberPass <= 0)
     355    majorIterations_ = static_cast< int >(2 + log10(static_cast< double >(numberColumns + 1)));
     356  else
     357    majorIterations_ = numberPass;
     358  // If mu not changed then compute
     359  if (mu_ == 1e-4)
     360    mu_ = CoinMax(1.0e-3, sum * 1.0e-5);
     361  if (maxIts2_ == 100) {
     362    if (!lightWeight_) {
     363      maxIts2_ = 105;
     364    } else if (lightWeight_ == 1) {
     365      mu_ *= 1000.0;
     366      maxIts2_ = 23;
     367    } else if (lightWeight_ == 2) {
     368      maxIts2_ = 11;
     369    } else {
     370      maxIts2_ = 23;
     371    }
     372  }
     373  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
     374  if (numberColumns)
     375    solve2(handler, messages);
    378376#ifndef OSI_IDIOT
    379      if (doCrossover) {
    380           double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
    381           if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0)
    382                crossOver(16 + 1);
    383           else
    384                crossOver(majorIterations_ < 1000000 ? 3 : 2);
    385      }
     377  if (doCrossover) {
     378    double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast< double >(model_->numberRows());
     379    if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0)
     380      crossOver(16 + 1);
     381    else
     382      crossOver(majorIterations_ < 1000000 ? 3 : 2);
     383  }
    386384#endif
    387385}
    388386
    389 void
    390 Idiot::solve()
     387void Idiot::solve()
    391388{
    392      CoinMessages dummy;
    393      solve2(NULL, &dummy);
     389  CoinMessages dummy;
     390  solve2(NULL, &dummy);
    394391}
    395 void
    396 Idiot::solve2(CoinMessageHandler * handler, const CoinMessages * messages)
     392void Idiot::solve2(CoinMessageHandler *handler, const CoinMessages *messages)
    397393{
    398      int strategy = 0;
    399      double d2;
    400      int i, n;
    401      int allOnes = 1;
    402      int iteration = 0;
    403      int iterationTotal = 0;
    404      int nTry = 0; /* number of tries at same weight */
    405      double fixTolerance = IDIOT_FIX_TOLERANCE;
    406      int maxBigIts = maxBigIts_;
    407      int maxIts = maxIts_;
    408      int logLevel = logLevel_;
    409      int saveMajorIterations = majorIterations_;
    410      majorIterations_ = majorIterations_ % 1000000;
    411      if (handler) {
    412           if (handler->logLevel() > 0 && handler->logLevel() < 3)
    413                logLevel = 1;
    414           else if (!handler->logLevel())
    415                logLevel = 0;
    416           else
    417                logLevel = 7;
    418      }
    419      double djExit = djTolerance_;
    420      double djFlag = 1.0 + 100.0 * djExit;
    421      double djTol = 0.00001;
    422      double mu = mu_;
    423      double drop = drop_;
    424      int maxIts2 = maxIts2_;
    425      double factor = muFactor_;
    426      double smallInfeas = smallInfeas_;
    427      double reasonableInfeas = reasonableInfeas_;
    428      double stopMu = stopMu_;
    429      double maxmin, offset;
    430      double lastWeighted = 1.0e50;
    431      double exitDrop = exitDrop_;
    432      double fakeSmall = smallInfeas;
    433      double firstInfeas;
    434      int badIts = 0;
    435      int slackStart, ordStart, ordEnd;
    436      int checkIteration = 0;
    437      int lambdaIteration = 0;
    438      int belowReasonable = 0; /* set if ever gone below reasonable infeas */
    439      double bestWeighted = 1.0e60;
    440      double bestFeasible = 1.0e60; /* best solution while feasible */
    441      IdiotResult result, lastResult;
    442      int saveStrategy = strategy_;
    443      const int strategies[] = {0, 2, 128};
    444      double saveLambdaScale = 0.0;
    445      if ((saveStrategy & 128) != 0) {
    446           fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
    447      }
     394  int strategy = 0;
     395  double d2;
     396  int i, n;
     397  int allOnes = 1;
     398  int iteration = 0;
     399  int iterationTotal = 0;
     400  int nTry = 0; /* number of tries at same weight */
     401  double fixTolerance = IDIOT_FIX_TOLERANCE;
     402  int maxBigIts = maxBigIts_;
     403  int maxIts = maxIts_;
     404  int logLevel = logLevel_;
     405  int saveMajorIterations = majorIterations_;
     406  majorIterations_ = majorIterations_ % 1000000;
     407  if (handler) {
     408    if (handler->logLevel() > 0 && handler->logLevel() < 3)
     409      logLevel = 1;
     410    else if (!handler->logLevel())
     411      logLevel = 0;
     412    else
     413      logLevel = 7;
     414  }
     415  double djExit = djTolerance_;
     416  double djFlag = 1.0 + 100.0 * djExit;
     417  double djTol = 0.00001;
     418  double mu = mu_;
     419  double drop = drop_;
     420  int maxIts2 = maxIts2_;
     421  double factor = muFactor_;
     422  double smallInfeas = smallInfeas_;
     423  double reasonableInfeas = reasonableInfeas_;
     424  double stopMu = stopMu_;
     425  double maxmin, offset;
     426  double lastWeighted = 1.0e50;
     427  double exitDrop = exitDrop_;
     428  double fakeSmall = smallInfeas;
     429  double firstInfeas;
     430  int badIts = 0;
     431  int slackStart, ordStart, ordEnd;
     432  int checkIteration = 0;
     433  int lambdaIteration = 0;
     434  int belowReasonable = 0; /* set if ever gone below reasonable infeas */
     435  double bestWeighted = 1.0e60;
     436  double bestFeasible = 1.0e60; /* best solution while feasible */
     437  IdiotResult result, lastResult;
     438  int saveStrategy = strategy_;
     439  const int strategies[] = { 0, 2, 128 };
     440  double saveLambdaScale = 0.0;
     441  if ((saveStrategy & 128) != 0) {
     442    fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
     443  }
    448444#ifdef OSI_IDIOT
    449      const CoinPackedMatrix * matrix = model_->getMatrixByCol();
     445  const CoinPackedMatrix *matrix = model_->getMatrixByCol();
    450446#else
    451      // safer for odd matrices
    452      const CoinPackedMatrix * matrix = model_->matrix();
    453      //ClpMatrixBase * matrix = model_->clpMatrix();
    454 #endif
    455      const int * row = matrix->getIndices();
    456      const CoinBigIndex * columnStart = matrix->getVectorStarts();
    457      const int * columnLength = matrix->getVectorLengths();
    458      const double * element = matrix->getElements();
    459      int nrows = model_->getNumRows();
    460      int ncols = model_->getNumCols();
    461      double * rowsol, * colsol;
    462      double * pi, * dj;
     447  // safer for odd matrices
     448  const CoinPackedMatrix *matrix = model_->matrix();
     449  //ClpMatrixBase * matrix = model_->clpMatrix();
     450#endif
     451  const int *row = matrix->getIndices();
     452  const CoinBigIndex *columnStart = matrix->getVectorStarts();
     453  const int *columnLength = matrix->getVectorLengths();
     454  const double *element = matrix->getElements();
     455  int nrows = model_->getNumRows();
     456  int ncols = model_->getNumCols();
     457  double *rowsol, *colsol;
     458  double *pi, *dj;
    463459#ifndef OSI_IDIOT
    464      double * cost = model_->objective();
    465      double * lower = model_->columnLower();
    466      double * upper = model_->columnUpper();
     460  double *cost = model_->objective();
     461  double *lower = model_->columnLower();
     462  double *upper = model_->columnUpper();
    467463#else
    468      double * cost = new double [ncols];
    469      CoinMemcpyN( model_->getObjCoefficients(), ncols, cost);
    470      const double * lower = model_->getColLower();
    471      const double * upper = model_->getColUpper();
    472 #endif
    473      const double *elemXX;
    474      double * saveSol;
    475      double * rowupper = new double[nrows]; // not const as modified
    476      CoinMemcpyN(model_->getRowUpper(), nrows, rowupper);
    477      double * rowlower = new double[nrows]; // not const as modified
    478      CoinMemcpyN(model_->getRowLower(), nrows, rowlower);
    479      CoinThreadRandom * randomNumberGenerator = model_->randomNumberGenerator();
    480      int * whenUsed;
    481      double * lambda;
    482      saveSol = new double[ncols];
    483      lambda = new double [nrows];
    484      rowsol = new double[nrows];
    485      colsol = new double [ncols];
    486      CoinMemcpyN(model_->getColSolution(), ncols, colsol);
    487      pi = new double[nrows];
    488      dj = new double[ncols];
     464  double *cost = new double[ncols];
     465  CoinMemcpyN(model_->getObjCoefficients(), ncols, cost);
     466  const double *lower = model_->getColLower();
     467  const double *upper = model_->getColUpper();
     468#endif
     469  const double *elemXX;
     470  double *saveSol;
     471  double *rowupper = new double[nrows]; // not const as modified
     472  CoinMemcpyN(model_->getRowUpper(), nrows, rowupper);
     473  double *rowlower = new double[nrows]; // not const as modified
     474  CoinMemcpyN(model_->getRowLower(), nrows, rowlower);
     475  CoinThreadRandom *randomNumberGenerator = model_->randomNumberGenerator();
     476  int *whenUsed;
     477  double *lambda;
     478  saveSol = new double[ncols];
     479  lambda = new double[nrows];
     480  rowsol = new double[nrows];
     481  colsol = new double[ncols];
     482  CoinMemcpyN(model_->getColSolution(), ncols, colsol);
     483  pi = new double[nrows];
     484  dj = new double[ncols];
    489485#ifndef OSI_IDIOT
    490      bool fixAfterSome = false; //(model_->specialOptions()&8388608)!=0;
    491      int exitAfter = 50; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
    492      {
    493         int numberColumns = model_->numberColumns();
    494         for (int i=0;i<numberColumns;i++) {
    495           if (upper[i]==lower[i])
    496             model_->setColumnStatus(i,ClpSimplex::isFixed);
    497         }
    498      }
    499 #endif
    500      delete [] whenUsed_;
    501      bool oddSlacks = false;
    502      // See if any costed slacks
    503      int numberSlacks = 0;
    504      for (i = 0; i < ncols; i++) {
    505           if (columnLength[i] == 1)
    506                numberSlacks++;
    507      }
    508      if (!numberSlacks || (strategy_&524288)!=0) {
    509           whenUsed_ = new int[ncols];
    510      } else {
     486  bool fixAfterSome = false; //(model_->specialOptions()&8388608)!=0;
     487  int exitAfter = 50; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
     488  {
     489    int numberColumns = model_->numberColumns();
     490    for (int i = 0; i < numberColumns; i++) {
     491      if (upper[i] == lower[i])
     492        model_->setColumnStatus(i, ClpSimplex::isFixed);
     493    }
     494  }
     495#endif
     496  delete[] whenUsed_;
     497  bool oddSlacks = false;
     498  // See if any costed slacks
     499  int numberSlacks = 0;
     500  for (i = 0; i < ncols; i++) {
     501    if (columnLength[i] == 1)
     502      numberSlacks++;
     503  }
     504  if (!numberSlacks || (strategy_ & 524288) != 0) {
     505    whenUsed_ = new int[ncols];
     506  } else {
    511507#ifdef COIN_DEVELOP
    512           printf("%d slacks\n", numberSlacks);
    513 #endif
    514           oddSlacks = true;
    515           int extra = static_cast<int> (nrows * sizeof(double) / sizeof(int));
    516           whenUsed_ = new int[2*ncols+2*nrows+extra];
    517           int * posSlack = whenUsed_ + ncols;
    518           int * negSlack = posSlack + nrows;
    519           int * nextSlack = negSlack + nrows;
    520           for (i = 0; i < nrows; i++) {
    521                posSlack[i] = -1;
    522                negSlack[i] = -1;
    523           }
    524           for (i = 0; i < ncols; i++)
    525                nextSlack[i] = -1;
    526           for (i = 0; i < ncols; i++) {
    527                if (columnLength[i] == 1) {
    528                     CoinBigIndex j = columnStart[i];
    529                     int iRow = row[j];
    530                     if (element[j] > 0.0) {
    531                          if (posSlack[iRow] == -1) {
    532                               posSlack[iRow] = i;
    533                          } else {
    534                               int iCol = posSlack[iRow];
    535                               while (nextSlack[iCol] >= 0)
    536                                    iCol = nextSlack[iCol];
    537                               nextSlack[iCol] = i;
    538                          }
    539                     } else {
    540                          if (negSlack[iRow] == -1) {
    541                               negSlack[iRow] = i;
    542                          } else {
    543                               int iCol = negSlack[iRow];
    544                               while (nextSlack[iCol] >= 0)
    545                                    iCol = nextSlack[iCol];
    546                               nextSlack[iCol] = i;
    547                          }
    548                     }
    549                }
    550           }
    551           // now sort
    552           for (i = 0; i < nrows; i++) {
    553                int iCol;
    554                iCol = posSlack[i];
    555                if (iCol >= 0) {
    556                     CoinBigIndex j = columnStart[iCol];
     508    printf("%d slacks\n", numberSlacks);
     509#endif
     510    oddSlacks = true;
     511    int extra = static_cast< int >(nrows * sizeof(double) / sizeof(int));
     512    whenUsed_ = new int[2 * ncols + 2 * nrows + extra];
     513    int *posSlack = whenUsed_ + ncols;
     514    int *negSlack = posSlack + nrows;
     515    int *nextSlack = negSlack + nrows;
     516    for (i = 0; i < nrows; i++) {
     517      posSlack[i] = -1;
     518      negSlack[i] = -1;
     519    }
     520    for (i = 0; i < ncols; i++)
     521      nextSlack[i] = -1;
     522    for (i = 0; i < ncols; i++) {
     523      if (columnLength[i] == 1) {
     524        CoinBigIndex j = columnStart[i];
     525        int iRow = row[j];
     526        if (element[j] > 0.0) {
     527          if (posSlack[iRow] == -1) {
     528            posSlack[iRow] = i;
     529          } else {
     530            int iCol = posSlack[iRow];
     531            while (nextSlack[iCol] >= 0)
     532              iCol = nextSlack[iCol];
     533            nextSlack[iCol] = i;
     534          }
     535        } else {
     536          if (negSlack[iRow] == -1) {
     537            negSlack[iRow] = i;
     538          } else {
     539            int iCol = negSlack[iRow];
     540            while (nextSlack[iCol] >= 0)
     541              iCol = nextSlack[iCol];
     542            nextSlack[iCol] = i;
     543          }
     544        }
     545      }
     546    }
     547    // now sort
     548    for (i = 0; i < nrows; i++) {
     549      int iCol;
     550      iCol = posSlack[i];
     551      if (iCol >= 0) {
     552        CoinBigIndex j = columnStart[iCol];
    557553#ifndef NDEBUG
    558                     int iRow = row[j];
    559 #endif
    560                     assert (element[j] > 0.0);
    561                     assert (iRow == i);
    562                     dj[0] = cost[iCol] / element[j];
    563                     whenUsed_[0] = iCol;
    564                     int n = 1;
    565                     while (nextSlack[iCol] >= 0) {
    566                          iCol = nextSlack[iCol];
    567                          CoinBigIndex j = columnStart[iCol];
     554        int iRow = row[j];
     555#endif
     556        assert(element[j] > 0.0);
     557        assert(iRow == i);
     558        dj[0] = cost[iCol] / element[j];
     559        whenUsed_[0] = iCol;
     560        int n = 1;
     561        while (nextSlack[iCol] >= 0) {
     562          iCol = nextSlack[iCol];
     563          CoinBigIndex j = columnStart[iCol];
    568564#ifndef NDEBUG
    569                          int iRow = row[j];
    570 #endif
    571                          assert (element[j] > 0.0);
    572                          assert (iRow == i);
    573                          dj[n] = cost[iCol] / element[j];
    574                          whenUsed_[n++] = iCol;
    575                     }
    576                     for (j = 0; j < n; j++) {
    577                          int jCol = whenUsed_[j];
    578                          nextSlack[jCol] = -2;
    579                     }
    580                     CoinSort_2(dj, dj + n, whenUsed_);
    581                     // put back
    582                     iCol = whenUsed_[0];
    583                     posSlack[i] = iCol;
    584                     for (j = 1; j < n; j++) {
    585                          int jCol = whenUsed_[j];
    586                          nextSlack[iCol] = jCol;
    587                          iCol = jCol;
    588                     }
    589                }
    590                iCol = negSlack[i];
    591                if (iCol >= 0) {
    592                     CoinBigIndex j = columnStart[iCol];
     565          int iRow = row[j];
     566#endif
     567          assert(element[j] > 0.0);
     568          assert(iRow == i);
     569          dj[n] = cost[iCol] / element[j];
     570          whenUsed_[n++] = iCol;
     571        }
     572        for (j = 0; j < n; j++) {
     573          int jCol = whenUsed_[j];
     574          nextSlack[jCol] = -2;
     575        }
     576        CoinSort_2(dj, dj + n, whenUsed_);
     577        // put back
     578        iCol = whenUsed_[0];
     579        posSlack[i] = iCol;
     580        for (j = 1; j < n; j++) {
     581          int jCol = whenUsed_[j];
     582          nextSlack[iCol] = jCol;
     583          iCol = jCol;
     584        }
     585      }
     586      iCol = negSlack[i];
     587      if (iCol >= 0) {
     588        CoinBigIndex j = columnStart[iCol];
    593589#ifndef NDEBUG
    594                     int iRow = row[j];
    595 #endif
    596                     assert (element[j] < 0.0);
    597                     assert (iRow == i);
    598                     dj[0] = -cost[iCol] / element[j];
    599                     whenUsed_[0] = iCol;
    600                     int n = 1;
    601                     while (nextSlack[iCol] >= 0) {
    602                          iCol = nextSlack[iCol];
    603                          CoinBigIndex j = columnStart[iCol];
     590        int iRow = row[j];
     591#endif
     592        assert(element[j] < 0.0);
     593        assert(iRow == i);
     594        dj[0] = -cost[iCol] / element[j];
     595        whenUsed_[0] = iCol;
     596        int n = 1;
     597        while (nextSlack[iCol] >= 0) {
     598          iCol = nextSlack[iCol];
     599          CoinBigIndex j = columnStart[iCol];
    604600#ifndef NDEBUG
    605                          int iRow = row[j];
    606 #endif
    607                          assert (element[j] < 0.0);
    608                          assert (iRow == i);
    609                          dj[n] = -cost[iCol] / element[j];
    610                          whenUsed_[n++] = iCol;
    611                     }
    612                     for (j = 0; j < n; j++) {
    613                          int jCol = whenUsed_[j];
    614                          nextSlack[jCol] = -2;
    615                     }
    616                     CoinSort_2(dj, dj + n, whenUsed_);
    617                     // put back
    618                     iCol = whenUsed_[0];
    619                     negSlack[i] = iCol;
    620                     for (j = 1; j < n; j++) {
    621                          int jCol = whenUsed_[j];
    622                          nextSlack[iCol] = jCol;
    623                          iCol = jCol;
    624                     }
    625                }
    626           }
    627      }
    628      whenUsed = whenUsed_;
    629      if (model_->getObjSense() == -1.0) {
    630           maxmin = -1.0;
    631      } else {
    632           maxmin = 1.0;
    633      }
    634      model_->getDblParam(OsiObjOffset, offset);
    635      if (!maxIts2) maxIts2 = maxIts;
    636      strategy = strategy_;
    637      strategy &= 3;
    638      memset(lambda, 0, nrows * sizeof(double));
    639      slackStart = countCostedSlacks(model_);
    640      // redo in case odd matrix
    641      row = matrix->getIndices();
    642      columnStart = matrix->getVectorStarts();
    643      columnLength = matrix->getVectorLengths();
    644      element = matrix->getElements();
    645      if (slackStart >= 0) {
    646        COIN_DETAIL_PRINT(printf("This model has costed slacks\n"));
    647           if (slackStart) {
    648                ordStart = 0;
    649                ordEnd = slackStart;
    650           } else {
    651                ordStart = nrows;
    652                ordEnd = ncols;
    653           }
    654      } else {
    655           ordStart = 0;
    656           ordEnd = ncols;
    657      }
    658      if (offset && logLevel > 2) {
    659           printf("** Objective offset is %g\n", offset);
    660      }
    661      /* compute reasonable solution cost */
    662      for (i = 0; i < nrows; i++) {
    663           rowsol[i] = 1.0e31;
    664      }
    665      for (i = 0; i < ncols; i++) {
    666           CoinBigIndex j;
    667           for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    668                if (element[j] != 1.0) {
    669                     allOnes = 0;
    670                     break;
    671                }
    672           }
    673      }
    674      if (allOnes) {
    675           elemXX = NULL;
    676      } else {
    677           elemXX = element;
    678      }
    679      // Do scaling if wanted
    680      bool scaled = false;
     601          int iRow = row[j];
     602#endif
     603          assert(element[j] < 0.0);
     604          assert(iRow == i);
     605          dj[n] = -cost[iCol] / element[j];
     606          whenUsed_[n++] = iCol;
     607        }
     608        for (j = 0; j < n; j++) {
     609          int jCol = whenUsed_[j];
     610          nextSlack[jCol] = -2;
     611        }
     612        CoinSort_2(dj, dj + n, whenUsed_);
     613        // put back
     614        iCol = whenUsed_[0];
     615        negSlack[i] = iCol;
     616        for (j = 1; j < n; j++) {
     617          int jCol = whenUsed_[j];
     618          nextSlack[iCol] = jCol;
     619          iCol = jCol;
     620        }
     621      }
     622    }
     623  }
     624  whenUsed = whenUsed_;
     625  if (model_->getObjSense() == -1.0) {
     626    maxmin = -1.0;
     627  } else {
     628    maxmin = 1.0;
     629  }
     630  model_->getDblParam(OsiObjOffset, offset);
     631  if (!maxIts2)
     632    maxIts2 = maxIts;
     633  strategy = strategy_;
     634  strategy &= 3;
     635  memset(lambda, 0, nrows * sizeof(double));
     636  slackStart = countCostedSlacks(model_);
     637  // redo in case odd matrix
     638  row = matrix->getIndices();
     639  columnStart = matrix->getVectorStarts();
     640  columnLength = matrix->getVectorLengths();
     641  element = matrix->getElements();
     642  if (slackStart >= 0) {
     643    COIN_DETAIL_PRINT(printf("This model has costed slacks\n"));
     644    if (slackStart) {
     645      ordStart = 0;
     646      ordEnd = slackStart;
     647    } else {
     648      ordStart = nrows;
     649      ordEnd = ncols;
     650    }
     651  } else {
     652    ordStart = 0;
     653    ordEnd = ncols;
     654  }
     655  if (offset && logLevel > 2) {
     656    printf("** Objective offset is %g\n", offset);
     657  }
     658  /* compute reasonable solution cost */
     659  for (i = 0; i < nrows; i++) {
     660    rowsol[i] = 1.0e31;
     661  }
     662  for (i = 0; i < ncols; i++) {
     663    CoinBigIndex j;
     664    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     665      if (element[j] != 1.0) {
     666        allOnes = 0;
     667        break;
     668      }
     669    }
     670  }
     671  if (allOnes) {
     672    elemXX = NULL;
     673  } else {
     674    elemXX = element;
     675  }
     676  // Do scaling if wanted
     677  bool scaled = false;
    681678#ifndef OSI_IDIOT
    682      if ((strategy_ & 32) != 0 && !allOnes) {
    683           if (model_->scalingFlag() > 0)
    684                scaled = model_->clpMatrix()->scale(model_) == 0;
    685           if (scaled) {
     679  if ((strategy_ & 32) != 0 && !allOnes) {
     680    if (model_->scalingFlag() > 0)
     681      scaled = model_->clpMatrix()->scale(model_) == 0;
     682    if (scaled) {
    686683#define IDIOT_SCALE 2
    687684#ifndef IDIOT_SCALE
    688                const double * rowScale = model_->rowScale();
     685      const double *rowScale = model_->rowScale();
    689686#else
    690                double * rowScale = model_->mutableRowScale();
    691 #endif
    692                const double * columnScale = model_->columnScale();
    693                double * oldLower = lower;
    694                double * oldUpper = upper;
    695                double * oldCost = cost;
    696                lower = new double[ncols];
    697                upper = new double[ncols];
    698                cost = new double[ncols];
    699                CoinMemcpyN(oldLower, ncols, lower);
    700                CoinMemcpyN(oldUpper, ncols, upper);
    701                CoinMemcpyN(oldCost, ncols, cost);
    702                int icol, irow;
    703 #if IDIOT_SCALE<0
    704                for (irow = 0; irow < nrows; irow++) {
    705                  rowlower[irow]=1.0e100;
    706                  rowupper[irow]=1.0e-100;
    707                }
    708 #endif
    709                for (icol = 0; icol < ncols; icol++) {
    710                     double multiplier = 1.0 / columnScale[icol];
    711                     if (lower[icol] > -1.0e50)
    712                          lower[icol] *= multiplier;
    713                     if (upper[icol] < 1.0e50)
    714                          upper[icol] *= multiplier;
    715                     colsol[icol] *= multiplier;
    716                     cost[icol] *= columnScale[icol];
    717 #if IDIOT_SCALE<0
    718                     CoinBigIndex j;
    719                     double scale = columnScale[i];
    720                     for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    721                          int jrow = row[j];
    722                          double scaledValue = fabs(scale*element[j]);
    723                          rowlower[jrow]=CoinMin(rowlower[jrow],scaledValue);
    724                          rowupper[jrow]=CoinMax(rowupper[jrow],scaledValue);
    725                     }
    726 #endif
    727                }
     687      double *rowScale = model_->mutableRowScale();
     688#endif
     689      const double *columnScale = model_->columnScale();
     690      double *oldLower = lower;
     691      double *oldUpper = upper;
     692      double *oldCost = cost;
     693      lower = new double[ncols];
     694      upper = new double[ncols];
     695      cost = new double[ncols];
     696      CoinMemcpyN(oldLower, ncols, lower);
     697      CoinMemcpyN(oldUpper, ncols, upper);
     698      CoinMemcpyN(oldCost, ncols, cost);
     699      int icol, irow;
     700#if IDIOT_SCALE < 0
     701      for (irow = 0; irow < nrows; irow++) {
     702        rowlower[irow] = 1.0e100;
     703        rowupper[irow] = 1.0e-100;
     704      }
     705#endif
     706      for (icol = 0; icol < ncols; icol++) {
     707        double multiplier = 1.0 / columnScale[icol];
     708        if (lower[icol] > -1.0e50)
     709          lower[icol] *= multiplier;
     710        if (upper[icol] < 1.0e50)
     711          upper[icol] *= multiplier;
     712        colsol[icol] *= multiplier;
     713        cost[icol] *= columnScale[icol];
     714#if IDIOT_SCALE < 0
     715        CoinBigIndex j;
     716        double scale = columnScale[i];
     717        for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     718          int jrow = row[j];
     719          double scaledValue = fabs(scale * element[j]);
     720          rowlower[jrow] = CoinMin(rowlower[jrow], scaledValue);
     721          rowupper[jrow] = CoinMax(rowupper[jrow], scaledValue);
     722        }
     723#endif
     724      }
    728725#ifdef IDIOT_SCALE
    729 #if IDIOT_SCALE>1||IDIOT_SCALE<-1
    730                const double * rowLower=model_->rowLower();
    731                const double * rowUpper=model_->rowUpper();
    732 #endif
    733                for (irow = 0; irow < nrows; irow++) {
    734 #if IDIOT_SCALE<0
    735                  double multiplier = 1.0/sqrt(rowlower[irow]*rowupper[irow]);
     726#if IDIOT_SCALE > 1 || IDIOT_SCALE < -1
     727      const double *rowLower = model_->rowLower();
     728      const double *rowUpper = model_->rowUpper();
     729#endif
     730      for (irow = 0; irow < nrows; irow++) {
     731#if IDIOT_SCALE < 0
     732        double multiplier = 1.0 / sqrt(rowlower[irow] * rowupper[irow]);
    736733#else
    737                 double multiplier = rowScale[irow];
    738 #endif
    739 #if IDIOT_SCALE>1||IDIOT_SCALE<-1
     734        double multiplier = rowScale[irow];
     735#endif
     736#if IDIOT_SCALE > 1 || IDIOT_SCALE < -1
    740737#define EQUALITY_MULTIPLIER 2
    741                  if (rowLower[irow]==rowUpper[irow])
    742                    multiplier *= EQUALITY_MULTIPLIER;
    743 #if IDIOT_SCALE>2||IDIOT_SCALE<-2
    744                  if (rowLower[irow]==rowUpper[irow]&&!rowlower[irow])
    745                    multiplier *= EQUALITY_MULTIPLIER;
    746 #endif
    747 #endif
    748                  rowScale[irow]=multiplier;
    749                }
    750                CoinMemcpyN(model_->rowUpper(), nrows, rowupper);
    751 #endif
    752                CoinMemcpyN(model_->rowLower(), nrows, rowlower);
    753                for (irow = 0; irow < nrows; irow++) {
    754                     double multiplier = rowScale[irow];
    755                     if (rowlower[irow] > -1.0e50)
    756                          rowlower[irow] *= multiplier;
    757                     if (rowupper[irow] < 1.0e50)
    758                          rowupper[irow] *= multiplier;
    759                     rowsol[irow] *= multiplier;
    760                }
    761                CoinBigIndex length = columnStart[ncols-1] + columnLength[ncols-1];
    762                double * elemYY = new double[length];
    763                for (i = 0; i < ncols; i++) {
    764                     CoinBigIndex j;
    765                     double scale = columnScale[i];
    766                     for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    767                          int irow = row[j];
    768                          elemYY[j] = element[j] * scale * rowScale[irow];
    769                     }
    770                }
    771                elemXX = elemYY;
    772           }
    773      }
    774 #endif
    775      if ((strategy_&131072)!=0) {
    776        // going to mess with cost
    777        if (cost==model_->objective())
    778          cost=CoinCopyOfArray(cost,ncols);
    779      }
    780      for (i = 0; i < ncols; i++) {
    781           CoinBigIndex j;
    782           double dd = columnLength[i];
    783           dd = cost[i] / dd;
    784           for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    785                int irow = row[j];
    786                if (dd < rowsol[irow]) {
    787                     rowsol[irow] = dd;
    788                }
    789           }
    790      }
    791      d2 = 0.0;
    792      for (i = 0; i < nrows; i++) {
    793           d2 += rowsol[i];
    794      }
    795      d2 *= 2.0; /* for luck */
     738        if (rowLower[irow] == rowUpper[irow])
     739          multiplier *= EQUALITY_MULTIPLIER;
     740#if IDIOT_SCALE > 2 || IDIOT_SCALE < -2
     741        if (rowLower[irow] == rowUpper[irow] && !rowlower[irow])
     742          multiplier *= EQUALITY_MULTIPLIER;
     743#endif
     744#endif
     745        rowScale[irow] = multiplier;
     746      }
     747      CoinMemcpyN(model_->rowUpper(), nrows, rowupper);
     748#endif
     749      CoinMemcpyN(model_->rowLower(), nrows, rowlower);
     750      for (irow = 0; irow < nrows; irow++) {
     751        double multiplier = rowScale[irow];
     752        if (rowlower[irow] > -1.0e50)
     753          rowlower[irow] *= multiplier;
     754        if (rowupper[irow] < 1.0e50)
     755          rowupper[irow] *= multiplier;
     756        rowsol[irow] *= multiplier;
     757      }
     758      CoinBigIndex length = columnStart[ncols - 1] + columnLength[ncols - 1];
     759      double *elemYY = new double[length];
     760      for (i = 0; i < ncols; i++) {
     761        CoinBigIndex j;
     762        double scale = columnScale[i];
     763        for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     764          int irow = row[j];
     765          elemYY[j] = element[j] * scale * rowScale[irow];
     766        }
     767      }
     768      elemXX = elemYY;
     769    }
     770  }
     771#endif
     772  if ((strategy_ & 131072) != 0) {
     773    // going to mess with cost
     774    if (cost == model_->objective())
     775      cost = CoinCopyOfArray(cost, ncols);
     776  }
     777  for (i = 0; i < ncols; i++) {
     778    CoinBigIndex j;
     779    double dd = columnLength[i];
     780    dd = cost[i] / dd;
     781    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     782      int irow = row[j];
     783      if (dd < rowsol[irow]) {
     784        rowsol[irow] = dd;
     785      }
     786    }
     787  }
     788  d2 = 0.0;
     789  for (i = 0; i < nrows; i++) {
     790    d2 += rowsol[i];
     791  }
     792  d2 *= 2.0; /* for luck */
    796793
    797      d2 = d2 / static_cast<double> (4 * nrows + 8000);
    798      d2 *= 0.5; /* halve with more flexible method */
    799      if (d2 < 5.0) d2 = 5.0;
    800      if (djExit == 0.0) {
    801           djExit = d2;
    802      }
    803      if ((saveStrategy & 4) != 0) {
    804           /* go to relative tolerances - first small */
    805           djExit = 1.0e-10;
    806           djFlag = 1.0e-5;
    807           drop = 1.0e-10;
    808      }
    809      memset(whenUsed, 0, ncols * sizeof(int));
    810      strategy = strategies[strategy];
    811      if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */
    812      CoinMemcpyN(colsol, ncols, saveSol);
     794  d2 = d2 / static_cast< double >(4 * nrows + 8000);
     795  d2 *= 0.5; /* halve with more flexible method */
     796  if (d2 < 5.0)
     797    d2 = 5.0;
     798  if (djExit == 0.0) {
     799    djExit = d2;
     800  }
     801  if ((saveStrategy & 4) != 0) {
     802    /* go to relative tolerances - first small */
     803    djExit = 1.0e-10;
     804    djFlag = 1.0e-5;
     805    drop = 1.0e-10;
     806  }
     807  memset(whenUsed, 0, ncols * sizeof(int));
     808  strategy = strategies[strategy];
     809  if ((saveStrategy & 8) != 0)
     810    strategy |= 64; /* don't allow large theta */
     811  CoinMemcpyN(colsol, ncols, saveSol);
    813812
    814      lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi,
    815                            dj, cost, rowlower, rowupper,
    816                            lower, upper, elemXX, row, columnStart, columnLength, lambda,
    817                            0, mu, drop,
    818                            maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
    819      // update whenUsed_
    820      double maxInfeasibility=COIN_DBL_MAX;
    821      n = cleanIteration(iteration, ordStart, ordEnd,
    822                         colsol,  lower, upper,
    823                         rowlower, rowupper,
    824                         cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
    825      if ((strategy_ & 16384) != 0) {
    826           int * posSlack = whenUsed_ + ncols;
    827           int * negSlack = posSlack + nrows;
    828           int * nextSlack = negSlack + nrows;
    829           double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
    830           for (i = 0; i < nrows; i++)
    831                rowsol[i] += rowsol2[i];
    832      }
    833      if ((logLevel_ & 1) != 0) {
     813  lastResult = IdiSolve(nrows, ncols, rowsol, colsol, pi,
     814    dj, cost, rowlower, rowupper,
     815    lower, upper, elemXX, row, columnStart, columnLength, lambda,
     816    0, mu, drop,
     817    maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
     818  // update whenUsed_
     819  double maxInfeasibility = COIN_DBL_MAX;
     820  n = cleanIteration(iteration, ordStart, ordEnd,
     821    colsol, lower, upper,
     822    rowlower, rowupper,
     823    cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
     824  if ((strategy_ & 16384) != 0) {
     825    int *posSlack = whenUsed_ + ncols;
     826    int *negSlack = posSlack + nrows;
     827    int *nextSlack = negSlack + nrows;
     828    double *rowsol2 = reinterpret_cast< double * >(nextSlack + ncols);
     829    for (i = 0; i < nrows; i++)
     830      rowsol[i] += rowsol2[i];
     831  }
     832  if ((logLevel_ & 1) != 0) {
    834833#ifndef OSI_IDIOT
    835           if (!handler) {
    836 #endif
    837                printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
    838                       iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n);
     834    if (!handler) {
     835#endif
     836      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
     837        iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n);
    839838#ifndef OSI_IDIOT
    840           } else {
    841                handler->message(CLP_IDIOT_ITERATION, *messages)
    842                          << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n
    843                          << CoinMessageEol;
    844           }
    845 #endif
    846      }
    847      int numberBaseTrys = 0; // for first time
    848      int numberAway = -1;
    849      iterationTotal = lastResult.iteration;
    850      firstInfeas = lastResult.infeas;
    851      if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas;
    852      if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas;
    853      double keepinfeas = 1.0e31;
    854      double lastInfeas = 1.0e31;
    855      double bestInfeas = 1.0e31;
    856      while ((mu > stopMu && lastResult.infeas > smallInfeas) ||
    857                (lastResult.infeas <= smallInfeas &&
    858                 dropping(lastResult, exitDrop, smallInfeas, &badIts)) ||
    859                checkIteration < 2 || lambdaIteration < lambdaIterations_) {
    860           if (lastResult.infeas <= exitFeasibility_)
    861                break;
    862           iteration++;
    863           //if (iteration>=exitAfter)
    864           //break;
    865           checkIteration++;
    866           if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
    867                bestFeasible = lastResult.objval;
    868           }
    869           if (lastResult.infeas + mu * lastResult.objval < bestWeighted) {
    870                bestWeighted = lastResult.objval + mu * lastResult.objval;
    871           }
    872           if ((saveStrategy & 4096)) strategy |= 256;
    873           if ((saveStrategy & 4) != 0 && iteration > 2) {
    874                /* go to relative tolerances */
    875                double weighted = 10.0 + fabs(lastWeighted);
    876                djExit = djTolerance_ * weighted;
    877                djFlag = 2.0 * djExit;
    878                drop = drop_ * weighted;
    879                djTol = 0.01 * djExit;
    880           }
    881           CoinMemcpyN(colsol, ncols, saveSol);
    882           result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj,
    883                             cost, rowlower, rowupper,
    884                             lower, upper, elemXX, row, columnStart, columnLength, lambda,
    885                             maxIts, mu, drop,
    886                             maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
    887           n = cleanIteration(iteration, ordStart, ordEnd,
    888                              colsol,  lower, upper,
    889                              rowlower, rowupper,
    890                              cost, elemXX, fixTolerance, result.objval, result.infeas, maxInfeasibility);
    891           if ((strategy_ & 16384) != 0) {
    892                int * posSlack = whenUsed_ + ncols;
    893                int * negSlack = posSlack + nrows;
    894                int * nextSlack = negSlack + nrows;
    895                double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
    896                for (i = 0; i < nrows; i++)
    897                     rowsol[i] += rowsol2[i];
    898           } else {
    899             maxInfeasibility=0.0;
    900             for (i = 0; i < nrows; i++) {
    901               double infeasibility = CoinMax(CoinMax(0.0, rowlower[i] - rowsol[i]), rowsol[i] - rowupper[i]);
    902               maxInfeasibility=CoinMax(maxInfeasibility,infeasibility);
    903             }
    904           }
    905           if ((logLevel_ & 1) != 0) {
     839    } else {
     840      handler->message(CLP_IDIOT_ITERATION, *messages)
     841        << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n
     842        << CoinMessageEol;
     843    }
     844#endif
     845  }
     846  int numberBaseTrys = 0; // for first time
     847  int numberAway = -1;
     848  iterationTotal = lastResult.iteration;
     849  firstInfeas = lastResult.infeas;
     850  if ((strategy_ & 1024) != 0)
     851    reasonableInfeas = 0.5 * firstInfeas;
     852  if (lastResult.infeas < reasonableInfeas)
     853    lastResult.infeas = reasonableInfeas;
     854  double keepinfeas = 1.0e31;
     855  double lastInfeas = 1.0e31;
     856  double bestInfeas = 1.0e31;
     857  while ((mu > stopMu && lastResult.infeas > smallInfeas) || (lastResult.infeas <= smallInfeas && dropping(lastResult, exitDrop, smallInfeas, &badIts)) || checkIteration < 2 || lambdaIteration < lambdaIterations_) {
     858    if (lastResult.infeas <= exitFeasibility_)
     859      break;
     860    iteration++;
     861    //if (iteration>=exitAfter)
     862    //break;
     863    checkIteration++;
     864    if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
     865      bestFeasible = lastResult.objval;
     866    }
     867    if (lastResult.infeas + mu * lastResult.objval < bestWeighted) {
     868      bestWeighted = lastResult.objval + mu * lastResult.objval;
     869    }
     870    if ((saveStrategy & 4096))
     871      strategy |= 256;
     872    if ((saveStrategy & 4) != 0 && iteration > 2) {
     873      /* go to relative tolerances */
     874      double weighted = 10.0 + fabs(lastWeighted);
     875      djExit = djTolerance_ * weighted;
     876      djFlag = 2.0 * djExit;
     877      drop = drop_ * weighted;
     878      djTol = 0.01 * djExit;
     879    }
     880    CoinMemcpyN(colsol, ncols, saveSol);
     881    result = IdiSolve(nrows, ncols, rowsol, colsol, pi, dj,
     882      cost, rowlower, rowupper,
     883      lower, upper, elemXX, row, columnStart, columnLength, lambda,
     884      maxIts, mu, drop,
     885      maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
     886    n = cleanIteration(iteration, ordStart, ordEnd,
     887      colsol, lower, upper,
     888      rowlower, rowupper,
     889      cost, elemXX, fixTolerance, result.objval, result.infeas, maxInfeasibility);
     890    if ((strategy_ & 16384) != 0) {
     891      int *posSlack = whenUsed_ + ncols;
     892      int *negSlack = posSlack + nrows;
     893      int *nextSlack = negSlack + nrows;
     894      double *rowsol2 = reinterpret_cast< double * >(nextSlack + ncols);
     895      for (i = 0; i < nrows; i++)
     896        rowsol[i] += rowsol2[i];
     897    } else {
     898      maxInfeasibility = 0.0;
     899      for (i = 0; i < nrows; i++) {
     900        double infeasibility = CoinMax(CoinMax(0.0, rowlower[i] - rowsol[i]), rowsol[i] - rowupper[i]);
     901        maxInfeasibility = CoinMax(maxInfeasibility, infeasibility);
     902      }
     903    }
     904    if ((logLevel_ & 1) != 0) {
    906905#ifndef OSI_IDIOT
    907                if (!handler) {
    908 #endif
    909                     printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
    910                            iteration, result.infeas, result.objval, mu, result.iteration, n);
     906      if (!handler) {
     907#endif
     908        printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
     909          iteration, result.infeas, result.objval, mu, result.iteration, n);
    911910#ifndef OSI_IDIOT
    912                } else {
    913                     handler->message(CLP_IDIOT_ITERATION, *messages)
    914                               << iteration << result.infeas << result.objval << mu << result.iteration << n
    915                               << CoinMessageEol;
    916                }
    917 #endif
    918           }
     911      } else {
     912        handler->message(CLP_IDIOT_ITERATION, *messages)
     913          << iteration << result.infeas << result.objval << mu << result.iteration << n
     914          << CoinMessageEol;
     915      }
     916#endif
     917    }
    919918#ifndef OSI_IDIOT
    920           if (fixAfterSome) {
    921             if (result.infeas<0.01*nrows&&iteration>10&&(3*n>2*nrows||4*n>2*ncols)) {
    922                 // flag
    923                 int numberColumns = model_->numberColumns();
    924                 printf("Flagging satisfied\n");
    925                 fixAfterSome=false;
    926                 for (int i=0;i<numberColumns;i++) {
    927                   if (colsol[i]>upper[i]-1.0e-7||
    928                       colsol[i]<lower[i]+1.0e-7)
    929                     model_->setColumnStatus(i,ClpSimplex::isFixed);
    930                 }
    931               }
    932               }
    933 #endif
    934           if (iteration > exitAfter ) {
    935             if((result.infeas < 1.0e-4 && majorIterations_<200 && n == numberAway)||result.infeas<1.0e-8||maxInfeasibility<1.0e-8) {
     919    if (fixAfterSome) {
     920      if (result.infeas < 0.01 * nrows && iteration > 10 && (3 * n > 2 * nrows || 4 * n > 2 * ncols)) {
     921        // flag
     922        int numberColumns = model_->numberColumns();
     923        printf("Flagging satisfied\n");
     924        fixAfterSome = false;
     925        for (int i = 0; i < numberColumns; i++) {
     926          if (colsol[i] > upper[i] - 1.0e-7 || colsol[i] < lower[i] + 1.0e-7)
     927            model_->setColumnStatus(i, ClpSimplex::isFixed);
     928        }
     929      }
     930    }
     931#endif
     932    if (iteration > exitAfter) {
     933      if ((result.infeas < 1.0e-4 && majorIterations_ < 200 && n == numberAway) || result.infeas < 1.0e-8 || maxInfeasibility < 1.0e-8) {
    936934#ifdef CLP_INVESTIGATE
    937                printf("infeas small %g\n", result.infeas);
    938 #endif
    939                if ((strategy_&131072)==0) {
    940                  break; // not much happening
    941                } else {
    942                  int n=0;
    943                  for (int i=0;i<ncols;i++) {
    944                    if (cost[i])
    945                      n++;
    946                  }
    947                  if (n*10<ncols) {
    948                    // fix ones with costs
    949                    exitAfter=100;
    950                    for (int i=0;i<ncols;i++) {
    951                      if (cost[i]||colsol[i]<lower[i]+1.0e-9||
    952                          colsol[i]>upper[i]-1.0e-9) {
    953                        cost[i]=0.0;
    954                        model_->setColumnStatus(i,ClpSimplex::isFixed);
    955                      } else {
    956                        cost[i]=0.0;
    957                        double gap = upper[i]-lower[i];
    958                        if (colsol[i]>upper[i]-0.25*gap) {
    959                          cost[i]=gap/(colsol[i]-upper[i]);
    960                        } else if (colsol[i]<lower[i]+0.25*gap) {
    961                          cost[i]=gap/(colsol[i]-lower[i]);
    962                        }
    963                      }
    964                    }
    965                  }
    966                }
    967             }
    968           }
    969           if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
    970                if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas)
    971                     majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible
    972           }
    973           lastInfeas = result.infeas;
    974           numberAway = n;
    975           keepinfeas = result.infeas;
    976           lastWeighted = result.weighted;
    977           iterationTotal += result.iteration;
    978           if (iteration == 1) {
    979                if ((strategy_ & 1024) != 0 && mu < 1.0e-10)
    980                     result.infeas = firstInfeas * 0.8;
    981                if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0)
    982                     result.infeas *= 0.8;
    983                if (result.infeas > firstInfeas * 0.9
    984                          && result.infeas > reasonableInfeas) {
    985                     iteration--;
    986                     if (majorIterations_ < 50)
    987                          mu *= 1.0e-1;
    988                     else
    989                          mu *= 0.7;
    990                     bestFeasible = 1.0e31;
    991                     bestWeighted = 1.0e60;
    992                     numberBaseTrys++;
    993                     if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) {
    994                          // back to all slack basis
    995                          lightWeight_ = 2;
    996                          break;
    997                     }
    998                     CoinMemcpyN(saveSol, ncols, colsol);
    999                } else {
    1000                     maxIts = maxIts2;
    1001                     checkIteration = 0;
    1002                     if ((strategy_ & 1024) != 0) mu *= 1.0e-1;
    1003                }
    1004           } else {
    1005           }
    1006           bestInfeas = CoinMin(bestInfeas, result.infeas);
    1007           if (majorIterations_>100&&majorIterations_<200) {
    1008             if (iteration==majorIterations_-100) {
    1009               // redo
    1010               double muX=mu*10.0;
    1011               bestInfeas=1.0e3;
    1012               mu=muX;
    1013               nTry=0;
    1014             }
    1015           }
    1016           if (iteration) {
    1017                /* this code is in to force it to terminate sometime */
    1018                double changeMu = factor;
    1019                if ((saveStrategy & 64) != 0) {
    1020                     keepinfeas = 0.0; /* switch off ranga's increase */
    1021                     fakeSmall = smallInfeas;
    1022                } else {
    1023                     fakeSmall = -1.0;
    1024                }
    1025                saveLambdaScale = 0.0;
    1026                if (result.infeas > reasonableInfeas ||
    1027                          (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) {
    1028                     if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) ||
    1029                               nTry + 1 == maxBigIts ||
    1030                               (result.infeas > lastResult.infeas * 0.9
    1031                                && result.weighted > lastResult.weighted
    1032                                - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) {
    1033                          mu *= changeMu;
    1034                          if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) {
    1035                               reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu));
    1036                               COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n", reasonableInfeas));
    1037                          }
    1038                          result.weighted = 1.0e60;
    1039                          nTry = 0;
    1040                          bestFeasible = 1.0e31;
    1041                          bestWeighted = 1.0e60;
    1042                          checkIteration = 0;
    1043                          lambdaIteration = 0;
     935        printf("infeas small %g\n", result.infeas);
     936#endif
     937        if ((strategy_ & 131072) == 0) {
     938          break; // not much happening
     939        } else {
     940          int n = 0;
     941          for (int i = 0; i < ncols; i++) {
     942            if (cost[i])
     943              n++;
     944          }
     945          if (n * 10 < ncols) {
     946            // fix ones with costs
     947            exitAfter = 100;
     948            for (int i = 0; i < ncols; i++) {
     949              if (cost[i] || colsol[i] < lower[i] + 1.0e-9 || colsol[i] > upper[i] - 1.0e-9) {
     950                cost[i] = 0.0;
     951                model_->setColumnStatus(i, ClpSimplex::isFixed);
     952              } else {
     953                cost[i] = 0.0;
     954                double gap = upper[i] - lower[i];
     955                if (colsol[i] > upper[i] - 0.25 * gap) {
     956                  cost[i] = gap / (colsol[i] - upper[i]);
     957                } else if (colsol[i] < lower[i] + 0.25 * gap) {
     958                  cost[i] = gap / (colsol[i] - lower[i]);
     959                }
     960              }
     961            }
     962          }
     963        }
     964      }
     965    }
     966    if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
     967      if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas)
     968        majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible
     969    }
     970    lastInfeas = result.infeas;
     971    numberAway = n;
     972    keepinfeas = result.infeas;
     973    lastWeighted = result.weighted;
     974    iterationTotal += result.iteration;
     975    if (iteration == 1) {
     976      if ((strategy_ & 1024) != 0 && mu < 1.0e-10)
     977        result.infeas = firstInfeas * 0.8;
     978      if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0)
     979        result.infeas *= 0.8;
     980      if (result.infeas > firstInfeas * 0.9
     981        && result.infeas > reasonableInfeas) {
     982        iteration--;
     983        if (majorIterations_ < 50)
     984          mu *= 1.0e-1;
     985        else
     986          mu *= 0.7;
     987        bestFeasible = 1.0e31;
     988        bestWeighted = 1.0e60;
     989        numberBaseTrys++;
     990        if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) {
     991          // back to all slack basis
     992          lightWeight_ = 2;
     993          break;
     994        }
     995        CoinMemcpyN(saveSol, ncols, colsol);
     996      } else {
     997        maxIts = maxIts2;
     998        checkIteration = 0;
     999        if ((strategy_ & 1024) != 0)
     1000          mu *= 1.0e-1;
     1001      }
     1002    } else {
     1003    }
     1004    bestInfeas = CoinMin(bestInfeas, result.infeas);
     1005    if (majorIterations_ > 100 && majorIterations_ < 200) {
     1006      if (iteration == majorIterations_ - 100) {
     1007        // redo
     1008        double muX = mu * 10.0;
     1009        bestInfeas = 1.0e3;
     1010        mu = muX;
     1011        nTry = 0;
     1012      }
     1013    }
     1014    if (iteration) {
     1015      /* this code is in to force it to terminate sometime */
     1016      double changeMu = factor;
     1017      if ((saveStrategy & 64) != 0) {
     1018        keepinfeas = 0.0; /* switch off ranga's increase */
     1019        fakeSmall = smallInfeas;
     1020      } else {
     1021        fakeSmall = -1.0;
     1022      }
     1023      saveLambdaScale = 0.0;
     1024      if (result.infeas > reasonableInfeas || (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) {
     1025        if (result.infeas > lastResult.infeas * (1.0 - dropEnoughFeasibility_) || nTry + 1 == maxBigIts || (result.infeas > lastResult.infeas * 0.9 && result.weighted > lastResult.weighted - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) {
     1026          mu *= changeMu;
     1027          if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) {
     1028            reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu));
     1029            COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n", reasonableInfeas));
     1030          }
     1031          result.weighted = 1.0e60;
     1032          nTry = 0;
     1033          bestFeasible = 1.0e31;
     1034          bestWeighted = 1.0e60;
     1035          checkIteration = 0;
     1036          lambdaIteration = 0;
    10441037#define LAMBDA
    10451038#ifdef LAMBDA
    1046                          if ((saveStrategy & 2048) == 0) {
    1047                               memset(lambda, 0, nrows * sizeof(double));
    1048                          }
     1039          if ((saveStrategy & 2048) == 0) {
     1040            memset(lambda, 0, nrows * sizeof(double));
     1041          }
    10491042#else
    1050                          memset(lambda, 0, nrows * sizeof(double));
    1051 #endif
    1052                     } else {
    1053                          nTry++;
    1054                     }
    1055                } else if (lambdaIterations_ >= 0) {
    1056                     /* update lambda  */
    1057                     double scale = 1.0 / mu;
    1058                     int i, nnz = 0;
    1059                     saveLambdaScale = scale;
    1060                     lambdaIteration++;
    1061                     if ((saveStrategy & 4) == 0) drop = drop_ / 50.0;
    1062                     if (lambdaIteration > 4 &&
    1063                               (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) ||
    1064                                ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) {
    1065                          //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
    1066                          smallInfeas *= 1.5;
    1067                     }
    1068                     if ((saveStrategy & 2048) == 0) {
    1069                          for (i = 0; i < nrows; i++) {
    1070                               if (lambda[i]) nnz++;
    1071                               lambda[i] += scale * rowsol[i];
    1072                          }
    1073                     } else {
    1074                          nnz = 1;
     1043          memset(lambda, 0, nrows * sizeof(double));
     1044#endif
     1045        } else {
     1046          nTry++;
     1047        }
     1048      } else if (lambdaIterations_ >= 0) {
     1049        /* update lambda  */
     1050        double scale = 1.0 / mu;
     1051        int i, nnz = 0;
     1052        saveLambdaScale = scale;
     1053        lambdaIteration++;
     1054        if ((saveStrategy & 4) == 0)
     1055          drop = drop_ / 50.0;
     1056        if (lambdaIteration > 4 && (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) || ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) {
     1057          //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
     1058          smallInfeas *= 1.5;
     1059        }
     1060        if ((saveStrategy & 2048) == 0) {
     1061          for (i = 0; i < nrows; i++) {
     1062            if (lambda[i])
     1063              nnz++;
     1064            lambda[i] += scale * rowsol[i];
     1065          }
     1066        } else {
     1067          nnz = 1;
    10751068#ifdef LAMBDA
    1076                          for (i = 0; i < nrows; i++) {
    1077                               lambda[i] += scale * rowsol[i];
    1078                          }
     1069          for (i = 0; i < nrows; i++) {
     1070            lambda[i] += scale * rowsol[i];
     1071          }
    10791072#else
    1080                          for (i = 0; i < nrows; i++) {
    1081                               lambda[i] = scale * rowsol[i];
    1082                          }
    1083                          for (i = 0; i < ncols; i++) {
    1084                               CoinBigIndex j;
    1085                               double value = cost[i] * maxmin;
    1086                               for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1087                                    int irow = row[j];
    1088                                    value += element[j] * lambda[irow];
    1089                               }
    1090                               cost[i] = value * maxmin;
    1091                          }
    1092                          for (i = 0; i < nrows; i++) {
    1093                               offset += lambda[i] * rowupper[i];
    1094                               lambda[i] = 0.0;
    1095                          }
     1073          for (i = 0; i < nrows; i++) {
     1074            lambda[i] = scale * rowsol[i];
     1075          }
     1076          for (i = 0; i < ncols; i++) {
     1077            CoinBigIndex j;
     1078            double value = cost[i] * maxmin;
     1079            for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     1080              int irow = row[j];
     1081              value += element[j] * lambda[irow];
     1082            }
     1083            cost[i] = value * maxmin;
     1084          }
     1085          for (i = 0; i < nrows; i++) {
     1086            offset += lambda[i] * rowupper[i];
     1087            lambda[i] = 0.0;
     1088          }
    10961089#ifdef DEBUG
    1097                          printf("offset %g\n", offset);
    1098 #endif
    1099                          model_->setDblParam(OsiObjOffset, offset);
    1100 #endif
    1101                     }
    1102                     nTry++;
    1103                     if (!nnz) {
    1104                          bestFeasible = 1.0e32;
    1105                          bestWeighted = 1.0e60;
    1106                          checkIteration = 0;
    1107                          result.weighted = 1.0e31;
    1108                     }
     1090          printf("offset %g\n", offset);
     1091#endif
     1092          model_->setDblParam(OsiObjOffset, offset);
     1093#endif
     1094        }
     1095        nTry++;
     1096        if (!nnz) {
     1097          bestFeasible = 1.0e32;
     1098          bestWeighted = 1.0e60;
     1099          checkIteration = 0;
     1100          result.weighted = 1.0e31;
     1101        }
    11091102#ifdef DEBUG
    1110                     double trueCost = 0.0;
    1111                     for (i = 0; i < ncols; i++) {
    1112                          int j;
    1113                          trueCost += cost[i] * colsol[i];
    1114                     }
    1115                     printf("True objective %g\n", trueCost - offset);
    1116 #endif
    1117                } else {
    1118                     nTry++;
    1119                }
    1120                lastResult = result;
    1121                if (result.infeas < reasonableInfeas && !belowReasonable) {
    1122                     belowReasonable = 1;
    1123                     bestFeasible = 1.0e32;
    1124                     bestWeighted = 1.0e60;
    1125                     checkIteration = 0;
    1126                     result.weighted = 1.0e31;
    1127                }
    1128           }
    1129           if (iteration >= majorIterations_) {
    1130                // If not feasible and crash then dive dive dive
    1131                if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) {
    1132                     mu = 1.0e-30;
    1133                     majorIterations_ = iteration + 1;
    1134                     stopMu = 0.0;
    1135                } else {
    1136                     if (logLevel > 2)
    1137                          printf("Exiting due to number of major iterations\n");
    1138                     break;
    1139                }
    1140           }
    1141      }
    1142      majorIterations_ = saveMajorIterations;
     1103        double trueCost = 0.0;
     1104        for (i = 0; i < ncols; i++) {
     1105          int j;
     1106          trueCost += cost[i] * colsol[i];
     1107        }
     1108        printf("True objective %g\n", trueCost - offset);
     1109#endif
     1110      } else {
     1111        nTry++;
     1112      }
     1113      lastResult = result;
     1114      if (result.infeas < reasonableInfeas && !belowReasonable) {
     1115        belowReasonable = 1;
     1116        bestFeasible = 1.0e32;
     1117        bestWeighted = 1.0e60;
     1118        checkIteration = 0;
     1119        result.weighted = 1.0e31;
     1120      }
     1121    }
     1122    if (iteration >= majorIterations_) {
     1123      // If not feasible and crash then dive dive dive
     1124      if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) {
     1125        mu = 1.0e-30;
     1126        majorIterations_ = iteration + 1;
     1127        stopMu = 0.0;
     1128      } else {
     1129        if (logLevel > 2)
     1130          printf("Exiting due to number of major iterations\n");
     1131        break;
     1132      }
     1133    }
     1134  }
     1135  majorIterations_ = saveMajorIterations;
    11431136#ifndef OSI_IDIOT
    1144      if (scaled) {
    1145           // Scale solution and free arrays
    1146           const double * rowScale = model_->rowScale();
    1147           const double * columnScale = model_->columnScale();
    1148           int icol, irow;
    1149           for (icol = 0; icol < ncols; icol++) {
    1150                colsol[icol] *= columnScale[icol];
    1151                saveSol[icol] *= columnScale[icol];
    1152                dj[icol] /= columnScale[icol];
    1153           }
    1154           for (irow = 0; irow < nrows; irow++) {
    1155                rowsol[irow] /= rowScale[irow];
    1156                pi[irow] *= rowScale[irow];
    1157           }
    1158           // Don't know why getting Microsoft problems
    1159 #if defined (_MSC_VER)
    1160           delete [] ( double *) elemXX;
     1137  if (scaled) {
     1138    // Scale solution and free arrays
     1139    const double *rowScale = model_->rowScale();
     1140    const double *columnScale = model_->columnScale();
     1141    int icol, irow;
     1142    for (icol = 0; icol < ncols; icol++) {
     1143      colsol[icol] *= columnScale[icol];
     1144      saveSol[icol] *= columnScale[icol];
     1145      dj[icol] /= columnScale[icol];
     1146    }
     1147    for (irow = 0; irow < nrows; irow++) {
     1148      rowsol[irow] /= rowScale[irow];
     1149      pi[irow] *= rowScale[irow];
     1150    }
     1151    // Don't know why getting Microsoft problems
     1152#if defined(_MSC_VER)
     1153    delete[](double *) elemXX;
    11611154#else
    1162           delete [] elemXX;
    1163 #endif
    1164           model_->setRowScale(NULL);
    1165           model_->setColumnScale(NULL);
    1166           delete [] lower;
    1167           delete [] upper;
    1168           delete [] cost;
    1169           lower = model_->columnLower();
    1170           upper = model_->columnUpper();
    1171           cost = model_->objective();
    1172           //rowlower = model_->rowLower();
    1173      } else if (cost != model_->objective()) {
    1174        delete [] cost;
    1175        cost = model_->objective();
    1176      }
     1155    delete[] elemXX;
     1156#endif
     1157    model_->setRowScale(NULL);
     1158    model_->setColumnScale(NULL);
     1159    delete[] lower;
     1160    delete[] upper;
     1161    delete[] cost;
     1162    lower = model_->columnLower();
     1163    upper = model_->columnUpper();
     1164    cost = model_->objective();
     1165    //rowlower = model_->rowLower();
     1166  } else if (cost != model_->objective()) {
     1167    delete[] cost;
     1168    cost = model_->objective();
     1169  }
    11771170#endif
    11781171#define TRYTHIS
    11791172#ifdef TRYTHIS
    1180      if ((saveStrategy & 2048) != 0) {
    1181           double offset;
    1182           model_->getDblParam(OsiObjOffset, offset);
    1183           for (i = 0; i < ncols; i++) {
    1184                CoinBigIndex j;
    1185                double djval = cost[i] * maxmin;
    1186                for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1187                     int irow = row[j];
    1188                     djval -= element[j] * lambda[irow];
    1189                }
    1190                cost[i] = djval;
    1191           }
    1192           for (i = 0; i < nrows; i++) {
    1193                offset += lambda[i] * rowupper[i];
    1194           }
    1195           model_->setDblParam(OsiObjOffset, offset);
    1196      }
    1197 #endif
    1198      if (saveLambdaScale) {
    1199           /* back off last update */
    1200           for (i = 0; i < nrows; i++) {
    1201                lambda[i] -= saveLambdaScale * rowsol[i];
    1202           }
    1203      }
    1204      muAtExit_ = mu;
    1205      // For last iteration make as feasible as possible
    1206      if (oddSlacks)
    1207           strategy_ |= 16384;
    1208      // not scaled
    1209      n = cleanIteration(iteration, ordStart, ordEnd,
    1210                         colsol,  lower, upper,
    1211                         model_->rowLower(), model_->rowUpper(),
    1212                         cost, element, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
     1173  if ((saveStrategy & 2048) != 0) {
     1174    double offset;
     1175    model_->getDblParam(OsiObjOffset, offset);
     1176    for (i = 0; i < ncols; i++) {
     1177      CoinBigIndex j;
     1178      double djval = cost[i] * maxmin;
     1179      for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     1180        int irow = row[j];
     1181        djval -= element[j] * lambda[irow];
     1182      }
     1183      cost[i] = djval;
     1184    }
     1185    for (i = 0; i < nrows; i++) {
     1186      offset += lambda[i] * rowupper[i];
     1187    }
     1188    model_->setDblParam(OsiObjOffset, offset);
     1189  }
     1190#endif
     1191  if (saveLambdaScale) {
     1192    /* back off last update */
     1193    for (i = 0; i < nrows; i++) {
     1194      lambda[i] -= saveLambdaScale * rowsol[i];
     1195    }
     1196  }
     1197  muAtExit_ = mu;
     1198  // For last iteration make as feasible as possible
     1199  if (oddSlacks)
     1200    strategy_ |= 16384;
     1201  // not scaled
     1202  n = cleanIteration(iteration, ordStart, ordEnd,
     1203    colsol, lower, upper,
     1204    model_->rowLower(), model_->rowUpper(),
     1205    cost, element, fixTolerance, lastResult.objval, lastResult.infeas, maxInfeasibility);
    12131206#if 0
    12141207     if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) {
     
    12191212#endif
    12201213#ifndef OSI_IDIOT
    1221      model_->setSumPrimalInfeasibilities(lastResult.infeas);
    1222 #endif
    1223      // Put back more feasible solution
    1224      double saveInfeas[] = {0.0, 0.0};
    1225      for (int iSol = 0; iSol < 3; iSol++) {
    1226           const double * solution = iSol ? colsol : saveSol;
    1227           if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) {
    1228                // put back best solution
    1229                CoinMemcpyN(saveSol, ncols, colsol);
    1230           }
    1231           double large = 0.0;
    1232           int i;
    1233           memset(rowsol, 0, nrows * sizeof(double));
    1234           for (i = 0; i < ncols; i++) {
    1235                CoinBigIndex j;
    1236                double value = solution[i];
    1237                for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1238                     int irow = row[j];
    1239                     rowsol[irow] += element[j] * value;
    1240                }
    1241           }
    1242           for (i = 0; i < nrows; i++) {
    1243                if (rowsol[i] > rowupper[i]) {
    1244                     double diff = rowsol[i] - rowupper[i];
    1245                     if (diff > large)
    1246                          large = diff;
    1247                } else if (rowsol[i] < rowlower[i]) {
    1248                     double diff = rowlower[i] - rowsol[i];
    1249                     if (diff > large)
    1250                          large = diff;
    1251                }
    1252           }
    1253           if (iSol < 2)
    1254                saveInfeas[iSol] = large;
    1255           if (logLevel > 2)
    1256                printf("largest infeasibility is %g\n", large);
    1257      }
    1258      /* subtract out lambda */
    1259      for (i = 0; i < nrows; i++) {
    1260           pi[i] -= lambda[i];
    1261      }
    1262      for (i = 0; i < ncols; i++) {
    1263           CoinBigIndex j;
    1264           double djval = cost[i] * maxmin;
    1265           for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1266                int irow = row[j];
    1267                djval -= element[j] * pi[irow];
    1268           }
    1269           dj[i] = djval;
    1270      }
    1271      if ((strategy_ & 1024) != 0) {
    1272           double ratio = static_cast<double> (ncols) / static_cast<double> (nrows);
    1273           COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n", ratio, lastResult.infeas / firstInfeas));
    1274           if (lastResult.infeas > 0.01 * firstInfeas * ratio) {
    1275                strategy_ &= (~1024);
    1276                COIN_DETAIL_PRINT(printf(" - layer off\n"));
    1277           } else {
    1278             COIN_DETAIL_PRINT(printf(" - layer on\n"));
    1279           }
    1280      }
    1281      delete [] saveSol;
    1282      delete [] lambda;
    1283      // save solution
    1284      // duals not much use - but save anyway
     1214  model_->setSumPrimalInfeasibilities(lastResult.infeas);
     1215#endif
     1216  // Put back more feasible solution
     1217  double saveInfeas[] = { 0.0, 0.0 };
     1218  for (int iSol = 0; iSol < 3; iSol++) {
     1219    const double *solution = iSol ? colsol : saveSol;
     1220    if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) {
     1221      // put back best solution
     1222      CoinMemcpyN(saveSol, ncols, colsol);
     1223    }
     1224    double large = 0.0;
     1225    int i;
     1226    memset(rowsol, 0, nrows * sizeof(double));
     1227    for (i = 0; i < ncols; i++) {
     1228      CoinBigIndex j;
     1229      double value = solution[i];
     1230      for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     1231        int irow = row[j];
     1232        rowsol[irow] += element[j] * value;
     1233      }
     1234    }
     1235    for (i = 0; i < nrows; i++) {
     1236      if (rowsol[i] > rowupper[i]) {
     1237        double diff = rowsol[i] - rowupper[i];
     1238        if (diff > large)
     1239          large = diff;
     1240      } else if (rowsol[i] < rowlower[i]) {
     1241        double diff = rowlower[i] - rowsol[i];
     1242        if (diff > large)
     1243          large = diff;
     1244      }
     1245    }
     1246    if (iSol < 2)
     1247      saveInfeas[iSol] = large;
     1248    if (logLevel > 2)
     1249      printf("largest infeasibility is %g\n", large);
     1250  }
     1251  /* subtract out lambda */
     1252  for (i = 0; i < nrows; i++) {
     1253    pi[i] -= lambda[i];
     1254  }
     1255  for (i = 0; i < ncols; i++) {
     1256    CoinBigIndex j;
     1257    double djval = cost[i] * maxmin;
     1258    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     1259      int irow = row[j];
     1260      djval -= element[j] * pi[irow];
     1261    }
     1262    dj[i] = djval;
     1263  }
     1264  if ((strategy_ & 1024) != 0) {
     1265    double ratio = static_cast< double >(ncols) / static_cast< double >(nrows);
     1266    COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n", ratio, lastResult.infeas / firstInfeas));
     1267    if (lastResult.infeas > 0.01 * firstInfeas * ratio) {
     1268      strategy_ &= (~1024);
     1269      COIN_DETAIL_PRINT(printf(" - layer off\n"));
     1270    } else {
     1271      COIN_DETAIL_PRINT(printf(" - layer on\n"));
     1272    }
     1273  }
     1274  delete[] saveSol;
     1275  delete[] lambda;
     1276  // save solution
     1277  // duals not much use - but save anyway
    12851278#ifndef OSI_IDIOT
    1286      CoinMemcpyN(rowsol, nrows, model_->primalRowSolution());
    1287      CoinMemcpyN(colsol, ncols, model_->primalColumnSolution());
    1288      CoinMemcpyN(pi, nrows, model_->dualRowSolution());
    1289      CoinMemcpyN(dj, ncols, model_->dualColumnSolution());
     1279  CoinMemcpyN(rowsol, nrows, model_->primalRowSolution());
     1280  CoinMemcpyN(colsol, ncols, model_->primalColumnSolution());
     1281  CoinMemcpyN(pi, nrows, model_->dualRowSolution());
     1282  CoinMemcpyN(dj, ncols, model_->dualColumnSolution());
    12901283#else
    1291      model_->setColSolution(colsol);
    1292      model_->setRowPrice(pi);
    1293      delete [] cost;
    1294 #endif
    1295      delete [] rowsol;
    1296      delete [] colsol;
    1297      delete [] pi;
    1298      delete [] dj;
    1299      delete [] rowlower;
    1300      delete [] rowupper;
    1301      return ;
     1284  model_->setColSolution(colsol);
     1285  model_->setRowPrice(pi);
     1286  delete[] cost;
     1287#endif
     1288  delete[] rowsol;
     1289  delete[] colsol;
     1290  delete[] pi;
     1291  delete[] dj;
     1292  delete[] rowlower;
     1293  delete[] rowupper;
     1294  return;
    13021295}
    13031296#ifndef OSI_IDIOT
    1304 void
    1305 Idiot::crossOver(int mode)
     1297void Idiot::crossOver(int mode)
    13061298{
    1307      if (lightWeight_ == 2) {
    1308           // total failure
    1309           model_->allSlackBasis();
    1310           return;
    1311      }
    1312      double fixTolerance = IDIOT_FIX_TOLERANCE;
     1299  if (lightWeight_ == 2) {
     1300    // total failure
     1301    model_->allSlackBasis();
     1302    return;
     1303  }
     1304  double fixTolerance = IDIOT_FIX_TOLERANCE;
    13131305#ifdef COIN_DEVELOP
    1314      double startTime = CoinCpuTime();
    1315 #endif
    1316      ClpSimplex * saveModel = NULL;
    1317      ClpMatrixBase * matrix = model_->clpMatrix();
    1318      const int * row = matrix->getIndices();
    1319      const CoinBigIndex * columnStart = matrix->getVectorStarts();
    1320      const int * columnLength = matrix->getVectorLengths();
    1321      const double * element = matrix->getElements();
    1322      const double * rowupper = model_->getRowUpper();
    1323      model_->eventHandler()->event(ClpEventHandler::startOfCrossover);
    1324      int nrows = model_->getNumRows();
    1325      int ncols = model_->getNumCols();
    1326      double * rowsol, * colsol;
    1327      // different for Osi
    1328      double * lower = model_->columnLower();
    1329      double * upper = model_->columnUpper();
    1330      const double * rowlower = model_->getRowLower();
    1331      int * whenUsed = whenUsed_;
    1332      rowsol = model_->primalRowSolution();
    1333      colsol = model_->primalColumnSolution();;
    1334      double * cost = model_->objective();
    1335      int slackEnd, ordStart, ordEnd;
    1336      int slackStart = countCostedSlacks(model_);
     1306  double startTime = CoinCpuTime();
     1307#endif
     1308  ClpSimplex *saveModel = NULL;
     1309  ClpMatrixBase *matrix = model_->clpMatrix();
     1310  const int *row = matrix->getIndices();
     1311  const CoinBigIndex *columnStart = matrix->getVectorStarts();
     1312  const int *columnLength = matrix->getVectorLengths();
     1313  const double *element = matrix->getElements();
     1314  const double *rowupper = model_->getRowUpper();
     1315  model_->eventHandler()->event(ClpEventHandler::startOfCrossover);
     1316  int nrows = model_->getNumRows();
     1317  int ncols = model_->getNumCols();
     1318  double *rowsol, *colsol;
     1319  // different for Osi
     1320  double *lower = model_->columnLower();
     1321  double *upper = model_->columnUpper();
     1322  const double *rowlower = model_->getRowLower();
     1323  int *whenUsed = whenUsed_;
     1324  rowsol = model_->primalRowSolution();
     1325  colsol = model_->primalColumnSolution();
     1326  ;
     1327  double *cost = model_->objective();
     1328  int slackEnd, ordStart, ordEnd;
     1329  int slackStart = countCostedSlacks(model_);
    13371330
    1338      int addAll = mode & 7;
    1339      int presolve = 0;
     1331  int addAll = mode & 7;
     1332  int presolve = 0;
    13401333
    1341      double djTolerance = djTolerance_;
    1342      if (djTolerance > 0.0 && djTolerance < 1.0)
    1343           djTolerance = 1.0;
    1344      int iteration;
    1345      int i, n = 0;
    1346      double ratio = 1.0;
    1347      double objValue = 0.0;
    1348      if ((strategy_ & 128) != 0) {
    1349           fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
    1350      }
    1351      if ((mode & 16) != 0 && addAll < 3) presolve = 1;
    1352      double * saveUpper = NULL;
    1353      double * saveLower = NULL;
    1354      double * saveRowUpper = NULL;
    1355      double * saveRowLower = NULL;
    1356      bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000);
    1357      if (addAll < 3) {
    1358           saveUpper = new double [ncols];
    1359           saveLower = new double [ncols];
    1360           CoinMemcpyN(upper, ncols, saveUpper);
    1361           CoinMemcpyN(lower, ncols, saveLower);
    1362           if (allowInfeasible) {
    1363                saveRowUpper = new double [nrows];
    1364                saveRowLower = new double [nrows];
    1365                CoinMemcpyN(rowupper, nrows, saveRowUpper);
    1366                CoinMemcpyN(rowlower, nrows, saveRowLower);
    1367                double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
    1368                fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas);
    1369           }
    1370      }
    1371      if (slackStart >= 0) {
    1372           slackEnd = slackStart + nrows;
    1373           if (slackStart) {
    1374                ordStart = 0;
    1375                ordEnd = slackStart;
    1376           } else {
    1377                ordStart = nrows;
    1378                ordEnd = ncols;
    1379           }
    1380      } else {
    1381           slackEnd = slackStart;
    1382           ordStart = 0;
    1383           ordEnd = ncols;
    1384      }
    1385      /* get correct rowsol (without known slacks) */
    1386      memset(rowsol, 0, nrows * sizeof(double));
    1387      for (i = ordStart; i < ordEnd; i++) {
    1388           CoinBigIndex j;
    1389           double value = colsol[i];
    1390           if (value < lower[i] + fixTolerance) {
    1391                value = lower[i];
    1392                colsol[i] = value;
    1393           }
    1394           for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
    1395                int irow = row[j];
    1396                rowsol[irow] += value * element[j];
    1397           }
    1398      }
    1399      if (slackStart >= 0) {
    1400           for (i = 0; i < nrows; i++) {
    1401                if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) {
    1402                     ratio = rowlower[i] / rowsol[i];
    1403                }
    1404           }
    1405           for (i = 0; i < nrows; i++) {
    1406                rowsol[i] *= ratio;
    1407           }
    1408           for (i = ordStart; i < ordEnd; i++) {
    1409                double value = colsol[i] * ratio;
    1410                colsol[i] = value;
    1411                objValue += value * cost[i];
    1412           }
    1413           for (i = 0; i < nrows; i++) {
    1414                double value = rowlower[i] - rowsol[i];
    1415                colsol[i+slackStart] = value;
    1416                objValue += value * cost[i+slackStart];
    1417           }
    1418           COIN_DETAIL_PRINT(printf("New objective after scaling %g\n", objValue));
    1419      }
     1334  double djTolerance = djTolerance_;
     1335  if (djTolerance > 0.0 && djTolerance < 1.0)
     1336    djTolerance = 1.0;
     1337  int iteration;
     1338  int i, n = 0;
     1339  double ratio = 1.0;
     1340  double objValue = 0.0;
     1341  if ((strategy_ & 128) != 0) {
     1342    fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
     1343  }
     1344  if ((mode & 16) != 0 && addAll < 3)
     1345    presolve = 1;
     1346  double *saveUpper = NULL;
     1347  double *saveLower = NULL;
     1348  double *saveRowUpper = NULL;
     1349  double *saveRowLower = NULL;
     1350  bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000);
     1351  if (addAll < 3) {
     1352    saveUpper = new double[ncols];
     1353    saveLower = new double[ncols];
     1354    CoinMemcpyN(upper, ncols, saveUpper);
     1355    CoinMemcpyN(lower, ncols, saveLower);
     1356    if (allowInfeasible) {
     1357      saveRowUpper = new double[nrows];
     1358      saveRowLower = new double[nrows];
     1359      CoinMemcpyN(rowupper, nrows, saveRowUpper);
     1360      CoinMemcpyN(rowlower, nrows, saveRowLower);
     1361      double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast< double >(model_->numberRows());
     1362      fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas);
     1363    }
     1364  }
     1365  if (slackStart >= 0) {
     1366    slackEnd = slackStart + nrows;
     1367    if (slackStart) {
     1368      ordStart = 0;
     1369      ordEnd = slackStart;
     1370    } else {
     1371      ordStart = nrows;
     1372      ordEnd = ncols;
     1373    }
     1374  } else {
     1375    slackEnd = slackStart;
     1376    ordStart = 0;
     1377    ordEnd = ncols;
     1378  }
     1379  /* get correct rowsol (without known slacks) */
     1380  memset(rowsol, 0, nrows * sizeof(double));
     1381  for (i = ordStart; i < ordEnd; i++) {
     1382    CoinBigIndex j;
     1383    double value = colsol[i];
     1384    if (value < lower[i] + fixTolerance) {
     1385      value = lower[i];
     1386      colsol[i] = value;
     1387    }
     1388    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
     1389      int irow = row[j];
     1390      rowsol[irow] += value * element[j];
     1391    }
     1392  }
     1393  if (slackStart >= 0) {
     1394    for (i = 0; i < nrows; i++) {
     1395      if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) {
     1396        ratio = rowlower[i] / rowsol[i];
     1397      }
     1398    }
     1399    for (i = 0; i < nrows; i++) {
     1400      rowsol[i] *= ratio;
     1401    }
     1402    for (i = ordStart; i < ordEnd; i++) {
     1403      double value = colsol[i] * ratio;
     1404      colsol[i] = value;
     1405      objValue += value * cost[i];
     1406    }
     1407    for (i = 0; i < nrows; i++) {
     1408      double value = rowlower[i] - rowsol[i];
     1409      colsol[i + slackStart] = value;
     1410      objValue += value * cost[i + slackStart];
     1411    }
     1412    COIN_DETAIL_PRINT(printf("New objective after scaling %g\n", objValue));
     1413  }
    14201414#if 0
    14211415     //maybe put back - but just get feasible ?
     
    14351429#define FEB_TRY
    14361430#ifdef FEB_TRY
    1437      int savePerturbation = model_->perturbation();
    1438      int saveOptions = model_->specialOptions();
    1439      model_->setSpecialOptions(saveOptions | 8192);
    1440      //if (savePerturbation_ == 50)
    1441      //   model_->setPerturbation(56);
    1442 #endif
    1443      model_->createStatus();
    1444      /* addAll
     1431  int savePerturbation = model_->perturbation();
     1432  int saveOptions = model_->specialOptions();
     1433  model_->setSpecialOptions(saveOptions | 8192);
     1434  //if (savePerturbation_ == 50)
     1435  //   model_->setPerturbation(56);
     1436#endif
     1437  model_->createStatus();
     1438  /* addAll
    14451439        0 - chosen,all used, all
    14461440        1 - chosen, all
     
    14481442        3 - do not do anything  - maybe basis
    14491443     */
    1450      for (i = ordStart; i < ordEnd; i++) {
    1451           if (addAll < 2) {
    1452                if (colsol[i] < lower[i] + fixTolerance) {
    1453                     upper[i] = lower[i];
    1454                     colsol[i] = lower[i];
    1455                } else if (colsol[i] > upper[i] - fixTolerance) {
    1456                     lower[i] = upper[i];
    1457                     colsol[i] = upper[i];
    1458                }
    1459           }
    1460           model_->setColumnStatus(i, ClpSimplex::superBasic);
    1461      }
    1462      if ((strategy_ & 16384) != 0) {
    1463           // put in basis
    1464           int * posSlack = whenUsed_ + ncols;
    1465           int * negSlack = posSlack + nrows;
    1466           int * nextSlack = negSlack + nrows;
    1467           /* Laci - try both ways - to see what works -
     1444  for (i = ordStart; i < ordEnd; i++) {
     1445    if (addAll < 2) {
     1446      if (colsol[i] < lower[i] + fixTolerance) {
     1447        upper[i] = lower[i];
     1448        colsol[i] = lower[i];
     1449      } else if (colsol[i] > upper[i] - fixTolerance) {
     1450        lower[i] = upper[i];
     1451        colsol[i] = upper[i];
     1452      }
     1453    }
     1454    model_->setColumnStatus(i, ClpSimplex::superBasic);
     1455  }
     1456  if ((strategy_ & 16384) != 0) {
     1457    // put in basis
     1458    int *posSlack = whenUsed_ + ncols;
     1459    int *negSlack = posSlack + nrows;
     1460    int *nextSlack = negSlack + nrows;
     1461    /* Laci - try both ways - to see what works -
    14681462             you can change second part as much as you want */
    1469 #ifndef LACI_TRY  // was #if 1
    1470           // Array for sorting out slack values
    1471           double * ratio = new double [ncols];
    1472           int * which = new int [ncols];
    1473           for (i = 0; i < nrows; i++) {
    1474                if (posSlack[i] >= 0 || negSlack[i] >= 0) {
    1475                     int iCol;
    1476                     int nPlus = 0;
    1477                     int nMinus = 0;
    1478                     bool possible = true;
    1479                     // Get sum
    1480                     double sum = 0.0;
    1481                     iCol = posSlack[i];
    1482                     while (iCol >= 0) {
    1483                          double value = element[columnStart[iCol]];
    1484                          sum += value * colsol[iCol];
    1485                          if (lower[iCol]) {
    1486                               possible = false;
    1487                               break;
    1488                          } else {
    1489                               nPlus++;
    1490                          }
    1491                          iCol = nextSlack[iCol];
    1492                     }
    1493                     iCol = negSlack[i];
    1494                     while (iCol >= 0) {
    1495                          double value = -element[columnStart[iCol]];
    1496                          sum -= value * colsol[iCol];
    1497                          if (lower[iCol]) {
    1498                               possible = false;
    1499                               break;
    1500                          } else {
    1501                               nMinus++;
    1502                          }
    1503                          iCol = nextSlack[iCol];
    1504                     }
    1505                     //printf("%d plus, %d minus",nPlus,nMinus);
    1506                     //printf("\n");
    1507                     if ((rowsol[i] - rowlower[i] < 1.0e-7 ||
    1508                               rowupper[i] - rowsol[i] < 1.0e-7) &&
    1509                               nPlus + nMinus < 2)
    1510                          possible = false;
    1511                     if (possible) {
    1512                          // Amount contributed by other varaibles
    1513                          sum = rowsol[i] - sum;
    1514                          double lo = rowlower[i];
    1515                          if (lo > -1.0e20)
    1516                               lo -= sum;
    1517                          double up = rowupper[i];
    1518                          if (up < 1.0e20)
    1519                               up -= sum;
    1520                          //printf("row bounds %g %g\n",lo,up);
    1521                          if (0) {
    1522                               double sum = 0.0;
    1523                               double x = 0.0;
    1524                               for (int k = 0; k < ncols; k++) {
    1525                                    CoinBigIndex j;
    1526                                    double value = colsol[k];
    1527                                    x += value * cost[k];
    1528                                    for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
    1529                                         int irow = row[j];
    1530                                         if (irow == i)
    1531                                              sum += element[j] * value;
    1532                                    }
    1533                               }
    1534                               printf("Before sum %g <= %g <= %g cost %.18g\n",
    1535                                      rowlower[i], sum, rowupper[i], x);
    1536                          }
    1537                          // set all to zero
    1538                          iCol = posSlack[i];
    1539                          while (iCol >= 0) {
    1540                               colsol[iCol] = 0.0;
    1541                               iCol = nextSlack[iCol];
    1542                          }
    1543                          iCol = negSlack[i];
    1544                          while (iCol >= 0) {
    1545                               colsol[iCol] = 0.0;
    1546                               iCol = nextSlack[iCol];
    1547                          }
    1548                          {
    1549                               int iCol;
    1550                               iCol = posSlack[i];
    1551                               while (iCol >= 0) {
    1552                                    //printf("col %d el %g sol %g bounds %g %g cost %g\n",
    1553                                    //     iCol,element[columnStart[iCol]],
    1554                                    //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
    1555                                    iCol = nextSlack[iCol];
    1556                               }
    1557                               iCol = negSlack[i];
    1558                               while (iCol >= 0) {
    1559                                    //printf("col %d el %g sol %g bounds %g %g cost %g\n",
    1560                                    //     iCol,element[columnStart[iCol]],
    1561                                    //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
    1562                                    iCol = nextSlack[iCol];
    1563                               }
    1564                          }
    1565                          //printf("now what?\n");
    1566                          int n = 0;
    1567                          bool basic = false;
    1568                          if (lo > 0.0) {
    1569                               // Add in positive
    1570                               iCol = posSlack[i];
    1571                               while (iCol >= 0) {
    1572                                    double value = element[columnStart[iCol]];
    1573                                    ratio[n] = cost[iCol] / value;
    1574                                    which[n++] = iCol;
    1575                                    iCol = nextSlack[iCol];
    1576                               }
    1577                               CoinSort_2(ratio, ratio + n, which);
    1578                               for (int i = 0; i < n; i++) {
    1579                                    iCol = which[i];
    1580                                    double value = element[columnStart[iCol]];
    1581                                    if (lo >= upper[iCol]*value) {
    1582                                         value *= upper[iCol];
    1583                                         sum += value;
    1584                                         lo -= value;
    1585                                         colsol[iCol] = upper[iCol];
    1586                                    } else {
    1587                                         value = lo / value;
    1588                                         sum += lo;
    1589                                         lo = 0.0;
    1590                                         colsol[iCol] = value;
    1591                                         model_->setColumnStatus(iCol, ClpSimplex::basic);
    1592                                         basic = true;
    1593                                    }
    1594                                    if (lo < 1.0e-7)
    1595                                         break;
    1596                               }
    1597                          } else if (up < 0.0) {
    1598                               // Use lo so coding is more similar
    1599                               lo = -up;
    1600                               // Add in negative
    1601                               iCol = negSlack[i];
    1602                               while (iCol >= 0) {
    1603                                    double value = -element[columnStart[iCol]];
    1604                                    ratio[n] = cost[iCol] / value;
    1605                                    which[n++] = iCol;
    1606                                    iCol = nextSlack[iCol];
    1607                               }
    1608                               CoinSort_2(ratio, ratio + n, which);
    1609                               for (int i = 0; i < n; i++) {
    1610                                    iCol = which[i];
    1611                                    double value = -element[columnStart[iCol]];
    1612                                    if (lo >= upper[iCol]*value) {
    1613                                         value *= upper[iCol];
    1614                                         sum += value;
    1615                                         lo -= value;
    1616                                         colsol[iCol] = upper[iCol];
    1617                                    } else {
    1618                                         value = lo / value;
    1619                                         sum += lo;
    1620                                         lo = 0.0;
    1621                                         colsol[iCol] = value;
    1622                                         model_->setColumnStatus(iCol, ClpSimplex::basic);
    1623                                         basic = true;
    1624                                    }
    1625                                    if (lo < 1.0e-7)
    1626                                         break;
    1627                               }
    1628                          }
    1629                          if (0) {
    1630                               double sum2 = 0.0;
    1631                               double x = 0.0;
    1632                               for (int k = 0; k < ncols; k++) {
    1633                                    CoinBigIndex j;
    1634                                    double value = colsol[k];
    1635                                    x += value * cost[k];
    1636                                    for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
    1637                                         int irow = row[j];
    1638                                         if (irow == i)
    1639                                              sum2 += element[j] * value;
    1640                                    }
    1641                               }
    1642                               printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n",
    1643                                      rowlower[i], sum2, rowupper[i], x, sum);
    1644                          }
    1645                          rowsol[i] = sum;
    1646                          if (basic) {
    1647                               if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
    1648                                    model_->setRowStatus(i, ClpSimplex::atLowerBound);
    1649                               else
    1650                                    model_->setRowStatus(i, ClpSimplex::atUpperBound);
    1651                          }
    1652                     } else {
    1653                          int n = 0;
    1654                          int iCol;
    1655                          iCol = posSlack[i];
    1656                          while (iCol >= 0) {
    1657                               if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
    1658                                         colsol[iCol] < upper[iCol] - 1.0e-8) {
    1659                                    model_->setColumnStatus(iCol, ClpSimplex::basic);
    1660                                    n++;
    1661                               }
    1662                               iCol = nextSlack[iCol];
    1663                          }
    1664                          iCol = negSlack[i];
    1665                          while (iCol >= 0) {
    1666                               if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
    1667                                         colsol[iCol] < upper[iCol] - 1.0e-8) {
    1668                                    model_->setColumnStatus(iCol, ClpSimplex::basic);
    1669                                    n++;
    1670                               }
    1671                               iCol = nextSlack[iCol];
    1672                          }
    1673                          if (n) {
    1674                               if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
    1675                                    model_->setRowStatus(i, ClpSimplex::atLowerBound);
    1676                               else
    1677                                    model_->setRowStatus(i, ClpSimplex::atUpperBound);
     1463#ifndef LACI_TRY // was #if 1
     1464    // Array for sorting out slack values
     1465    double *ratio = new double[ncols];
     1466    int *which = new int[ncols];
     1467    for (i = 0; i < nrows; i++) {
     1468      if (posSlack[i] >= 0 || negSlack[i] >= 0) {
     1469        int iCol;
     1470        int nPlus = 0;
     1471        int nMinus = 0;
     1472        bool possible = true;
     1473        // Get sum
     1474        double sum = 0.0;
     1475        iCol = posSlack[i];
     1476        while (iCol >= 0) {
     1477          double value = element[columnStart[iCol]];
     1478          sum += value * colsol[iCol];
     1479          if (lower[iCol]) {
     1480            possible = false;
     1481            break;
     1482          } else {
     1483            nPlus++;
     1484          }
     1485          iCol = nextSlack[iCol];
     1486        }
     1487        iCol = negSlack[i];
     1488        while (iCol >= 0) {
     1489          double value = -element[columnStart[iCol]];
     1490          sum -= value * colsol[iCol];
     1491          if (lower[iCol]) {
     1492            possible = false;
     1493            break;
     1494          } else {
     1495            nMinus++;
     1496          }
     1497          iCol = nextSlack[iCol];
     1498        }
     1499        //printf("%d plus, %d minus",nPlus,nMinus);
     1500        //printf("\n");
     1501        if ((rowsol[i] - rowlower[i] < 1.0e-7 || rowupper[i] - rowsol[i] < 1.0e-7) && nPlus + nMinus < 2)
     1502          possible = false;
     1503        if (possible) {
     1504          // Amount contributed by other varaibles
     1505          sum = rowsol[i] - sum;
     1506          double lo = rowlower[i];
     1507          if (lo > -1.0e20)
     1508            lo -= sum;
     1509          double up = rowupper[i];
     1510          if (up < 1.0e20)
     1511            up -= sum;
     1512          //printf("row bounds %g %g\n",lo,up);
     1513          if (0) {
     1514            double sum = 0.0;
     1515            double x = 0.0;
     1516            for (int k = 0; k < ncols; k++) {
     1517              CoinBigIndex j;
     1518              double value = colsol[k];
     1519              x += value * cost[k];
     1520              for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
     1521                int irow = row[j];
     1522                if (irow == i)
     1523                  sum += element[j] * value;
     1524              }
     1525            }
     1526            printf("Before sum %g <= %g <= %g cost %.18g\n",
     1527              rowlower[i], sum, rowupper[i], x);
     1528          }
     1529          // set all to zero
     1530          iCol = posSlack[i];
     1531          while (iCol >= 0) {
     1532            colsol[iCol] = 0.0;
     1533            iCol = nextSlack[iCol];
     1534          }
     1535          iCol = negSlack[i];
     1536          while (iCol >= 0) {
     1537            colsol[iCol] = 0.0;
     1538            iCol = nextSlack[iCol];
     1539          }
     1540          {
     1541            int iCol;
     1542            iCol = posSlack[i];
     1543            while (iCol >= 0) {
     1544              //printf("col %d el %g sol %g bounds %g %g cost %g\n",
     1545              //     iCol,element[columnStart[iCol]],
     1546              //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
     1547              iCol = nextSlack[iCol];
     1548            }
     1549            iCol = negSlack[i];
     1550            while (iCol >= 0) {
     1551              //printf("col %d el %g sol %g bounds %g %g cost %g\n",
     1552              //     iCol,element[columnStart[iCol]],
     1553              //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
     1554              iCol = nextSlack[iCol];
     1555            }
     1556          }
     1557          //printf("now what?\n");
     1558          int n = 0;
     1559          bool basic = false;
     1560          if (lo > 0.0) {
     1561            // Add in positive
     1562            iCol = posSlack[i];
     1563            while (iCol >= 0) {
     1564              double value = element[columnStart[iCol]];
     1565              ratio[n] = cost[iCol] / value;
     1566              which[n++] = iCol;
     1567              iCol = nextSlack[iCol];
     1568            }
     1569            CoinSort_2(ratio, ratio + n, which);
     1570            for (int i = 0; i < n; i++) {
     1571              iCol = which[i];
     1572              double value = element[columnStart[iCol]];
     1573              if (lo >= upper[iCol] * value) {
     1574                value *= upper[iCol];
     1575                sum += value;
     1576                lo -= value;
     1577                colsol[iCol] = upper[iCol];
     1578              } else {
     1579                value = lo / value;
     1580                sum += lo;
     1581                lo = 0.0;
     1582                colsol[iCol] = value;
     1583                model_->setColumnStatus(iCol, ClpSimplex::basic);
     1584                basic = true;
     1585              }
     1586              if (lo < 1.0e-7)
     1587                break;
     1588            }
     1589          } else if (up < 0.0) {
     1590            // Use lo so coding is more similar
     1591            lo = -up;
     1592            // Add in negative
     1593            iCol = negSlack[i];
     1594            while (iCol >= 0) {
     1595              double value = -element[columnStart[iCol]];
     1596              ratio[n] = cost[iCol] / value;
     1597              which[n++] = iCol;
     1598              iCol = nextSlack[iCol];
     1599            }
     1600            CoinSort_2(ratio, ratio + n, which);
     1601            for (int i = 0; i < n; i++) {
     1602              iCol = which[i];
     1603              double value = -element[columnStart[iCol]];
     1604              if (lo >= upper[iCol] * value) {
     1605                value *= upper[iCol];
     1606                sum += value;
     1607                lo -= value;
     1608                colsol[iCol] = upper[iCol];
     1609              } else {
     1610                value = lo / value;
     1611                sum += lo;
     1612                lo = 0.0;
     1613                colsol[iCol] = value;
     1614                model_->setColumnStatus(iCol, ClpSimplex::basic);
     1615                basic = true;
     1616              }
     1617              if (lo < 1.0e-7)
     1618                break;
     1619            }
     1620          }
     1621          if (0) {
     1622            double sum2 = 0.0;
     1623            double x = 0.0;
     1624            for (int k = 0; k < ncols; k++) {
     1625              CoinBigIndex j;
     1626              double value = colsol[k];
     1627              x += value * cost[k];
     1628              for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
     1629                int irow = row[j];
     1630                if (irow == i)
     1631                  sum2 += element[j] * value;
     1632              }
     1633            }
     1634            printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n",
     1635              rowlower[i], sum2, rowupper[i], x, sum);
     1636          }
     1637          rowsol[i] = sum;
     1638          if (basic) {
     1639            if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
     1640              model_->setRowStatus(i, ClpSimplex::atLowerBound);
     1641            else
     1642              model_->setRowStatus(i, ClpSimplex::atUpperBound);
     1643          }
     1644        } else {
     1645          int n = 0;
     1646          int iCol;
     1647          iCol = posSlack[i];
     1648          while (iCol >= 0) {
     1649            if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) {
     1650              model_->setColumnStatus(iCol, ClpSimplex::basic);
     1651              n++;
     1652            }
     1653            iCol = nextSlack[iCol];
     1654          }
     1655          iCol = negSlack[i];
     1656          while (iCol >= 0) {
     1657            if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) {
     1658              model_->setColumnStatus(iCol, ClpSimplex::basic);
     1659              n++;
     1660            }
     1661            iCol = nextSlack[iCol];
     1662          }
     1663          if (n) {
     1664            if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
     1665              model_->setRowStatus(i, ClpSimplex::atLowerBound);
     1666            else
     1667              model_->setRowStatus(i, ClpSimplex::atUpperBound);
    16781668#ifdef CLP_INVESTIGATE
    1679                               if (n > 1)
    1680                                    printf("%d basic on row %d!\n", n, i);
    1681 #endif
    1682                          }
    1683                     }
    1684                }
    1685           }
    1686           delete [] ratio;
    1687           delete [] which;
     1669            if (n > 1)
     1670              printf("%d basic on row %d!\n", n, i);
     1671#endif
     1672          }
     1673        }
     1674      }
     1675    }
     1676    delete[] ratio;
     1677    delete[] which;
    16881678#else
    1689           for (i = 0; i < nrows; i++) {
    1690                int n = 0;
    1691                int iCol;
    1692                iCol = posSlack[i];
    1693                while (iCol >= 0) {
    1694                     if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
    1695                               colsol[iCol] < upper[iCol] - 1.0e-8) {
    1696                          model_->setColumnStatus(iCol, ClpSimplex::basic);
    1697                          n++;
    1698                     }
    1699                     iCol = nextSlack[iCol];
    1700                }
    1701                iCol = negSlack[i];
    1702                while (iCol >= 0) {
    1703                     if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
    1704                               colsol[iCol] < upper[iCol] - 1.0e-8) {
    1705                          model_->setColumnStatus(iCol, ClpSimplex::basic);
    1706                          n++;
    1707                     }
    1708                     iCol = nextSlack[iCol];
    1709                }
    1710                if (n) {
    1711                     if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
    1712                          model_->setRowStatus(i, ClpSimplex::atLowerBound);
    1713                     else
    1714                          model_->setRowStatus(i, ClpSimplex::atUpperBound);
     1679    for (i = 0; i < nrows; i++) {
     1680      int n = 0;
     1681      int iCol;
     1682      iCol = posSlack[i];
     1683      while (iCol >= 0) {
     1684        if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) {
     1685          model_->setColumnStatus(iCol, ClpSimplex::basic);
     1686          n++;
     1687        }
     1688        iCol = nextSlack[iCol];
     1689      }
     1690      iCol = negSlack[i];
     1691      while (iCol >= 0) {
     1692        if (colsol[iCol] > lower[iCol] + 1.0e-8 && colsol[iCol] < upper[iCol] - 1.0e-8) {
     1693          model_->setColumnStatus(iCol, ClpSimplex::basic);
     1694          n++;
     1695        }
     1696        iCol = nextSlack[iCol];
     1697      }
     1698      if (n) {
     1699        if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
     1700          model_->setRowStatus(i, ClpSimplex::atLowerBound);
     1701        else
     1702          model_->setRowStatus(i, ClpSimplex::atUpperBound);
    17151703#ifdef CLP_INVESTIGATE
    1716                     if (n > 1)
    1717                          printf("%d basic on row %d!\n", n, i);
    1718 #endif
    1719                }
    1720           }
    1721 #endif
    1722      }
    1723      double maxmin;
    1724      if (model_->getObjSense() == -1.0) {
    1725           maxmin = -1.0;
    1726      } else {
    1727           maxmin = 1.0;
    1728      }
    1729      bool justValuesPass = majorIterations_ > 1000000;
    1730      if (slackStart >= 0) {
    1731           for (i = 0; i < nrows; i++) {
    1732                model_->setRowStatus(i, ClpSimplex::superBasic);
    1733           }
    1734           for (i = slackStart; i < slackEnd; i++) {
    1735                model_->setColumnStatus(i, ClpSimplex::basic);
    1736           }
    1737      } else {
    1738           /* still try and put singletons rather than artificials in basis */
    1739           for (i = 0; i < nrows; i++) {
    1740               model_->setRowStatus(i, ClpSimplex::basic);
    1741           }
    1742           int ninbas = 0;
    1743           for (i = 0; i < ncols; i++) {
    1744               if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
    1745                     CoinBigIndex j = columnStart[i];
    1746                     double value = element[j];
    1747                     int irow = row[j];
    1748                     double rlo = rowlower[irow];
    1749                     double rup = rowupper[irow];
    1750                     double clo = lower[i];
    1751                     double cup = upper[i];
    1752                     double csol = colsol[i];
    1753                     /* adjust towards feasibility */
    1754                     double move = 0.0;
    1755                     if (rowsol[irow] > rup) {
    1756                          move = (rup - rowsol[irow]) / value;
    1757                          if (value > 0.0) {
    1758                               /* reduce */
    1759                               if (csol + move < clo) move = CoinMin(0.0, clo - csol);
    1760                          } else {
    1761                               /* increase */
    1762                               if (csol + move > cup) move = CoinMax(0.0, cup - csol);
    1763                          }
    1764                     } else if (rowsol[irow] < rlo) {
    1765                          move = (rlo - rowsol[irow]) / value;
    1766                          if (value > 0.0) {
    1767                               /* increase */
    1768                               if (csol + move > cup) move = CoinMax(0.0, cup - csol);
    1769                          } else {
    1770                               /* reduce */
    1771                               if (csol + move < clo) move = CoinMin(0.0, clo - csol);
    1772                          }
    1773                     } else {
    1774                          /* move to improve objective */
    1775                          if (cost[i]*maxmin > 0.0) {
    1776                               if (value > 0.0) {
    1777                                    move = (rlo - rowsol[irow]) / value;
    1778                                    /* reduce */
    1779                                    if (csol + move < clo) move = CoinMin(0.0, clo - csol);
    1780                               } else {
    1781                                    move = (rup - rowsol[irow]) / value;
    1782                                    /* increase */
    1783                                    if (csol + move > cup) move = CoinMax(0.0, cup - csol);
    1784                               }
    1785                          } else if (cost[i]*maxmin < 0.0) {
    1786                               if (value > 0.0) {
    1787                                    move = (rup - rowsol[irow]) / value;
    1788                                    /* increase */
    1789                                    if (csol + move > cup) move = CoinMax(0.0, cup - csol);
    1790                               } else {
    1791                                    move = (rlo - rowsol[irow]) / value;
    1792                                    /* reduce */
    1793                                    if (csol + move < clo) move = CoinMin(0.0, clo - csol);
    1794                               }
    1795                          }
    1796                     }
    1797                     rowsol[irow] += move * value;
    1798                     colsol[i] += move;
    1799                     /* put in basis if row was artificial */
    1800                     if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) {
    1801                          model_->setRowStatus(irow, ClpSimplex::superBasic);
    1802                          model_->setColumnStatus(i, ClpSimplex::basic);
    1803                          ninbas++;
    1804                     }
    1805                }
    1806             }
    1807             /*printf("%d in basis\n",ninbas);*/
    1808      }
    1809      bool wantVector = false;
    1810      if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
    1811           // See if original wanted vector
    1812           ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
    1813           wantVector = clpMatrixO->wantsSpecialColumnCopy();
    1814      }
    1815      if ((strategy_&32768)!=0)
    1816        allowInfeasible=true;
    1817      if ((strategy_&65536)!=0)
    1818        justValuesPass=true;
    1819      //double * saveBounds=NULL;
    1820      if (addAll < 3) {
    1821           ClpPresolve pinfo;
    1822           if (presolve) {
    1823                double * rhs = new double[nrows];
    1824                double * saveBounds=new double[2*ncols];
    1825                char line[200];
    1826                memcpy(saveBounds,lower,ncols*sizeof(double));
    1827                memcpy(saveBounds+ncols,upper,ncols*sizeof(double));
    1828                if (allowInfeasible) {
    1829                     // fix up so will be feasible
    1830                     const double * dual = model_->dualRowSolution();
    1831                     for (i = 0; i < nrows; i++)
    1832                       rhs[i]=fabs(dual[i]);
    1833                     std::sort(rhs,rhs+nrows);
    1834                     int nSmall=nrows;
    1835                     int nMedium=nrows;
    1836                     double largest=rhs[nrows-1];
    1837                     double small=CoinMax(1.0e-4,1.0e-5*largest);
    1838                     small = CoinMin(small,1.0e-2);
    1839                     double medium=small*100.0;
    1840                     double * rowupper = model_->rowUpper();
    1841                     double * rowlower = model_->rowLower();
    1842                     // if tiny then drop row??
    1843                     for (i = 0; i < nrows; i++) {
    1844                       if (rhs[i]>=small) {
    1845                         nSmall=i-1;
    1846                         break;
    1847                       }
    1848                     }
    1849                     for (; i < nrows; i++) {
    1850                       if (rhs[i]>=medium) {
    1851                         nMedium=i-1;
    1852                         break;
    1853                       }
    1854                     }
    1855                     printf("%d < %g, %d < %g, %d <= %g\n",
    1856                            nSmall,small,nMedium-nSmall,medium,nrows-nMedium,largest);
    1857                     memset(rhs, 0, nrows * sizeof(double));
    1858                     int nFixed=0;
    1859                     for (i = 0; i < ncols; i++) {
    1860                       if (colsol[i]<lower[i]+1.0e-8) {
    1861                         upper[i]=lower[i];
    1862                         colsol[i]=lower[i];
    1863                         nFixed++;
    1864                       } else if (colsol[i]>upper[i]-1.0e-8) {
    1865                         lower[i]=upper[i];
    1866                         colsol[i]=lower[i];
    1867                         nFixed++;
    1868                       }
    1869                     }
    1870                     model_->clpMatrix()->times(1.0, colsol, rhs);
    1871                     saveRowUpper = CoinCopyOfArray(rowupper, nrows);
    1872                     saveRowLower = CoinCopyOfArray(rowlower, nrows);
    1873                     double sum = 0.0;
    1874                     for (i = 0; i < nrows; i++) {
    1875                          if (rhs[i] > rowupper[i]) {
    1876                               sum += rhs[i] - rowupper[i];
    1877                          }
    1878                          if (rhs[i] < rowlower[i]) {
    1879                               sum += rowlower[i] - rhs[i];
    1880                          }
    1881                     }
    1882                     double averageInfeasibility=sum/nrows;
    1883                     double check=CoinMin(1.0e-3,0.1*averageInfeasibility);
    1884                     int nFixedRows=0;
    1885                     int nFreed=0;
     1704        if (n > 1)
     1705          printf("%d basic on row %d!\n", n, i);
     1706#endif
     1707      }
     1708    }
     1709#endif
     1710  }
     1711  double maxmin;
     1712  if (model_->getObjSense() == -1.0) {
     1713    maxmin = -1.0;
     1714  } else {
     1715    maxmin = 1.0;
     1716  }
     1717  bool justValuesPass = majorIterations_ > 1000000;
     1718  if (slackStart >= 0) {
     1719    for (i = 0; i < nrows; i++) {
     1720      model_->setRowStatus(i, ClpSimplex::superBasic);
     1721    }
     1722    for (i = slackStart; i < slackEnd; i++) {
     1723      model_->setColumnStatus(i, ClpSimplex::basic);
     1724    }
     1725  } else {
     1726    /* still try and put singletons rather than artificials in basis */
     1727    for (i = 0; i < nrows; i++) {
     1728      model_->setRowStatus(i, ClpSimplex::basic);
     1729    }
     1730    int ninbas = 0;
     1731    for (i = 0; i < ncols; i++) {
     1732      if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
     1733        CoinBigIndex j = columnStart[i];
     1734        double value = element[j];
     1735        int irow = row[j];
     1736        double rlo = rowlower[irow];
     1737        double rup = rowupper[irow];
     1738        double clo = lower[i];
     1739        double cup = upper[i];
     1740        double csol = colsol[i];
     1741        /* adjust towards feasibility */
     1742        double move = 0.0;
     1743        if (rowsol[irow] > rup) {
     1744          move = (rup - rowsol[irow]) / value;
     1745          if (value > 0.0) {
     1746            /* reduce */
     1747            if (csol + move < clo)
     1748              move = CoinMin(0.0, clo - csol);
     1749          } else {
     1750            /* increase */
     1751            if (csol + move > cup)
     1752              move = CoinMax(0.0, cup - csol);
     1753          }
     1754        } else if (rowsol[irow] < rlo) {
     1755          move = (rlo - rowsol[irow]) / value;
     1756          if (value > 0.0) {
     1757            /* increase */
     1758            if (csol + move > cup)
     1759              move = CoinMax(0.0, cup - csol);
     1760          } else {
     1761            /* reduce */
     1762            if (csol + move < clo)
     1763              move = CoinMin(0.0, clo - csol);
     1764          }
     1765        } else {
     1766          /* move to improve objective */
     1767          if (cost[i] * maxmin > 0.0) {
     1768            if (value > 0.0) {
     1769              move = (rlo - rowsol[irow]) / value;
     1770              /* reduce */
     1771              if (csol + move < clo)
     1772                move = CoinMin(0.0, clo - csol);
     1773            } else {
     1774              move = (rup - rowsol[irow]) / value;
     1775              /* increase */
     1776              if (csol + move > cup)
     1777                move = CoinMax(0.0, cup - csol);
     1778            }
     1779          } else if (cost[i] * maxmin < 0.0) {
     1780            if (value > 0.0) {
     1781              move = (rup - rowsol[irow]) / value;
     1782              /* increase */
     1783              if (csol + move > cup)
     1784                move = CoinMax(0.0, cup - csol);
     1785            } else {
     1786              move = (rlo - rowsol[irow]) / value;
     1787              /* reduce */
     1788              if (csol + move < clo)
     1789                move = CoinMin(0.0, clo - csol);
     1790            }
     1791          }
     1792        }
     1793        rowsol[irow] += move * value;
     1794        colsol[i] += move;
     1795        /* put in basis if row was artificial */
     1796        if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) {
     1797          model_->setRowStatus(irow, ClpSimplex::superBasic);
     1798          model_->setColumnStatus(i, ClpSimplex::basic);
     1799          ninbas++;
     1800        }
     1801      }
     1802    }
     1803    /*printf("%d in basis\n",ninbas);*/
     1804  }
     1805  bool wantVector = false;
     1806  if (dynamic_cast< ClpPackedMatrix * >(model_->clpMatrix())) {
     1807    // See if original wanted vector
     1808    ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(model_->clpMatrix());
     1809    wantVector = clpMatrixO->wantsSpecialColumnCopy();
     1810  }
     1811  if ((strategy_ & 32768) != 0)
     1812    allowInfeasible = true;
     1813  if ((strategy_ & 65536) != 0)
     1814    justValuesPass = true;
     1815  //double * saveBounds=NULL;
     1816  if (addAll < 3) {
     1817    ClpPresolve pinfo;
     1818    if (presolve) {
     1819      double *rhs = new double[nrows];
     1820      double *saveBounds = new double[2 * ncols];
     1821      char line[200];
     1822      memcpy(saveBounds, lower, ncols * sizeof(double));
     1823      memcpy(saveBounds + ncols, upper, ncols * sizeof(double));
     1824      if (allowInfeasible) {
     1825        // fix up so will be feasible
     1826        const double *dual = model_->dualRowSolution();
     1827        for (i = 0; i < nrows; i++)
     1828          rhs[i] = fabs(dual[i]);
     1829        std::sort(rhs, rhs + nrows);
     1830        int nSmall = nrows;
     1831        int nMedium = nrows;
     1832        double largest = rhs[nrows - 1];
     1833        double small = CoinMax(1.0e-4, 1.0e-5 * largest);
     1834        small = CoinMin(small, 1.0e-2);
     1835        double medium = small * 100.0;
     1836        double *rowupper = model_->rowUpper();
     1837        double *rowlower = model_->rowLower();
     1838        // if tiny then drop row??
     1839        for (i = 0; i < nrows; i++) {
     1840          if (rhs[i] >= small) {
     1841            nSmall = i - 1;
     1842            break;
     1843          }
     1844        }
     1845        for (; i < nrows; i++) {
     1846          if (rhs[i] >= medium) {
     1847            nMedium = i - 1;
     1848            break;
     1849          }
     1850        }
     1851        printf("%d < %g, %d < %g, %d <= %g\n",
     1852          nSmall, small, nMedium - nSmall, medium, nrows - nMedium, largest);
     1853        memset(rhs, 0, nrows * sizeof(double));
     1854        int nFixed = 0;
     1855        for (i = 0; i < ncols; i++) {
     1856          if (colsol[i] < lower[i] + 1.0e-8) {
     1857            upper[i] = lower[i];
     1858            colsol[i] = lower[i];
     1859            nFixed++;
     1860          } else if (colsol[i] > upper[i] - 1.0e-8) {
     1861            lower[i] = upper[i];
     1862            colsol[i] = lower[i];
     1863            nFixed++;
     1864          }
     1865        }
     1866        model_->clpMatrix()->times(1.0, colsol, rhs);
     1867        saveRowUpper = CoinCopyOfArray(rowupper, nrows);
     1868        saveRowLower = CoinCopyOfArray(rowlower, nrows);
     1869        double sum = 0.0;
     1870        for (i = 0; i < nrows; i++) {
     1871          if (rhs[i] > rowupper[i]) {
     1872            sum += rhs[i] - rowupper[i];
     1873          }
     1874          if (rhs[i] < rowlower[i]) {
     1875            sum += rowlower[i] - rhs[i];
     1876          }
     1877        }
     1878        double averageInfeasibility = sum / nrows;
     1879        double check = CoinMin(1.0e-3, 0.1 * averageInfeasibility);
     1880        int nFixedRows = 0;
     1881        int nFreed = 0;
    18861882#define MESS_UP 0
    1887                     for (i = 0; i < nrows; i++) {
    1888                          if (rowupper[i]>rowlower[i]+check) {
    1889                            // look at distance and sign of dual
    1890                            if (dual[i]<-medium&&rowupper[i]-rhs[i]<check) {
    1891                              rowupper[i]=rhs[i];
    1892                              rowlower[i]=rowupper[i];
    1893                              nFixedRows++;
    1894                            } else if (dual[i]>medium&&rhs[i]-rowlower[i]<check) {
    1895                              rowlower[i]=rhs[i];
    1896                              rowupper[i]=rowlower[i];
    1897                              nFixedRows++;
    1898                            } else if (fabs(dual[i])<small&&
    1899                                       rhs[i]-rowlower[i]>check&&
    1900                                       rowupper[i]-rhs[i]>check) {
    1901                              nFreed++;
    1902 #if MESS_UP ==1 || MESS_UP ==2
    1903                              rowupper[i]=COIN_DBL_MAX;
    1904                              rowlower[i]=-COIN_DBL_MAX;
    1905 #endif
    1906                            }
    1907                          }
    1908                          if (rhs[i] > rowupper[i]) {
    1909                               rowupper[i] = rhs[i];
    1910                               // maybe make equality
    1911 #if MESS_UP ==2 || MESS_UP ==3
    1912                               rowlower[i] = rhs[i];
    1913 #endif
    1914                          }
    1915                          if (rhs[i] < rowlower[i]) {
    1916                               rowlower[i] = rhs[i];
    1917                               // maybe make equality
    1918 #if MESS_UP ==2 || MESS_UP ==3
    1919                               rowupper[i] = rhs[i];
    1920 #endif
    1921                          }
    1922                     }
    1923                     sprintf(line,"sum of infeasibilities %g - %d fixed rows, %d fixed columns - might free %d rows",
    1924                            sum,nFixedRows,nFixed,nFreed);
    1925                } else {
    1926                     memset(rhs, 0, nrows * sizeof(double));
    1927                     int nFixed=0;
    1928                     for (i = 0; i < ncols; i++) {
    1929                       if (colsol[i]<lower[i]+1.0e-8) {
    1930                         upper[i]=lower[i];
    1931                         colsol[i]=lower[i];
    1932                         nFixed++;
    1933                       } else if (colsol[i]>upper[i]-1.0e-8) {
    1934                         lower[i]=upper[i];
    1935                         colsol[i]=lower[i];
    1936                         nFixed++;
    1937                       }
    1938                     }
    1939                     model_->clpMatrix()->times(1.0, colsol, rhs);
    1940                     double sum = 0.0;
    1941                     for (i = 0; i < nrows; i++) {
    1942                          if (rhs[i] > rowupper[i]) {
    1943                               sum += rhs[i] - rowupper[i];
    1944                          }
    1945                          if (rhs[i] < rowlower[i]) {
    1946                               sum += rowlower[i] - rhs[i];
    1947                          }
    1948                     }
     1883        for (i = 0; i < nrows; i++) {
     1884          if (rowupper[i] > rowlower[i] + check) {
     1885            // look at distance and sign of dual
     1886            if (dual[i] < -medium && rowupper[i] - rhs[i] < check) {
     1887              rowupper[i] = rhs[i];
     1888              rowlower[i] = rowupper[i];
     1889              nFixedRows++;
     1890            } else if (dual[i] > medium && rhs[i] - rowlower[i] < check) {
     1891              rowlower[i] = rhs[i];
     1892              rowupper[i] = rowlower[i];
     1893              nFixedRows++;
     1894            } else if (fabs(dual[i]) < small && rhs[i] - rowlower[i] > check && rowupper[i] - rhs[i] > check) {
     1895              nFreed++;
     1896#if MESS_UP == 1 || MESS_UP == 2
     1897              rowupper[i] = COIN_DBL_MAX;
     1898              rowlower[i] = -COIN_DBL_MAX;
     1899#endif
     1900            }
     1901          }
     1902          if (rhs[i] > rowupper[i]) {
     1903            rowupper[i] = rhs[i];
     1904            // maybe make equality
     1905#if MESS_UP == 2 || MESS_UP == 3
     1906            rowlower[i] = rhs[i];
     1907#endif
     1908          }
     1909          if (rhs[i] < rowlower[i]) {
     1910            rowlower[i] = rhs[i];
     1911            // maybe make equality
     1912#if MESS_UP == 2 || MESS_UP == 3
     1913            rowupper[i] = rhs[i];
     1914#endif
     1915          }
     1916        }
     1917        sprintf(line, "sum of infeasibilities %g - %d fixed rows, %d fixed columns - might free %d rows",
     1918          sum, nFixedRows, nFixed, nFreed);
     1919      } else {
     1920        memset(rhs, 0, nrows * sizeof(double));
     1921        int nFixed = 0;
     1922        for (i = 0; i < ncols; i++) {
     1923          if (colsol[i] < lower[i] + 1.0e-8) {
     1924            upper[i] = lower[i];
     1925            colsol[i] = lower[i];
     1926            nFixed++;
     1927          } else if (colsol[i] > upper[i] - 1.0e-8) {
     1928            lower[i] = upper[i];
     1929            colsol[i] = lower[i];
     1930            nFixed++;
     1931          }
     1932        }
     1933        model_->clpMatrix()->times(1.0, colsol, rhs);
     1934        double sum = 0.0;
     1935        for (i = 0; i < nrows; i++) {
     1936          if (rhs[i] > rowupper[i]) {
     1937            sum += rhs[i] - rowupper[i];
     1938          }
     1939          if (rhs[i] < rowlower[i]) {
     1940            sum += rowlower[i] - rhs[i];
     1941          }
     1942        }
    19491943
    1950                     double averageInfeasibility=sum/nrows;
    1951                     sprintf(line,"sum of infeasibilities %g - average %g, %d fixed columns",
    1952                            sum,averageInfeasibility,nFixed);
    1953                }
    1954                const CoinMessages * messages = model_->messagesPointer();
    1955                model_->messageHandler()->message(CLP_GENERAL, *messages)
    1956                 << line
    1957                 << CoinMessageEol;
    1958                delete [] rhs;
    1959                saveModel = model_;
    1960                pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
    1961                model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
    1962                if (saveBounds) {
    1963                  memcpy(saveModel->columnLower(),saveBounds,ncols*sizeof(double));
    1964                  memcpy(saveModel->columnUpper(),saveBounds+ncols,ncols*sizeof(double));
    1965                  delete [] saveBounds;
    1966                }
    1967                if (model_&&(strategy_&262144)!=0) {
    1968                 int nrows = model_->getNumRows();
    1969                 int ncols = model_->getNumCols();
    1970                  double * lower = model_->columnLower();
    1971                  double * upper = model_->columnUpper();
    1972                  const double * rowlower = model_->getRowLower();
    1973                  const double * rowupper = model_->getRowUpper();
    1974                  double * rowsol = model_->primalRowSolution();
    1975                  double * colsol = model_->primalColumnSolution();;
    1976                  int ninbas = 0;
    1977                  int * which = new int[2*ncols+nrows];
    1978                  double * dj = model_->dualColumnSolution();
    1979                  for (int i=0;i<ncols;i++) {
    1980                    dj[i]=-CoinMin(upper[i]-colsol[i],colsol[i]-lower[i]);
    1981                    which[i]=i;
    1982                  }
    1983                  CoinSort_2(dj,dj+ncols,which);
    1984                  ninbas=CoinMin(ncols,nrows);
    1985                  int * columnIsBasic=which+ncols;
    1986                  int * rowIsBasic=columnIsBasic+ncols;
    1987                  for (int i=0;i<nrows+ncols;i++)
    1988                    columnIsBasic[i]=-1;
    1989                  for (int i=0;i<ninbas;i++) {
    1990                    int iColumn=which[i];
    1991                    columnIsBasic[iColumn]=i;
    1992                  }
    1993                  // factorize
    1994                  CoinFactorization factor;
    1995                  factor.pivotTolerance(0.1);
    1996                  factor.setDenseThreshold(0);
    1997                  int status=-1;
    1998                  // If initial is too dense - then all-slack may be better
    1999                  double areaFactor=1.0; // was 4.0
    2000                  const CoinPackedMatrix * matrix = model_->matrix();
    2001                  while (status) {
    2002                    status=factor.factorize(*matrix,rowIsBasic,columnIsBasic,areaFactor);
    2003                    if (status==-99) {
    2004                      // put all slacks in
    2005                      for (int i=0;i<nrows;i++)
    2006                        rowIsBasic[i]=i;
    2007                      for (int i=0;i<ncols;i++)
    2008                       columnIsBasic[i]=-1;
    2009                      break;
    2010                    } else if (status==-1) {
    2011                      factor.pivotTolerance(0.99);
    2012                      // put all slacks in
    2013                      for (int i=0;i<nrows;i++)
    2014                        rowIsBasic[i]=i;
    2015                      for (int i=0;i<ncols;i++) {
    2016                        int iRow=columnIsBasic[i];
    2017                        if (iRow>=0)
    2018                          rowIsBasic[iRow]=-1; // out
    2019                      }
    2020                    }
    2021                 }
    2022                  for (int i = 0; i < nrows; i++) {
    2023                    if (rowIsBasic[i]>=0) {
    2024                      model_->setRowStatus(i, ClpSimplex::basic);
    2025                    } else if (rowlower[i]==rowupper[i]) {
    2026                      model_->setRowStatus(i, ClpSimplex::isFixed);
    2027                    } else if (rowsol[i]-rowlower[i]<rowupper[i]-rowsol[i]) {
    2028                      model_->setRowStatus(i, ClpSimplex::atLowerBound);
    2029                    } else {
    2030                      model_->setRowStatus(i, ClpSimplex::atUpperBound);
    2031                    }
    2032                 }
    2033                  for (int i=0;i<ncols;i++) {
    2034                    if (colsol[i]>upper[i]-1.0e-7||
    2035                        colsol[i]<lower[i]+1.0e-7) {
    2036                      model_->setColumnStatus(i,ClpSimplex::isFixed);
    2037                    } else if (columnIsBasic[i]>=0) {
    2038                      model_->setColumnStatus(i,ClpSimplex::basic);
    2039                    } else {
    2040                      model_->setColumnStatus(i,ClpSimplex::superBasic);
    2041                    }
    2042                 }
    2043                  delete [] which;
    2044                }
    2045           }
    2046           if (model_) {
    2047                // See if we want to go all way
    2048                int oldSize=2*saveModel->numberRows()+saveModel->numberColumns();
    2049                int newSize=2*model_->numberRows()+model_->numberColumns();
    2050                if (oldSize*2>newSize*3)
    2051                  justValuesPass=false;
    2052                if (!wantVector) {
    2053                     //#define TWO_GOES
     1944        double averageInfeasibility = sum / nrows;
     1945        sprintf(line, "sum of infeasibilities %g - average %g, %d fixed columns",
     1946          sum, averageInfeasibility, nFixed);
     1947      }
     1948      const CoinMessages *messages = model_->messagesPointer();
     1949      model_->messageHandler()->message(CLP_GENERAL, *messages)
     1950        << line
     1951        << CoinMessageEol;
     1952      delete[] rhs;
     1953      saveModel = model_;
     1954      pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
     1955      model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
     1956      if (saveBounds) {
     1957        memcpy(saveModel->columnLower(), saveBounds, ncols * sizeof(double));
     1958        memcpy(saveModel->columnUpper(), saveBounds + ncols, ncols * sizeof(double));
     1959        delete[] saveBounds;
     1960      }
     1961      if (model_ && (strategy_ & 262144) != 0) {
     1962        int nrows = model_->getNumRows();
     1963        int ncols = model_->getNumCols();
     1964        double *lower = model_->columnLower();
     1965        double *upper = model_->columnUpper();
     1966        const double *rowlower = model_->getRowLower();
     1967        const double *rowupper = model_->getRowUpper();
     1968        double *rowsol = model_->primalRowSolution();
     1969        double *colsol = model_->primalColumnSolution();
     1970        ;
     1971        int ninbas = 0;
     1972        int *which = new int[2 * ncols + nrows];
     1973        double *dj = model_->dualColumnSolution();
     1974        for (int i = 0; i < ncols; i++) {
     1975          dj[i] = -CoinMin(upper[i] - colsol[i], colsol[i] - lower[i]);
     1976          which[i] = i;
     1977        }
     1978        CoinSort_2(dj, dj + ncols, which);
     1979        ninbas = CoinMin(ncols, nrows);
     1980        int *columnIsBasic = which + ncols;
     1981        int *rowIsBasic = columnIsBasic + ncols;
     1982        for (int i = 0; i < nrows + ncols; i++)
     1983          columnIsBasic[i] = -1;
     1984        for (int i = 0; i < ninbas; i++) {
     1985          int iColumn = which[i];
     1986          columnIsBasic[iColumn] = i;
     1987        }
     1988        // factorize
     1989        CoinFactorization factor;
     1990        factor.pivotTolerance(0.1);
     1991        factor.setDenseThreshold(0);
     1992        int status = -1;
     1993        // If initial is too dense - then all-slack may be better
     1994        double areaFactor = 1.0; // was 4.0
     1995        const CoinPackedMatrix *matrix = model_->matrix();
     1996        while (status) {
     1997          status = factor.factorize(*matrix, rowIsBasic, columnIsBasic, areaFactor);
     1998          if (status == -99) {
     1999            // put all slacks in
     2000            for (int i = 0; i < nrows; i++)
     2001              rowIsBasic[i] = i;
     2002            for (int i = 0; i < ncols; i++)
     2003              columnIsBasic[i] = -1;
     2004            break;
     2005          } else if (status == -1) {
     2006            factor.pivotTolerance(0.99);
     2007            // put all slacks in
     2008            for (int i = 0; i < nrows; i++)
     2009              rowIsBasic[i] = i;
     2010            for (int i = 0; i < ncols; i++) {
     2011              int iRow = columnIsBasic[i];
     2012              if (iRow >= 0)
     2013                rowIsBasic[iRow] = -1; // out
     2014            }
     2015          }
     2016        }
     2017        for (int i = 0; i < nrows; i++) {
     2018          if (rowIsBasic[i] >= 0) {
     2019            model_->setRowStatus(i, ClpSimplex::basic);
     2020          } else if (rowlower[i] == rowupper[i]) {
     2021            model_->setRowStatus(i, ClpSimplex::isFixed);
     2022          } else if (rowsol[i] - rowlower[i] < rowupper[i] - rowsol[i]) {
     2023            model_->setRowStatus(i, ClpSimplex::atLowerBound);
     2024          } else {
     2025            model_->setRowStatus(i, ClpSimplex::atUpperBound);
     2026          }
     2027        }
     2028        for (int i = 0; i < ncols; i++) {
     2029          if (colsol[i] > upper[i] - 1.0e-7 || colsol[i] < lower[i] + 1.0e-7) {
     2030            model_->setColumnStatus(i, ClpSimplex::isFixed);
     2031          } else if (columnIsBasic[i] >= 0) {
     2032            model_->setColumnStatus(i, ClpSimplex::basic);
     2033          } else {
     2034            model_->setColumnStatus(i, ClpSimplex::superBasic);
     2035          }
     2036        }
     2037        delete[] which;
     2038      }
     2039    }
     2040    if (model_) {
     2041      // See if we want to go all way
     2042      int oldSize = 2 * saveModel->numberRows() + saveModel->numberColumns();
     2043      int newSize = 2 * model_->numberRows() + model_->numberColumns();
     2044      if (oldSize * 2 > newSize * 3)
     2045        justValuesPass = false;
     2046      if (!wantVector) {
     2047        //#define TWO_GOES
    20542048#ifdef ABC_INHERIT
    20552049#ifndef TWO_GOES
    2056                     model_->dealWithAbc(1,justValuesPass ? 3 : 1);
     2050        model_->dealWithAbc(1, justValuesPass ? 3 : 1);
    20572051#else
    2058                     model_->dealWithAbc(1,1 + 11);
     2052        model_->dealWithAbc(1, 1 + 11);
    20592053#endif
    20602054#else
    20612055#ifndef TWO_GOES
    2062                     model_->primal(justValuesPass ? 2 : 1);
     2056        model_->primal(justValuesPass ? 2 : 1);
    20632057#else
    2064                     model_->primal(1 + 11);
    2065 #endif
    2066 #endif
    2067                } else {
    2068                     ClpMatrixBase * matrix = model_->clpMatrix();
    2069                     ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    2070                     assert (clpMatrix);
    2071                     clpMatrix->makeSpecialColumnCopy();
     2058        model_->primal(1 + 11);
     2059#endif
     2060#endif
     2061      } else {
     2062        ClpMatrixBase *matrix = model_->clpMatrix();
     2063        ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     2064        assert(clpMatrix);
     2065        clpMatrix->makeSpecialColumnCopy();
    20722066#ifdef ABC_INHERIT
    2073                     model_->dealWithAbc(1,1);
     2067        model_->dealWithAbc(1, 1);
    20742068#else
    2075                     model_->primal(1);
    2076 #endif
    2077                     clpMatrix->releaseSpecialColumnCopy();
    2078                }
    2079                if (presolve) {
    2080                 if (!justValuesPass)
    2081                    model_->primal(1);
    2082                     pinfo.postsolve(true);
    2083                     delete model_;
    2084                     model_ = saveModel;
    2085                     saveModel = NULL;
    2086                }
    2087           } else {
    2088                // not feasible
    2089                addAll = 1;
    2090                presolve = 0;
    2091                model_ = saveModel;
    2092                saveModel = NULL;
    2093                if (justValuesPass)
     2069        model_->primal(1);
     2070#endif
     2071        clpMatrix->releaseSpecialColumnCopy();
     2072      }
     2073      if (presolve) {
     2074        if (!justValuesPass)
     2075          model_->primal(1);
     2076        pinfo.postsolve(true);
     2077        delete model_;
     2078        model_ = saveModel;
     2079        saveModel = NULL;
     2080      }
     2081    } else {
     2082      // not feasible
     2083      addAll = 1;
     2084      presolve = 0;
     2085      model_ = saveModel;
     2086      saveModel = NULL;
     2087      if (justValuesPass)
    20942088#ifdef ABC_INHERIT
    2095                     model_->dealWithAbc(1,3);
     2089        model_->dealWithAbc(1, 3);
    20962090#else
    2097                     model_->primal(2);
    2098 #endif
    2099           }
    2100           if (allowInfeasible) {
    2101                CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper());
    2102                CoinMemcpyN(saveRowLower, nrows, model_->rowLower());
    2103                delete [] saveRowUpper;
    2104                delete [] saveRowLower;
    2105                saveRowUpper = NULL;
    2106                saveRowLower = NULL;
    2107           }
    2108           if (addAll < 2) {
    2109                n = 0;
    2110                if (!addAll ) {
    2111                     /* could do scans to get a good number */
    2112                     iteration = 1;
    2113                     for (i = ordStart; i < ordEnd; i++) {
    2114                          if (whenUsed[i] >= iteration) {
    2115                               if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
    2116                                    n++;
    2117                                    upper[i] = saveUpper[i];
    2118                                    lower[i] = saveLower[i];
    2119                               }
    2120                          }
    2121                     }
    2122                } else {
    2123                     for (i = ordStart; i < ordEnd; i++) {
    2124                          if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
    2125                               n++;
    2126                               upper[i] = saveUpper[i];
    2127                               lower[i] = saveLower[i];
    2128                          }
    2129                     }
    2130                     delete [] saveUpper;
    2131                     delete [] saveLower;
    2132                     saveUpper = NULL;
    2133                     saveLower = NULL;
    2134                }
     2091        model_->primal(2);
     2092#endif
     2093    }
     2094    if (allowInfeasible) {
     2095      CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper());
     2096      CoinMemcpyN(saveRowLower, nrows, model_->rowLower());
     2097      delete[] saveRowUpper;
     2098      delete[] saveRowLower;
     2099      saveRowUpper = NULL;
     2100      saveRowLower = NULL;
     2101    }
     2102    if (addAll < 2) {
     2103      n = 0;
     2104      if (!addAll) {
     2105        /* could do scans to get a good number */
     2106        iteration = 1;
     2107        for (i = ordStart; i < ordEnd; i++) {
     2108          if (whenUsed[i] >= iteration) {
     2109            if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
     2110              n++;
     2111              upper[i] = saveUpper[i];
     2112              lower[i] = saveLower[i];
     2113            }
     2114          }
     2115        }
     2116      } else {
     2117        for (i = ordStart; i < ordEnd; i++) {
     2118          if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
     2119            n++;
     2120            upper[i] = saveUpper[i];
     2121            lower[i] = saveLower[i];
     2122          }
     2123        }
     2124        delete[] saveUpper;
     2125        delete[] saveLower;
     2126        saveUpper = NULL;
     2127        saveLower = NULL;
     2128      }
    21352129#ifdef COIN_DEVELOP
    2136                printf("Time so far %g, %d now added from previous iterations\n",
    2137                       CoinCpuTime() - startTime, n);
    2138 #endif
    2139                if (justValuesPass)
    2140                     return;
    2141                if (addAll)
    2142                     presolve = 0;
    2143                if (presolve) {
    2144                     saveModel = model_;
    2145                     model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
    2146                } else {
    2147                     presolve = 0;
    2148                }
    2149                if (!wantVector) {
     2130      printf("Time so far %g, %d now added from previous iterations\n",
     2131        CoinCpuTime() - startTime, n);
     2132#endif
     2133      if (justValuesPass)
     2134        return;
     2135      if (addAll)
     2136        presolve = 0;
     2137      if (presolve) {
     2138        saveModel = model_;
     2139        model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
     2140      } else {
     2141        presolve = 0;
     2142      }
     2143      if (!wantVector) {
    21502144#ifdef ABC_INHERIT
    2151                     model_->dealWithAbc(1,1);
     2145        model_->dealWithAbc(1, 1);
    21522146#else
    2153                     model_->primal(1);
    2154 #endif
    2155                } else {
    2156                     ClpMatrixBase * matrix = model_->clpMatrix();
    2157                     ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    2158                     assert (clpMatrix);
    2159                     clpMatrix->makeSpecialColumnCopy();
     2147        model_->primal(1);
     2148#endif
     2149      } else {
     2150        ClpMatrixBase *matrix = model_->clpMatrix();
     2151        ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     2152        assert(clpMatrix);
     2153        clpMatrix->makeSpecialColumnCopy();
    21602154#ifdef ABC_INHERIT
    2161                     model_->dealWithAbc(1,1);
     2155        model_->dealWithAbc(1, 1);
    21622156#else
    2163                     model_->primal(1);
    2164 #endif
    2165                     clpMatrix->releaseSpecialColumnCopy();
    2166                }
    2167                if (presolve) {
    2168                     pinfo.postsolve(true);
    2169                     delete model_;
    2170                     model_ = saveModel;
    2171                     saveModel = NULL;
    2172                }
    2173                if (!addAll) {
    2174                     n = 0;
    2175                     for (i = ordStart; i < ordEnd; i++) {
    2176                          if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
    2177                               n++;
    2178                               upper[i] = saveUpper[i];
    2179                               lower[i] = saveLower[i];
    2180                          }
    2181                     }
    2182                     delete [] saveUpper;
    2183                     delete [] saveLower;
    2184                     saveUpper = NULL;
    2185                     saveLower = NULL;
     2157        model_->primal(1);
     2158#endif
     2159        clpMatrix->releaseSpecialColumnCopy();
     2160      }
     2161      if (presolve) {
     2162        pinfo.postsolve(true);
     2163        delete model_;
     2164        model_ = saveModel;
     2165        saveModel = NULL;
     2166      }
     2167      if (!addAll) {
     2168        n = 0;
     2169        for (i = ordStart; i < ordEnd; i++) {
     2170          if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
     2171            n++;
     2172            upper[i] = saveUpper[i];
     2173            lower[i] = saveLower[i];
     2174          }
     2175        }
     2176        delete[] saveUpper;
     2177        delete[] saveLower;
     2178        saveUpper = NULL;
     2179        saveLower = NULL;
    21862180#ifdef COIN_DEVELOP
    2187                     printf("Time so far %g, %d now added from previous iterations\n",
    2188                            CoinCpuTime() - startTime, n);
    2189 #endif
    2190                }
    2191                if (presolve) {
    2192                     saveModel = model_;
    2193                     model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
    2194                } else {
    2195                     presolve = 0;
    2196                }
    2197                if (!wantVector) {
     2181        printf("Time so far %g, %d now added from previous iterations\n",
     2182          CoinCpuTime() - startTime, n);
     2183#endif
     2184      }
     2185      if (presolve) {
     2186        saveModel = model_;
     2187        model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
     2188      } else {
     2189        presolve = 0;
     2190      }
     2191      if (!wantVector) {
    21982192#ifdef ABC_INHERIT
    2199                     model_->dealWithAbc(1,1);
     2193        model_->dealWithAbc(1, 1);
    22002194#else
    2201                     model_->primal(1);
    2202 #endif
    2203                } else {
    2204                     ClpMatrixBase * matrix = model_->clpMatrix();
    2205                     ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    2206                     assert (clpMatrix);
    2207                     clpMatrix->makeSpecialColumnCopy();
     2195        model_->primal(1);
     2196#endif
     2197      } else {
     2198        ClpMatrixBase *matrix = model_->clpMatrix();
     2199        ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     2200        assert(clpMatrix);
     2201        clpMatrix->makeSpecialColumnCopy();
    22082202#ifdef ABC_INHERIT
    2209                     model_->dealWithAbc(1,1);
     2203        model_->dealWithAbc(1, 1);
    22102204#else
    2211                     model_->primal(1);
    2212 #endif
    2213                     clpMatrix->releaseSpecialColumnCopy();
    2214                }
    2215                if (presolve) {
    2216                     pinfo.postsolve(true);
    2217                     delete model_;
    2218                     model_ = saveModel;
    2219                     saveModel = NULL;
    2220                }
    2221           }
     2205        model_->primal(1);
     2206#endif
     2207        clpMatrix->releaseSpecialColumnCopy();
     2208      }
     2209      if (presolve) {
     2210        pinfo.postsolve(true);
     2211        delete model_;
     2212        model_ = saveModel;
     2213        saveModel = NULL;
     2214      }
     2215    }
    22222216#ifdef COIN_DEVELOP
    2223           printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
    2224 #endif
    2225           delete [] saveUpper;
    2226           delete [] saveLower;
    2227      }
     2217    printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
     2218#endif
     2219    delete[] saveUpper;
     2220    delete[] saveLower;
     2221  }
    22282222#ifdef FEB_TRY
    2229      model_->setSpecialOptions(saveOptions);
    2230      model_->setPerturbation(savePerturbation);
    2231 #endif
    2232      return ;
     2223  model_->setSpecialOptions(saveOptions);
     2224  model_->setPerturbation(savePerturbation);
     2225#endif
     2226  return;
    22332227}
    22342228#endif
     
    22382232Idiot::Idiot()
    22392233{
    2240      model_ = NULL;
    2241      maxBigIts_ = 3;
    2242      maxIts_ = 5;
    2243      logLevel_ = 1;
    2244      logFreq_ = 100;
    2245      maxIts2_ = 100;
    2246      djTolerance_ = 1e-1;
    2247      mu_ = 1e-4;
    2248      drop_ = 5.0;
    2249      exitDrop_ = -1.0e20;
    2250      muFactor_ = 0.3333;
    2251      stopMu_ = 1e-12;
    2252      smallInfeas_ = 1e-1;
    2253      reasonableInfeas_ = 1e2;
    2254      muAtExit_ = 1.0e31;
    2255      strategy_ = 8;
    2256      lambdaIterations_ = 0;
    2257      checkFrequency_ = 100;
    2258      whenUsed_ = NULL;
    2259      majorIterations_ = 30;
    2260      exitFeasibility_ = -1.0;
    2261      dropEnoughFeasibility_ = 0.02;
    2262      dropEnoughWeighted_ = 0.01;
    2263      // adjust
    2264      double nrows = 10000.0;
    2265      int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
    2266      baseIts = baseIts / 10;
    2267      baseIts *= 10;
    2268      maxIts2_ = 200 + baseIts + 5;
    2269      maxIts2_ = 100;
    2270      reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
    2271      lightWeight_ = 0;
     2234  model_ = NULL;
     2235  maxBigIts_ = 3;
     2236  maxIts_ = 5;
     2237  logLevel_ = 1;
     2238  logFreq_ = 100;
     2239  maxIts2_ = 100;
     2240  djTolerance_ = 1e-1;
     2241  mu_ = 1e-4;
     2242  drop_ = 5.0;
     2243  exitDrop_ = -1.0e20;
     2244  muFactor_ = 0.3333;
     2245  stopMu_ = 1e-12;
     2246  smallInfeas_ = 1e-1;
     2247  reasonableInfeas_ = 1e2;
     2248  muAtExit_ = 1.0e31;
     2249  strategy_ = 8;
     2250  lambdaIterations_ = 0;
     2251  checkFrequency_ = 100;
     2252  whenUsed_ = NULL;
     2253  majorIterations_ = 30;
     2254  exitFeasibility_ = -1.0;
     2255  dropEnoughFeasibility_ = 0.02;
     2256  dropEnoughWeighted_ = 0.01;
     2257  // adjust
     2258  double nrows = 10000.0;
     2259  int baseIts = static_cast< int >(sqrt(static_cast< double >(nrows)));
     2260  baseIts = baseIts / 10;
     2261  baseIts *= 10;
     2262  maxIts2_ = 200 + baseIts + 5;
     2263  maxIts2_ = 100;
     2264  reasonableInfeas_ = static_cast< double >(nrows) * 0.05;
     2265  lightWeight_ = 0;
    22722266}
    22732267// Constructor from model
    22742268Idiot::Idiot(OsiSolverInterface &model)
    22752269{
    2276      model_ = & model;
    2277      maxBigIts_ = 3;
    2278      maxIts_ = 5;
    2279      logLevel_ = 1;
    2280      logFreq_ = 100;
    2281      maxIts2_ = 100;
    2282      djTolerance_ = 1e-1;
    2283      mu_ = 1e-4;
    2284      drop_ = 5.0;
    2285      exitDrop_ = -1.0e20;
    2286      muFactor_ = 0.3333;
    2287      stopMu_ = 1e-12;
    2288      smallInfeas_ = 1e-1;
    2289      reasonableInfeas_ = 1e2;
    2290      muAtExit_ = 1.0e31;
    2291      strategy_ = 8;
    2292      lambdaIterations_ = 0;
    2293      checkFrequency_ = 100;
    2294      whenUsed_ = NULL;
    2295      majorIterations_ = 30;
    2296      exitFeasibility_ = -1.0;
    2297      dropEnoughFeasibility_ = 0.02;
    2298      dropEnoughWeighted_ = 0.01;
    2299      // adjust
    2300      double nrows;
    2301      if (model_)
    2302           nrows = model_->getNumRows();
    2303      else
    2304           nrows = 10000.0;
    2305      int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
    2306      baseIts = baseIts / 10;
    2307      baseIts *= 10;
    2308      maxIts2_ = 200 + baseIts + 5;
    2309      maxIts2_ = 100;
    2310      reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
    2311      lightWeight_ = 0;
     2270  model_ = &model;
     2271  maxBigIts_ = 3;
     2272  maxIts_ = 5;
     2273  logLevel_ = 1;
     2274  logFreq_ = 100;
     2275  maxIts2_ = 100;
     2276  djTolerance_ = 1e-1;
     2277  mu_ = 1e-4;
     2278  drop_ = 5.0;
     2279  exitDrop_ = -1.0e20;
     2280  muFactor_ = 0.3333;
     2281  stopMu_ = 1e-12;
     2282  smallInfeas_ = 1e-1;
     2283  reasonableInfeas_ = 1e2;
     2284  muAtExit_ = 1.0e31;
     2285  strategy_ = 8;
     2286  lambdaIterations_ = 0;
     2287  checkFrequency_ = 100;
     2288  whenUsed_ = NULL;
     2289  majorIterations_ = 30;
     2290  exitFeasibility_ = -1.0;
     2291  dropEnoughFeasibility_ = 0.02;
     2292  dropEnoughWeighted_ = 0.01;
     2293  // adjust
     2294  double nrows;
     2295  if (model_)
     2296    nrows = model_->getNumRows();
     2297  else
     2298    nrows = 10000.0;
     2299  int baseIts = static_cast< int >(sqrt(static_cast< double >(nrows)));
     2300  baseIts = baseIts / 10;
     2301  baseIts *= 10;
     2302  maxIts2_ = 200 + baseIts + 5;
     2303  maxIts2_ = 100;
     2304  reasonableInfeas_ = static_cast< double >(nrows) * 0.05;
     2305  lightWeight_ = 0;
    23122306}
    23132307// Copy constructor.
    23142308Idiot::Idiot(const Idiot &rhs)
    23152309{
    2316      model_ = rhs.model_;
    2317      if (model_ && rhs.whenUsed_) {
    2318           int numberColumns = model_->getNumCols();
    2319           whenUsed_ = new int [numberColumns];
    2320           CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
    2321      } else {
    2322           whenUsed_ = NULL;
    2323      }
    2324      djTolerance_ = rhs.djTolerance_;
    2325      mu_ = rhs.mu_;
    2326      drop_ = rhs.drop_;
    2327      muFactor_ = rhs.muFactor_;
    2328      stopMu_ = rhs.stopMu_;
    2329      smallInfeas_ = rhs.smallInfeas_;
    2330      reasonableInfeas_ = rhs.reasonableInfeas_;
    2331      exitDrop_ = rhs.exitDrop_;
    2332      muAtExit_ = rhs.muAtExit_;
    2333      exitFeasibility_ = rhs.exitFeasibility_;
    2334      dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
    2335      dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
    2336      maxBigIts_ = rhs.maxBigIts_;
    2337      maxIts_ = rhs.maxIts_;
    2338      majorIterations_ = rhs.majorIterations_;
    2339      logLevel_ = rhs.logLevel_;
    2340      logFreq_ = rhs.logFreq_;
    2341      checkFrequency_ = rhs.checkFrequency_;
    2342      lambdaIterations_ = rhs.lambdaIterations_;
    2343      maxIts2_ = rhs.maxIts2_;
    2344      strategy_ = rhs.strategy_;
    2345      lightWeight_ = rhs.lightWeight_;
     2310  model_ = rhs.model_;
     2311  if (model_ && rhs.whenUsed_) {
     2312    int numberColumns = model_->getNumCols();
     2313    whenUsed_ = new int[numberColumns];
     2314    CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
     2315  } else {
     2316    whenUsed_ = NULL;
     2317  }
     2318  djTolerance_ = rhs.djTolerance_;
     2319  mu_ = rhs.mu_;
     2320  drop_ = rhs.drop_;
     2321  muFactor_ = rhs.muFactor_;
     2322  stopMu_ = rhs.stopMu_;
     2323  smallInfeas_ = rhs.smallInfeas_;
     2324  reasonableInfeas_ = rhs.reasonableInfeas_;
     2325  exitDrop_ = rhs.exitDrop_;
     2326  muAtExit_ = rhs.muAtExit_;
     2327  exitFeasibility_ = rhs.exitFeasibility_;
     2328  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
     2329  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
     2330  maxBigIts_ = rhs.maxBigIts_;
     2331  maxIts_ = rhs.maxIts_;
     2332  majorIterations_ = rhs.majorIterations_;
     2333  logLevel_ = rhs.logLevel_;
     2334  logFreq_ = rhs.logFreq_;
     2335  checkFrequency_ = rhs.checkFrequency_;
     2336  lambdaIterations_ = rhs.lambdaIterations_;
     2337  maxIts2_ = rhs.maxIts2_;
     2338  strategy_ = rhs.strategy_;
     2339  lightWeight_ = rhs.lightWeight_;
    23462340}
    23472341// Assignment operator. This copies the data
    23482342Idiot &
    2349 Idiot::operator=(const Idiot & rhs)
     2343Idiot::operator=(const Idiot &rhs)
    23502344{
    2351      if (this != &rhs) {
    2352           delete [] whenUsed_;
    2353           model_ = rhs.model_;
    2354           if (model_ && rhs.whenUsed_) {
    2355                int numberColumns = model_->getNumCols();
    2356                whenUsed_ = new int [numberColumns];
    2357                CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
    2358           } else {
    2359                whenUsed_ = NULL;
    2360           }
    2361           djTolerance_ = rhs.djTolerance_;
    2362           mu_ = rhs.mu_;
    2363           drop_ = rhs.drop_;
    2364           muFactor_ = rhs.muFactor_;
    2365           stopMu_ = rhs.stopMu_;
    2366           smallInfeas_ = rhs.smallInfeas_;
    2367           reasonableInfeas_ = rhs.reasonableInfeas_;
    2368           exitDrop_ = rhs.exitDrop_;
    2369           muAtExit_ = rhs.muAtExit_;
    2370           exitFeasibility_ = rhs.exitFeasibility_;
    2371           dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
    2372           dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
    2373           maxBigIts_ = rhs.maxBigIts_;
    2374           maxIts_ = rhs.maxIts_;
    2375           majorIterations_ = rhs.majorIterations_;
    2376           logLevel_ = rhs.logLevel_;
    2377           logFreq_ = rhs.logFreq_;
    2378           checkFrequency_ = rhs.checkFrequency_;
    2379           lambdaIterations_ = rhs.lambdaIterations_;
    2380           maxIts2_ = rhs.maxIts2_;
    2381           strategy_ = rhs.strategy_;
    2382           lightWeight_ = rhs.lightWeight_;
    2383      }
    2384      return *this;
     2345  if (this != &rhs) {
     2346    delete[] whenUsed_;
     2347    model_ = rhs.model_;
     2348    if (model_ && rhs.whenUsed_) {
     2349      int numberColumns = model_->getNumCols();
     2350      whenUsed_ = new int[numberColumns];
     2351      CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
     2352    } else {
     2353      whenUsed_ = NULL;
     2354    }
     2355    djTolerance_ = rhs.djTolerance_;
     2356    mu_ = rhs.mu_;
     2357    drop_ = rhs.drop_;
     2358    muFactor_ = rhs.muFactor_;
     2359    stopMu_ = rhs.stopMu_;
     2360    smallInfeas_ = rhs.smallInfeas_;
     2361    reasonableInfeas_ = rhs.reasonableInfeas_;
     2362    exitDrop_ = rhs.exitDrop_;
     2363    muAtExit_ = rhs.muAtExit_;
     2364    exitFeasibility_ = rhs.exitFeasibility_;
     2365    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
     2366    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
     2367    maxBigIts_ = rhs.maxBigIts_;
     2368    maxIts_ = rhs.maxIts_;
     2369    majorIterations_ = rhs.majorIterations_;
     2370    logLevel_ = rhs.logLevel_;
     2371    logFreq_ = rhs.logFreq_;
     2372    checkFrequency_ = rhs.checkFrequency_;
     2373    lambdaIterations_ = rhs.lambdaIterations_;
     2374    maxIts2_ = rhs.maxIts2_;
     2375    strategy_ = rhs.strategy_;
     2376    lightWeight_ = rhs.lightWeight_;
     2377  }
     2378  return *this;
    23852379}
    23862380Idiot::~Idiot()
    23872381{
    2388      delete [] whenUsed_;
     2382  delete[] whenUsed_;
    23892383}
     2384
     2385/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     2386*/
Note: See TracChangeset for help on using the changeset viewer.