source: trunk/Clp/src/Idiot.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 8 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 74.4 KB
Line 
1/* $Id: Idiot.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include <stdio.h>
8#include <stdarg.h>
9#include <stdlib.h>
10#include <math.h>
11#include "ClpPresolve.hpp"
12#include "Idiot.hpp"
13#include "CoinTime.hpp"
14#include "CoinSort.hpp"
15#include "CoinFactorization.hpp"
16#include "CoinMessageHandler.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "AbcCommon.hpp"
19#include "ClpEventHandler.hpp"
20// Redefine stuff for Clp
21#ifndef OSI_IDIOT
22#include "ClpMessage.hpp"
23#define OsiObjOffset ClpObjOffset
24#endif
25/**** strategy 4 - drop, exitDrop and djTolerance all relative:
26For first two major iterations these are small.  Then:
27
28drop - exit a major iteration if drop over 5*checkFrequency < this is
29used as info->drop*(10.0+fabs(last weighted objective))
30
31exitDrop - exit idiot if feasible and drop < this is
32used as info->exitDrop*(10.0+fabs(last objective))
33
34djExit - exit a major iteration if largest dj (averaged over 5 checks)
35drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective)
36
37djFlag - mostly skip variables with bad dj worse than this => 2*djExit
38
39djTol - only look at variables with dj better than this => 0.01*djExit
40****************/
41
42#define IDIOT_FIX_TOLERANCE 1e-6
43#define SMALL_IDIOT_FIX_TOLERANCE 1e-10
44int Idiot::dropping(IdiotResult result,
45  double tolerance,
46  double small,
47  int *nbad)
48{
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  }
66}
67// Deals with whenUsed and slacks
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)
73{
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 {
90#ifdef COIN_DEVELOP
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));
100#ifdef OSI_IDIOT
101    const CoinPackedMatrix *matrix = model_->getMatrixByCol();
102#else
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  }
284}
285
286/* returns -1 if none or start of costed slacks or -2 if
287   there are costed slacks but it is messy */
288static int countCostedSlacks(OsiSolverInterface *model)
289{
290#ifdef OSI_IDIOT
291  const CoinPackedMatrix *matrix = model->getMatrixByCol();
292#else
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;
307
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;
335}
336void Idiot::crash(int numberPass, CoinMessageHandler *handler,
337  const CoinMessages *messages, bool doCrossover)
338{
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);
376#ifndef OSI_IDIOT
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  }
384#endif
385}
386
387void Idiot::solve()
388{
389  CoinMessages dummy;
390  solve2(NULL, &dummy);
391}
392void Idiot::solve2(CoinMessageHandler *handler, const CoinMessages *messages)
393{
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  }
444#ifdef OSI_IDIOT
445  const CoinPackedMatrix *matrix = model_->getMatrixByCol();
446#else
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;
459#ifndef OSI_IDIOT
460  double *cost = model_->objective();
461  double *lower = model_->columnLower();
462  double *upper = model_->columnUpper();
463#else
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];
485#ifndef OSI_IDIOT
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 {
507#ifdef COIN_DEVELOP
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];
553#ifndef NDEBUG
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];
564#ifndef NDEBUG
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];
589#ifndef NDEBUG
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];
600#ifndef NDEBUG
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;
678#ifndef OSI_IDIOT
679  if ((strategy_ & 32) != 0 && !allOnes) {
680    if (model_->scalingFlag() > 0)
681      scaled = model_->clpMatrix()->scale(model_) == 0;
682    if (scaled) {
683#define IDIOT_SCALE 2
684#ifndef IDIOT_SCALE
685      const double *rowScale = model_->rowScale();
686#else
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      }
725#ifdef IDIOT_SCALE
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]);
733#else
734        double multiplier = rowScale[irow];
735#endif
736#if IDIOT_SCALE > 1 || IDIOT_SCALE < -1
737#define EQUALITY_MULTIPLIER 2
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 */
793
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);
812
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) {
833#ifndef OSI_IDIOT
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);
838#ifndef OSI_IDIOT
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) {
905#ifndef OSI_IDIOT
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);
910#ifndef OSI_IDIOT
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    }
918#ifndef OSI_IDIOT
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) {
934#ifdef CLP_INVESTIGATE
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;
1037#define LAMBDA
1038#ifdef LAMBDA
1039          if ((saveStrategy & 2048) == 0) {
1040            memset(lambda, 0, nrows * sizeof(double));
1041          }
1042#else
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;
1068#ifdef LAMBDA
1069          for (i = 0; i < nrows; i++) {
1070            lambda[i] += scale * rowsol[i];
1071          }
1072#else
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          }
1089#ifdef DEBUG
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        }
1102#ifdef DEBUG
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;
1136#ifndef OSI_IDIOT
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;
1154#else
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  }
1170#endif
1171#define TRYTHIS
1172#ifdef TRYTHIS
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);
1206#if 0
1207     if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) {
1208          printf(
1209               "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
1210               iteration, mu, lastResult.infeas, lastResult.objval, n);
1211     }
1212#endif
1213#ifndef OSI_IDIOT
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
1278#ifndef OSI_IDIOT
1279  CoinMemcpyN(rowsol, nrows, model_->primalRowSolution());
1280  CoinMemcpyN(colsol, ncols, model_->primalColumnSolution());
1281  CoinMemcpyN(pi, nrows, model_->dualRowSolution());
1282  CoinMemcpyN(dj, ncols, model_->dualColumnSolution());
1283#else
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;
1295}
1296#ifndef OSI_IDIOT
1297void Idiot::crossOver(int mode)
1298{
1299  if (lightWeight_ == 2) {
1300    // total failure
1301    model_->allSlackBasis();
1302    return;
1303  }
1304  double fixTolerance = IDIOT_FIX_TOLERANCE;
1305#ifdef COIN_DEVELOP
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_);
1330
1331  int addAll = mode & 7;
1332  int presolve = 0;
1333
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  }
1414#if 0
1415     //maybe put back - but just get feasible ?
1416     // If not many fixed then just exit
1417     int numberFixed = 0;
1418     for (i = ordStart; i < ordEnd; i++) {
1419          if (colsol[i] < lower[i] + fixTolerance)
1420               numberFixed++;
1421          else if (colsol[i] > upper[i] - fixTolerance)
1422               numberFixed++;
1423     }
1424     if (numberFixed < ncols / 3) {
1425          addAll = 3;
1426          presolve = 0;
1427     }
1428#endif
1429#define FEB_TRY
1430#ifdef FEB_TRY
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
1439        0 - chosen,all used, all
1440        1 - chosen, all
1441        2 - all
1442        3 - do not do anything  - maybe basis
1443     */
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 -
1462             you can change second part as much as you want */
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);
1668#ifdef CLP_INVESTIGATE
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;
1678#else
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);
1703#ifdef CLP_INVESTIGATE
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;
1882#define MESS_UP 0
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        }
1943
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
2048#ifdef ABC_INHERIT
2049#ifndef TWO_GOES
2050        model_->dealWithAbc(1, justValuesPass ? 3 : 1);
2051#else
2052        model_->dealWithAbc(1, 1 + 11);
2053#endif
2054#else
2055#ifndef TWO_GOES
2056        model_->primal(justValuesPass ? 2 : 1);
2057#else
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();
2066#ifdef ABC_INHERIT
2067        model_->dealWithAbc(1, 1);
2068#else
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)
2088#ifdef ABC_INHERIT
2089        model_->dealWithAbc(1, 3);
2090#else
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      }
2129#ifdef COIN_DEVELOP
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) {
2144#ifdef ABC_INHERIT
2145        model_->dealWithAbc(1, 1);
2146#else
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();
2154#ifdef ABC_INHERIT
2155        model_->dealWithAbc(1, 1);
2156#else
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;
2180#ifdef COIN_DEVELOP
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) {
2192#ifdef ABC_INHERIT
2193        model_->dealWithAbc(1, 1);
2194#else
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();
2202#ifdef ABC_INHERIT
2203        model_->dealWithAbc(1, 1);
2204#else
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    }
2216#ifdef COIN_DEVELOP
2217    printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
2218#endif
2219    delete[] saveUpper;
2220    delete[] saveLower;
2221  }
2222#ifdef FEB_TRY
2223  model_->setSpecialOptions(saveOptions);
2224  model_->setPerturbation(savePerturbation);
2225#endif
2226  return;
2227}
2228#endif
2229/*****************************************************************************/
2230
2231// Default contructor
2232Idiot::Idiot()
2233{
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;
2266}
2267// Constructor from model
2268Idiot::Idiot(OsiSolverInterface &model)
2269{
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;
2306}
2307// Copy constructor.
2308Idiot::Idiot(const Idiot &rhs)
2309{
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_;
2340}
2341// Assignment operator. This copies the data
2342Idiot &
2343Idiot::operator=(const Idiot &rhs)
2344{
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;
2379}
2380Idiot::~Idiot()
2381{
2382  delete[] whenUsed_;
2383}
2384
2385/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2386*/
Note: See TracBrowser for help on using the repository browser.