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

Last change on this file since 1910 was 1910, checked in by stefan, 7 years ago
  • add configure option --enable-aboca={1,2,3,4,yes,no}
  • compile Aboca source only if --enable-aboca set (instead of compiling empty source files)
  • fix svn properties
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 80.3 KB
Line 
1/* $Id: Idiot.cpp 1910 2013-01-27 02:00:13Z 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 "CoinMessageHandler.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "AbcCommon.hpp"
18// Redefine stuff for Clp
19#ifndef OSI_IDIOT
20#include "ClpMessage.hpp"
21#define OsiObjOffset ClpObjOffset
22#endif
23/**** strategy 4 - drop, exitDrop and djTolerance all relative:
24For first two major iterations these are small.  Then:
25
26drop - exit a major iteration if drop over 5*checkFrequency < this is
27used as info->drop*(10.0+fabs(last weighted objective))
28
29exitDrop - exit idiot if feasible and drop < this is
30used as info->exitDrop*(10.0+fabs(last objective))
31
32djExit - exit a major iteration if largest dj (averaged over 5 checks)
33drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective)
34
35djFlag - mostly skip variables with bad dj worse than this => 2*djExit
36
37djTol - only look at variables with dj better than this => 0.01*djExit
38****************/
39
40#define IDIOT_FIX_TOLERANCE 1e-6
41#define SMALL_IDIOT_FIX_TOLERANCE 1e-10
42int
43Idiot::dropping(IdiotResult result,
44                double tolerance,
45                double small,
46                int *nbad)
47{
48     if (result.infeas <= small) {
49          double value = CoinMax(fabs(result.objval), fabs(result.dropThis)) + 1.0;
50          if (result.dropThis > tolerance * value) {
51               *nbad = 0;
52               return 1;
53          } else {
54               (*nbad)++;
55               if (*nbad > 4) {
56                    return 0;
57               } else {
58                    return 1;
59               }
60          }
61     } else {
62          *nbad = 0;
63          return 1;
64     }
65}
66// Deals with whenUsed and slacks
67int
68Idiot::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)
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          ClpMatrixBase * matrix = model_->clpMatrix();
104#endif
105          const int * row = matrix->getIndices();
106          const CoinBigIndex * columnStart = matrix->getVectorStarts();
107          const int * columnLength = matrix->getVectorLengths();
108          //const double * element = matrix->getElements();
109          int i;
110          objValue = 0.0;
111          infValue = 0.0;
112          for ( i = 0; i < ncols; i++) {
113               if (nextSlack[i] == -1) {
114                    // not a slack
115                    if (colsol[i] > lower[i] + fixTolerance) {
116                         if (colsol[i] < upper[i] - fixTolerance) {
117                              n++;
118                              whenUsed_[i] = iteration;
119                         } else {
120                              colsol[i] = upper[i];
121                         }
122                         whenUsed_[i] = iteration;
123                    } else {
124                         colsol[i] = lower[i];
125                    }
126                    double value = colsol[i];
127                    if (value) {
128                         objValue += cost[i] * value;
129                         CoinBigIndex j;
130                         for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
131                              int iRow = row[j];
132                              rowsol[iRow] += value * element[j];
133                         }
134                    }
135               }
136          }
137          // temp fix for infinite lbs - just limit to -1000
138          for (i = 0; i < nrows; i++) {
139               double rowSave = rowsol[i];
140               int iCol;
141               iCol = posSlack[i];
142               if (iCol >= 0) {
143                    // slide all slack down
144                    double rowValue = rowsol[i];
145                    CoinBigIndex j = columnStart[iCol];
146                    double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
147                    rowSave += (colsol[iCol] - lowerValue) * element[j];
148                    colsol[iCol] = lowerValue;
149                    while (nextSlack[iCol] >= 0) {
150                         iCol = nextSlack[iCol];
151                         double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
152                         j = columnStart[iCol];
153                         rowSave += (colsol[iCol] - lowerValue) * element[j];
154                         colsol[iCol] = lowerValue;
155                    }
156                    iCol = posSlack[i];
157                    while (rowValue < rowLower[i] && iCol >= 0) {
158                         // want to increase
159                         double distance = rowLower[i] - rowValue;
160                         double value = element[columnStart[iCol]];
161                         double thisCost = cost[iCol];
162                         if (distance <= value*(upper[iCol] - colsol[iCol])) {
163                              // can get there
164                              double movement = distance / value;
165                              objValue += movement * thisCost;
166                              rowValue = rowLower[i];
167                              colsol[iCol] += movement;
168                         } else {
169                              // can't get there
170                              double movement = upper[iCol] - colsol[iCol];
171                              objValue += movement * thisCost;
172                              rowValue += movement * value;
173                              colsol[iCol] = upper[iCol];
174                              iCol = nextSlack[iCol];
175                         }
176                    }
177                    if (iCol >= 0) {
178                         // may want to carry on - because of cost?
179                         while (iCol >= 0 && cost[iCol] < 0 && rowValue < rowUpper[i]) {
180                              // want to increase
181                              double distance = rowUpper[i] - rowValue;
182                              double value = element[columnStart[iCol]];
183                              double thisCost = cost[iCol];
184                              if (distance <= value*(upper[iCol] - colsol[iCol])) {
185                                   // can get there
186                                   double movement = distance / value;
187                                   objValue += movement * thisCost;
188                                   rowValue = rowUpper[i];
189                                   colsol[iCol] += movement;
190                                   iCol = -1;
191                              } else {
192                                   // can't get there
193                                   double movement = upper[iCol] - colsol[iCol];
194                                   objValue += movement * thisCost;
195                                   rowValue += movement * value;
196                                   colsol[iCol] = upper[iCol];
197                                   iCol = nextSlack[iCol];
198                              }
199                         }
200                         if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance &&
201                                   colsol[iCol] < upper[iCol] - fixTolerance) {
202                              whenUsed_[i] = iteration;
203                              n++;
204                         }
205                    }
206                    rowsol[i] = rowValue;
207               }
208               iCol = negSlack[i];
209               if (iCol >= 0) {
210                    // slide all slack down
211                    double rowValue = rowsol[i];
212                    CoinBigIndex j = columnStart[iCol];
213                    double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
214                    rowSave += (colsol[iCol] - lowerValue) * element[j];
215                    colsol[iCol] = lowerValue;
216                    while (nextSlack[iCol] >= 0) {
217                         iCol = nextSlack[iCol];
218                         j = columnStart[iCol];
219                         double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
220                         rowSave += (colsol[iCol] - lowerValue) * element[j];
221                         colsol[iCol] = lowerValue;
222                    }
223                    iCol = negSlack[i];
224                    while (rowValue > rowUpper[i] && iCol >= 0) {
225                         // want to increase
226                         double distance = -(rowUpper[i] - rowValue);
227                         double value = -element[columnStart[iCol]];
228                         double thisCost = cost[iCol];
229                         if (distance <= value*(upper[iCol] - lower[iCol])) {
230                              // can get there
231                              double movement = distance / value;
232                              objValue += movement * thisCost;
233                              rowValue = rowUpper[i];
234                              colsol[iCol] += movement;
235                         } else {
236                              // can't get there
237                              double movement = upper[iCol] - lower[iCol];
238                              objValue += movement * thisCost;
239                              rowValue -= movement * value;
240                              colsol[iCol] = upper[iCol];
241                              iCol = nextSlack[iCol];
242                         }
243                    }
244                    if (iCol >= 0) {
245                         // may want to carry on - because of cost?
246                         while (iCol >= 0 && cost[iCol] < 0 && rowValue > rowLower[i]) {
247                              // want to increase
248                              double distance = -(rowLower[i] - rowValue);
249                              double value = -element[columnStart[iCol]];
250                              double thisCost = cost[iCol];
251                              if (distance <= value*(upper[iCol] - colsol[iCol])) {
252                                   // can get there
253                                   double movement = distance / value;
254                                   objValue += movement * thisCost;
255                                   rowValue = rowLower[i];
256                                   colsol[iCol] += movement;
257                                   iCol = -1;
258                              } else {
259                                   // can't get there
260                                   double movement = upper[iCol] - colsol[iCol];
261                                   objValue += movement * thisCost;
262                                   rowValue -= movement * value;
263                                   colsol[iCol] = upper[iCol];
264                                   iCol = nextSlack[iCol];
265                              }
266                         }
267                         if (iCol >= 0 && colsol[iCol] > lower[iCol] + fixTolerance &&
268                                   colsol[iCol] < upper[iCol] - fixTolerance) {
269                              whenUsed_[i] = iteration;
270                              n++;
271                         }
272                    }
273                    rowsol[i] = rowValue;
274               }
275               infValue += CoinMax(CoinMax(0.0, rowLower[i] - rowsol[i]), rowsol[i] - rowUpper[i]);
276               // just change
277               rowsol[i] -= rowSave;
278          }
279          return n;
280     }
281}
282
283/* returns -1 if none or start of costed slacks or -2 if
284   there are costed slacks but it is messy */
285static int countCostedSlacks(OsiSolverInterface * model)
286{
287#ifdef OSI_IDIOT
288     const CoinPackedMatrix * matrix = model->getMatrixByCol();
289#else
290     ClpMatrixBase * matrix = model->clpMatrix();
291#endif
292     const int * row = matrix->getIndices();
293     const CoinBigIndex * columnStart = matrix->getVectorStarts();
294     const int * columnLength = matrix->getVectorLengths();
295     const double * element = matrix->getElements();
296     const double * rowupper = model->getRowUpper();
297     int nrows = model->getNumRows();
298     int ncols = model->getNumCols();
299     int slackStart = ncols - nrows;
300     int nSlacks = nrows;
301     int i;
302
303     if (ncols <= nrows) return -1;
304     while (1) {
305          for (i = 0; i < nrows; i++) {
306               int j = i + slackStart;
307               CoinBigIndex k = columnStart[j];
308               if (columnLength[j] == 1) {
309                    if (row[k] != i || element[k] != 1.0) {
310                         nSlacks = 0;
311                         break;
312                    }
313               } else {
314                    nSlacks = 0;
315                    break;
316               }
317               if (rowupper[i] <= 0.0) {
318                    nSlacks = 0;
319                    break;
320               }
321          }
322          if (nSlacks || !slackStart) break;
323          slackStart = 0;
324     }
325     if (!nSlacks) slackStart = -1;
326     return slackStart;
327}
328void
329Idiot::crash(int numberPass, CoinMessageHandler * handler,
330             const CoinMessages *messages, bool doCrossover)
331{
332     // lightweight options
333     int numberColumns = model_->getNumCols();
334     const double * objective = model_->getObjCoefficients();
335     int nnzero = 0;
336     double sum = 0.0;
337     int i;
338     for (i = 0; i < numberColumns; i++) {
339          if (objective[i]) {
340               sum += fabs(objective[i]);
341               nnzero++;
342          }
343     }
344     sum /= static_cast<double> (nnzero + 1);
345     if (maxIts_ == 5)
346          maxIts_ = 2;
347     if (numberPass <= 0)
348          majorIterations_ = static_cast<int>(2 + log10(static_cast<double>(numberColumns + 1)));
349     else
350          majorIterations_ = numberPass;
351     // If mu not changed then compute
352     if (mu_ == 1e-4)
353          mu_ = CoinMax(1.0e-3, sum * 1.0e-5);
354     if (maxIts2_ == 100) {
355          if (!lightWeight_) {
356               maxIts2_ = 105;
357          } else if (lightWeight_ == 1) {
358               mu_ *= 1000.0;
359               maxIts2_ = 23;
360          } else if (lightWeight_ == 2) {
361               maxIts2_ = 11;
362          } else {
363               maxIts2_ = 23;
364          }
365     }
366     //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
367     solve2(handler, messages);
368#ifndef OSI_IDIOT
369     if (doCrossover) {
370          double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
371          if ((averageInfeas < 0.01 && (strategy_ & 512) != 0) || (strategy_ & 8192) != 0)
372               crossOver(16 + 1);
373          else
374               crossOver(majorIterations_ < 1000000 ? 3 : 2);
375     }
376#endif
377}
378
379void
380Idiot::solve()
381{
382     CoinMessages dummy;
383     solve2(NULL, &dummy);
384}
385void
386Idiot::solve2(CoinMessageHandler * handler, const CoinMessages * messages)
387{
388     int strategy = 0;
389     double d2;
390     int i, n;
391     int allOnes = 1;
392     int iteration = 0;
393     int iterationTotal = 0;
394     int nTry = 0; /* number of tries at same weight */
395     double fixTolerance = IDIOT_FIX_TOLERANCE;
396     int maxBigIts = maxBigIts_;
397     int maxIts = maxIts_;
398     int logLevel = logLevel_;
399     int saveMajorIterations = majorIterations_;
400     majorIterations_ = majorIterations_ % 1000000;
401     if (handler) {
402          if (handler->logLevel() > 0 && handler->logLevel() < 3)
403               logLevel = 1;
404          else if (!handler->logLevel())
405               logLevel = 0;
406          else
407               logLevel = 7;
408     }
409     double djExit = djTolerance_;
410     double djFlag = 1.0 + 100.0 * djExit;
411     double djTol = 0.00001;
412     double mu = mu_;
413     double drop = drop_;
414     int maxIts2 = maxIts2_;
415     double factor = muFactor_;
416     double smallInfeas = smallInfeas_;
417     double reasonableInfeas = reasonableInfeas_;
418     double stopMu = stopMu_;
419     double maxmin, offset;
420     double lastWeighted = 1.0e50;
421     double exitDrop = exitDrop_;
422     double fakeSmall = smallInfeas;
423     double firstInfeas;
424     int badIts = 0;
425     int slackStart, slackEnd, ordStart, ordEnd;
426     int checkIteration = 0;
427     int lambdaIteration = 0;
428     int belowReasonable = 0; /* set if ever gone below reasonable infeas */
429     double bestWeighted = 1.0e60;
430     double bestFeasible = 1.0e60; /* best solution while feasible */
431     IdiotResult result, lastResult;
432     int saveStrategy = strategy_;
433     const int strategies[] = {0, 2, 128};
434     double saveLambdaScale = 0.0;
435     if ((saveStrategy & 128) != 0) {
436          fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
437     }
438#ifdef OSI_IDIOT
439     const CoinPackedMatrix * matrix = model_->getMatrixByCol();
440#else
441     ClpMatrixBase * matrix = model_->clpMatrix();
442#endif
443     const int * row = matrix->getIndices();
444     const CoinBigIndex * columnStart = matrix->getVectorStarts();
445     const int * columnLength = matrix->getVectorLengths();
446     const double * element = matrix->getElements();
447     int nrows = model_->getNumRows();
448     int ncols = model_->getNumCols();
449     double * rowsol, * colsol;
450     double * pi, * dj;
451#ifndef OSI_IDIOT
452     double * cost = model_->objective();
453     double * lower = model_->columnLower();
454     double * upper = model_->columnUpper();
455#else
456     double * cost = new double [ncols];
457     CoinMemcpyN( model_->getObjCoefficients(), ncols, cost);
458     const double * lower = model_->getColLower();
459     const double * upper = model_->getColUpper();
460#endif
461     const double *elemXX;
462     double * saveSol;
463     double * rowupper = new double[nrows]; // not const as modified
464     CoinMemcpyN(model_->getRowUpper(), nrows, rowupper);
465     double * rowlower = new double[nrows]; // not const as modified
466     CoinMemcpyN(model_->getRowLower(), nrows, rowlower);
467     CoinThreadRandom * randomNumberGenerator = model_->randomNumberGenerator();
468     int * whenUsed;
469     double * lambda;
470     saveSol = new double[ncols];
471     lambda = new double [nrows];
472     rowsol = new double[nrows];
473     colsol = new double [ncols];
474     CoinMemcpyN(model_->getColSolution(), ncols, colsol);
475     pi = new double[nrows];
476     dj = new double[ncols];
477     delete [] whenUsed_;
478     bool oddSlacks = false;
479     // See if any costed slacks
480     int numberSlacks = 0;
481     for (i = 0; i < ncols; i++) {
482          if (columnLength[i] == 1)
483               numberSlacks++;
484     }
485     if (!numberSlacks) {
486          whenUsed_ = new int[ncols];
487     } else {
488#ifdef COIN_DEVELOP
489          printf("%d slacks\n", numberSlacks);
490#endif
491          oddSlacks = true;
492          int extra = static_cast<int> (nrows * sizeof(double) / sizeof(int));
493          whenUsed_ = new int[2*ncols+2*nrows+extra];
494          int * posSlack = whenUsed_ + ncols;
495          int * negSlack = posSlack + nrows;
496          int * nextSlack = negSlack + nrows;
497          for (i = 0; i < nrows; i++) {
498               posSlack[i] = -1;
499               negSlack[i] = -1;
500          }
501          for (i = 0; i < ncols; i++)
502               nextSlack[i] = -1;
503          for (i = 0; i < ncols; i++) {
504               if (columnLength[i] == 1) {
505                    CoinBigIndex j = columnStart[i];
506                    int iRow = row[j];
507                    if (element[j] > 0.0) {
508                         if (posSlack[iRow] == -1) {
509                              posSlack[iRow] = i;
510                         } else {
511                              int iCol = posSlack[iRow];
512                              while (nextSlack[iCol] >= 0)
513                                   iCol = nextSlack[iCol];
514                              nextSlack[iCol] = i;
515                         }
516                    } else {
517                         if (negSlack[iRow] == -1) {
518                              negSlack[iRow] = i;
519                         } else {
520                              int iCol = negSlack[iRow];
521                              while (nextSlack[iCol] >= 0)
522                                   iCol = nextSlack[iCol];
523                              nextSlack[iCol] = i;
524                         }
525                    }
526               }
527          }
528          // now sort
529          for (i = 0; i < nrows; i++) {
530               int iCol;
531               iCol = posSlack[i];
532               if (iCol >= 0) {
533                    CoinBigIndex j = columnStart[iCol];
534#ifndef NDEBUG
535                    int iRow = row[j];
536#endif
537                    assert (element[j] > 0.0);
538                    assert (iRow == i);
539                    dj[0] = cost[iCol] / element[j];
540                    whenUsed_[0] = iCol;
541                    int n = 1;
542                    while (nextSlack[iCol] >= 0) {
543                         iCol = nextSlack[iCol];
544                         CoinBigIndex j = columnStart[iCol];
545#ifndef NDEBUG
546                         int iRow = row[j];
547#endif
548                         assert (element[j] > 0.0);
549                         assert (iRow == i);
550                         dj[n] = cost[iCol] / element[j];
551                         whenUsed_[n++] = iCol;
552                    }
553                    for (j = 0; j < n; j++) {
554                         int jCol = whenUsed_[j];
555                         nextSlack[jCol] = -2;
556                    }
557                    CoinSort_2(dj, dj + n, whenUsed_);
558                    // put back
559                    iCol = whenUsed_[0];
560                    posSlack[i] = iCol;
561                    for (j = 1; j < n; j++) {
562                         int jCol = whenUsed_[j];
563                         nextSlack[iCol] = jCol;
564                         iCol = jCol;
565                    }
566               }
567               iCol = negSlack[i];
568               if (iCol >= 0) {
569                    CoinBigIndex j = columnStart[iCol];
570#ifndef NDEBUG
571                    int iRow = row[j];
572#endif
573                    assert (element[j] < 0.0);
574                    assert (iRow == i);
575                    dj[0] = -cost[iCol] / element[j];
576                    whenUsed_[0] = iCol;
577                    int n = 1;
578                    while (nextSlack[iCol] >= 0) {
579                         iCol = nextSlack[iCol];
580                         CoinBigIndex j = columnStart[iCol];
581#ifndef NDEBUG
582                         int iRow = row[j];
583#endif
584                         assert (element[j] < 0.0);
585                         assert (iRow == i);
586                         dj[n] = -cost[iCol] / element[j];
587                         whenUsed_[n++] = iCol;
588                    }
589                    for (j = 0; j < n; j++) {
590                         int jCol = whenUsed_[j];
591                         nextSlack[jCol] = -2;
592                    }
593                    CoinSort_2(dj, dj + n, whenUsed_);
594                    // put back
595                    iCol = whenUsed_[0];
596                    negSlack[i] = iCol;
597                    for (j = 1; j < n; j++) {
598                         int jCol = whenUsed_[j];
599                         nextSlack[iCol] = jCol;
600                         iCol = jCol;
601                    }
602               }
603          }
604     }
605     whenUsed = whenUsed_;
606     if (model_->getObjSense() == -1.0) {
607          maxmin = -1.0;
608     } else {
609          maxmin = 1.0;
610     }
611     model_->getDblParam(OsiObjOffset, offset);
612     if (!maxIts2) maxIts2 = maxIts;
613     strategy = strategy_;
614     strategy &= 3;
615     memset(lambda, 0, nrows * sizeof(double));
616     slackStart = countCostedSlacks(model_);
617     if (slackStart >= 0) {
618       COIN_DETAIL_PRINT(printf("This model has costed slacks\n"));
619          slackEnd = slackStart + nrows;
620          if (slackStart) {
621               ordStart = 0;
622               ordEnd = slackStart;
623          } else {
624               ordStart = nrows;
625               ordEnd = ncols;
626          }
627     } else {
628          slackEnd = slackStart;
629          ordStart = 0;
630          ordEnd = ncols;
631     }
632     if (offset && logLevel > 2) {
633          printf("** Objective offset is %g\n", offset);
634     }
635     /* compute reasonable solution cost */
636     for (i = 0; i < nrows; i++) {
637          rowsol[i] = 1.0e31;
638     }
639     for (i = 0; i < ncols; i++) {
640          CoinBigIndex j;
641          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
642               if (element[j] != 1.0) {
643                    allOnes = 0;
644                    break;
645               }
646          }
647     }
648     if (allOnes) {
649          elemXX = NULL;
650     } else {
651          elemXX = element;
652     }
653     // Do scaling if wanted
654     bool scaled = false;
655#ifndef OSI_IDIOT
656     if ((strategy_ & 32) != 0 && !allOnes) {
657          if (model_->scalingFlag() > 0)
658               scaled = model_->clpMatrix()->scale(model_) == 0;
659          if (scaled) {
660               const double * rowScale = model_->rowScale();
661               const double * columnScale = model_->columnScale();
662               double * oldLower = lower;
663               double * oldUpper = upper;
664               double * oldCost = cost;
665               lower = new double[ncols];
666               upper = new double[ncols];
667               cost = new double[ncols];
668               CoinMemcpyN(oldLower, ncols, lower);
669               CoinMemcpyN(oldUpper, ncols, upper);
670               CoinMemcpyN(oldCost, ncols, cost);
671               int icol, irow;
672               for (icol = 0; icol < ncols; icol++) {
673                    double multiplier = 1.0 / columnScale[icol];
674                    if (lower[icol] > -1.0e50)
675                         lower[icol] *= multiplier;
676                    if (upper[icol] < 1.0e50)
677                         upper[icol] *= multiplier;
678                    colsol[icol] *= multiplier;
679                    cost[icol] *= columnScale[icol];
680               }
681               CoinMemcpyN(model_->rowLower(), nrows, rowlower);
682               for (irow = 0; irow < nrows; irow++) {
683                    double multiplier = rowScale[irow];
684                    if (rowlower[irow] > -1.0e50)
685                         rowlower[irow] *= multiplier;
686                    if (rowupper[irow] < 1.0e50)
687                         rowupper[irow] *= multiplier;
688                    rowsol[irow] *= multiplier;
689               }
690               int length = columnStart[ncols-1] + columnLength[ncols-1];
691               double * elemYY = new double[length];
692               for (i = 0; i < ncols; i++) {
693                    CoinBigIndex j;
694                    double scale = columnScale[i];
695                    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
696                         int irow = row[j];
697                         elemYY[j] = element[j] * scale * rowScale[irow];
698                    }
699               }
700               elemXX = elemYY;
701          }
702     }
703#endif
704     for (i = 0; i < ncols; i++) {
705          CoinBigIndex j;
706          double dd = columnLength[i];
707          dd = cost[i] / dd;
708          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
709               int irow = row[j];
710               if (dd < rowsol[irow]) {
711                    rowsol[irow] = dd;
712               }
713          }
714     }
715     d2 = 0.0;
716     for (i = 0; i < nrows; i++) {
717          d2 += rowsol[i];
718     }
719     d2 *= 2.0; /* for luck */
720
721     d2 = d2 / static_cast<double> (4 * nrows + 8000);
722     d2 *= 0.5; /* halve with more flexible method */
723     if (d2 < 5.0) d2 = 5.0;
724     if (djExit == 0.0) {
725          djExit = d2;
726     }
727     if ((saveStrategy & 4) != 0) {
728          /* go to relative tolerances - first small */
729          djExit = 1.0e-10;
730          djFlag = 1.0e-5;
731          drop = 1.0e-10;
732     }
733     memset(whenUsed, 0, ncols * sizeof(int));
734     strategy = strategies[strategy];
735     if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */
736     CoinMemcpyN(colsol, ncols, saveSol);
737
738     lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi,
739                           dj, cost, rowlower, rowupper,
740                           lower, upper, elemXX, row, columnStart, columnLength, lambda,
741                           0, mu, drop,
742                           maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
743     // update whenUsed_
744     n = cleanIteration(iteration, ordStart, ordEnd,
745                        colsol,  lower,  upper,
746                        rowlower, rowupper,
747                        cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas);
748     if ((strategy_ & 16384) != 0) {
749          int * posSlack = whenUsed_ + ncols;
750          int * negSlack = posSlack + nrows;
751          int * nextSlack = negSlack + nrows;
752          double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
753          for (i = 0; i < nrows; i++)
754               rowsol[i] += rowsol2[i];
755     }
756     if ((logLevel_ & 1) != 0) {
757#ifndef OSI_IDIOT
758          if (!handler) {
759#endif
760               printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
761                      iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n);
762#ifndef OSI_IDIOT
763          } else {
764               handler->message(CLP_IDIOT_ITERATION, *messages)
765                         << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n
766                         << CoinMessageEol;
767          }
768#endif
769     }
770     int numberBaseTrys = 0; // for first time
771     int numberAway = -1;
772     iterationTotal = lastResult.iteration;
773     firstInfeas = lastResult.infeas;
774     if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas;
775     if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas;
776     double keepinfeas = 1.0e31;
777     double lastInfeas = 1.0e31;
778     double bestInfeas = 1.0e31;
779     while ((mu > stopMu && lastResult.infeas > smallInfeas) ||
780               (lastResult.infeas <= smallInfeas &&
781                dropping(lastResult, exitDrop, smallInfeas, &badIts)) ||
782               checkIteration < 2 || lambdaIteration < lambdaIterations_) {
783          if (lastResult.infeas <= exitFeasibility_)
784               break;
785          iteration++;
786          checkIteration++;
787          if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
788               bestFeasible = lastResult.objval;
789          }
790          if (lastResult.infeas + mu * lastResult.objval < bestWeighted) {
791               bestWeighted = lastResult.objval + mu * lastResult.objval;
792          }
793          if ((saveStrategy & 4096)) strategy |= 256;
794          if ((saveStrategy & 4) != 0 && iteration > 2) {
795               /* go to relative tolerances */
796               double weighted = 10.0 + fabs(lastWeighted);
797               djExit = djTolerance_ * weighted;
798               djFlag = 2.0 * djExit;
799               drop = drop_ * weighted;
800               djTol = 0.01 * djExit;
801          }
802          CoinMemcpyN(colsol, ncols, saveSol);
803          result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj,
804                            cost, rowlower, rowupper,
805                            lower, upper, elemXX, row, columnStart, columnLength, lambda,
806                            maxIts, mu, drop,
807                            maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
808          n = cleanIteration(iteration, ordStart, ordEnd,
809                             colsol,  lower,  upper,
810                             rowlower, rowupper,
811                             cost, elemXX, fixTolerance, result.objval, result.infeas);
812          if ((strategy_ & 16384) != 0) {
813               int * posSlack = whenUsed_ + ncols;
814               int * negSlack = posSlack + nrows;
815               int * nextSlack = negSlack + nrows;
816               double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
817               for (i = 0; i < nrows; i++)
818                    rowsol[i] += rowsol2[i];
819          }
820          if ((logLevel_ & 1) != 0) {
821#ifndef OSI_IDIOT
822               if (!handler) {
823#endif
824                    printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
825                           iteration, result.infeas, result.objval, mu, result.iteration, n);
826#ifndef OSI_IDIOT
827               } else {
828                    handler->message(CLP_IDIOT_ITERATION, *messages)
829                              << iteration << result.infeas << result.objval << mu << result.iteration << n
830                              << CoinMessageEol;
831               }
832#endif
833          }
834          if (iteration > 50 && n == numberAway ) {
835            if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) {
836#ifdef CLP_INVESTIGATE
837               printf("infeas small %g\n", result.infeas);
838#endif
839               break; // not much happening
840            }
841          }
842          if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
843               if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas)
844                    majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible
845          }
846          lastInfeas = result.infeas;
847          numberAway = n;
848          keepinfeas = result.infeas;
849          lastWeighted = result.weighted;
850          iterationTotal += result.iteration;
851          if (iteration == 1) {
852               if ((strategy_ & 1024) != 0 && mu < 1.0e-10)
853                    result.infeas = firstInfeas * 0.8;
854               if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0)
855                    result.infeas *= 0.8;
856               if (result.infeas > firstInfeas * 0.9
857                         && result.infeas > reasonableInfeas) {
858                    iteration--;
859                    if (majorIterations_ < 50)
860                         mu *= 1.0e-1;
861                    else
862                         mu *= 0.7;
863                    bestFeasible = 1.0e31;
864                    bestWeighted = 1.0e60;
865                    numberBaseTrys++;
866                    if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) {
867                         // back to all slack basis
868                         lightWeight_ = 2;
869                         break;
870                    }
871                    CoinMemcpyN(saveSol, ncols, colsol);
872               } else {
873                    maxIts = maxIts2;
874                    checkIteration = 0;
875                    if ((strategy_ & 1024) != 0) mu *= 1.0e-1;
876               }
877          } else {
878          }
879          bestInfeas = CoinMin(bestInfeas, result.infeas);
880          if (majorIterations_>100&&majorIterations_<200) {
881            if (iteration==majorIterations_-100) {
882              // redo
883              double muX=mu*10.0;
884              bestInfeas=1.0e3;
885              mu=muX;
886              nTry=0;
887            }
888          }
889          if (iteration) {
890               /* this code is in to force it to terminate sometime */
891               double changeMu = factor;
892               if ((saveStrategy & 64) != 0) {
893                    keepinfeas = 0.0; /* switch off ranga's increase */
894                    fakeSmall = smallInfeas;
895               } else {
896                    fakeSmall = -1.0;
897               }
898               saveLambdaScale = 0.0;
899               if (result.infeas > reasonableInfeas ||
900                         (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) {
901                    if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) ||
902                              nTry + 1 == maxBigIts ||
903                              (result.infeas > lastResult.infeas * 0.9
904                               && result.weighted > lastResult.weighted
905                               - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) {
906                         mu *= changeMu;
907                         if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) {
908                              reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu));
909                              COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n", reasonableInfeas));
910                         }
911                         result.weighted = 1.0e60;
912                         nTry = 0;
913                         bestFeasible = 1.0e31;
914                         bestWeighted = 1.0e60;
915                         checkIteration = 0;
916                         lambdaIteration = 0;
917#define LAMBDA
918#ifdef LAMBDA
919                         if ((saveStrategy & 2048) == 0) {
920                              memset(lambda, 0, nrows * sizeof(double));
921                         }
922#else
923                         memset(lambda, 0, nrows * sizeof(double));
924#endif
925                    } else {
926                         nTry++;
927                    }
928               } else if (lambdaIterations_ >= 0) {
929                    /* update lambda  */
930                    double scale = 1.0 / mu;
931                    int i, nnz = 0;
932                    saveLambdaScale = scale;
933                    lambdaIteration++;
934                    if ((saveStrategy & 4) == 0) drop = drop_ / 50.0;
935                    if (lambdaIteration > 4 &&
936                              (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) ||
937                               ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) {
938                         //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
939                         smallInfeas *= 1.5;
940                    }
941                    if ((saveStrategy & 2048) == 0) {
942                         for (i = 0; i < nrows; i++) {
943                              if (lambda[i]) nnz++;
944                              lambda[i] += scale * rowsol[i];
945                         }
946                    } else {
947                         nnz = 1;
948#ifdef LAMBDA
949                         for (i = 0; i < nrows; i++) {
950                              lambda[i] += scale * rowsol[i];
951                         }
952#else
953                         for (i = 0; i < nrows; i++) {
954                              lambda[i] = scale * rowsol[i];
955                         }
956                         for (i = 0; i < ncols; i++) {
957                              CoinBigIndex j;
958                              double value = cost[i] * maxmin;
959                              for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
960                                   int irow = row[j];
961                                   value += element[j] * lambda[irow];
962                              }
963                              cost[i] = value * maxmin;
964                         }
965                         for (i = 0; i < nrows; i++) {
966                              offset += lambda[i] * rowupper[i];
967                              lambda[i] = 0.0;
968                         }
969#ifdef DEBUG
970                         printf("offset %g\n", offset);
971#endif
972                         model_->setDblParam(OsiObjOffset, offset);
973#endif
974                    }
975                    nTry++;
976                    if (!nnz) {
977                         bestFeasible = 1.0e32;
978                         bestWeighted = 1.0e60;
979                         checkIteration = 0;
980                         result.weighted = 1.0e31;
981                    }
982#ifdef DEBUG
983                    double trueCost = 0.0;
984                    for (i = 0; i < ncols; i++) {
985                         int j;
986                         trueCost += cost[i] * colsol[i];
987                    }
988                    printf("True objective %g\n", trueCost - offset);
989#endif
990               } else {
991                    nTry++;
992               }
993               lastResult = result;
994               if (result.infeas < reasonableInfeas && !belowReasonable) {
995                    belowReasonable = 1;
996                    bestFeasible = 1.0e32;
997                    bestWeighted = 1.0e60;
998                    checkIteration = 0;
999                    result.weighted = 1.0e31;
1000               }
1001          }
1002          if (iteration >= majorIterations_) {
1003               // If not feasible and crash then dive dive dive
1004               if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) {
1005                    mu = 1.0e-30;
1006                    majorIterations_ = iteration + 1;
1007                    stopMu = 0.0;
1008               } else {
1009                    if (logLevel > 2)
1010                         printf("Exiting due to number of major iterations\n");
1011                    break;
1012               }
1013          }
1014     }
1015     majorIterations_ = saveMajorIterations;
1016#ifndef OSI_IDIOT
1017     if (scaled) {
1018          // Scale solution and free arrays
1019          const double * rowScale = model_->rowScale();
1020          const double * columnScale = model_->columnScale();
1021          int icol, irow;
1022          for (icol = 0; icol < ncols; icol++) {
1023               colsol[icol] *= columnScale[icol];
1024               saveSol[icol] *= columnScale[icol];
1025               dj[icol] /= columnScale[icol];
1026          }
1027          for (irow = 0; irow < nrows; irow++) {
1028               rowsol[irow] /= rowScale[irow];
1029               pi[irow] *= rowScale[irow];
1030          }
1031          // Don't know why getting Microsoft problems
1032#if defined (_MSC_VER)
1033          delete [] ( double *) elemXX;
1034#else
1035          delete [] elemXX;
1036#endif
1037          model_->setRowScale(NULL);
1038          model_->setColumnScale(NULL);
1039          delete [] lower;
1040          delete [] upper;
1041          delete [] cost;
1042          lower = model_->columnLower();
1043          upper = model_->columnUpper();
1044          cost = model_->objective();
1045          //rowlower = model_->rowLower();
1046     }
1047#endif
1048#define TRYTHIS
1049#ifdef TRYTHIS
1050     if ((saveStrategy & 2048) != 0) {
1051          double offset;
1052          model_->getDblParam(OsiObjOffset, offset);
1053          for (i = 0; i < ncols; i++) {
1054               CoinBigIndex j;
1055               double djval = cost[i] * maxmin;
1056               for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1057                    int irow = row[j];
1058                    djval -= element[j] * lambda[irow];
1059               }
1060               cost[i] = djval;
1061          }
1062          for (i = 0; i < nrows; i++) {
1063               offset += lambda[i] * rowupper[i];
1064          }
1065          model_->setDblParam(OsiObjOffset, offset);
1066     }
1067#endif
1068     if (saveLambdaScale) {
1069          /* back off last update */
1070          for (i = 0; i < nrows; i++) {
1071               lambda[i] -= saveLambdaScale * rowsol[i];
1072          }
1073     }
1074     muAtExit_ = mu;
1075     // For last iteration make as feasible as possible
1076     if (oddSlacks)
1077          strategy_ |= 16384;
1078     // not scaled
1079     n = cleanIteration(iteration, ordStart, ordEnd,
1080                        colsol,  lower,  upper,
1081                        model_->rowLower(), model_->rowUpper(),
1082                        cost, element, fixTolerance, lastResult.objval, lastResult.infeas);
1083#if 0
1084     if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) {
1085          printf(
1086               "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
1087               iteration, mu, lastResult.infeas, lastResult.objval, n);
1088     }
1089#endif
1090#ifndef OSI_IDIOT
1091     model_->setSumPrimalInfeasibilities(lastResult.infeas);
1092#endif
1093     // Put back more feasible solution
1094     double saveInfeas[] = {0.0, 0.0};
1095     for (int iSol = 0; iSol < 3; iSol++) {
1096          const double * solution = iSol ? colsol : saveSol;
1097          if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) {
1098               // put back best solution
1099               CoinMemcpyN(saveSol, ncols, colsol);
1100          }
1101          double large = 0.0;
1102          int i;
1103          memset(rowsol, 0, nrows * sizeof(double));
1104          for (i = 0; i < ncols; i++) {
1105               CoinBigIndex j;
1106               double value = solution[i];
1107               for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1108                    int irow = row[j];
1109                    rowsol[irow] += element[j] * value;
1110               }
1111          }
1112          for (i = 0; i < nrows; i++) {
1113               if (rowsol[i] > rowupper[i]) {
1114                    double diff = rowsol[i] - rowupper[i];
1115                    if (diff > large)
1116                         large = diff;
1117               } else if (rowsol[i] < rowlower[i]) {
1118                    double diff = rowlower[i] - rowsol[i];
1119                    if (diff > large)
1120                         large = diff;
1121               }
1122          }
1123          if (iSol < 2)
1124               saveInfeas[iSol] = large;
1125          if (logLevel > 2)
1126               printf("largest infeasibility is %g\n", large);
1127     }
1128     /* subtract out lambda */
1129     for (i = 0; i < nrows; i++) {
1130          pi[i] -= lambda[i];
1131     }
1132     for (i = 0; i < ncols; i++) {
1133          CoinBigIndex j;
1134          double djval = cost[i] * maxmin;
1135          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1136               int irow = row[j];
1137               djval -= element[j] * pi[irow];
1138          }
1139          dj[i] = djval;
1140     }
1141     if ((strategy_ & 1024) != 0) {
1142          double ratio = static_cast<double> (ncols) / static_cast<double> (nrows);
1143          COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n", ratio, lastResult.infeas / firstInfeas));
1144          if (lastResult.infeas > 0.01 * firstInfeas * ratio) {
1145               strategy_ &= (~1024);
1146               COIN_DETAIL_PRINT(printf(" - layer off\n"));
1147          } else {
1148            COIN_DETAIL_PRINT(printf(" - layer on\n"));
1149          }
1150     }
1151     delete [] saveSol;
1152     delete [] lambda;
1153     // save solution
1154     // duals not much use - but save anyway
1155#ifndef OSI_IDIOT
1156     CoinMemcpyN(rowsol, nrows, model_->primalRowSolution());
1157     CoinMemcpyN(colsol, ncols, model_->primalColumnSolution());
1158     CoinMemcpyN(pi, nrows, model_->dualRowSolution());
1159     CoinMemcpyN(dj, ncols, model_->dualColumnSolution());
1160#else
1161     model_->setColSolution(colsol);
1162     model_->setRowPrice(pi);
1163     delete [] cost;
1164#endif
1165     delete [] rowsol;
1166     delete [] colsol;
1167     delete [] pi;
1168     delete [] dj;
1169     delete [] rowlower;
1170     delete [] rowupper;
1171     return ;
1172}
1173#ifndef OSI_IDIOT
1174void
1175Idiot::crossOver(int mode)
1176{
1177     if (lightWeight_ == 2) {
1178          // total failure
1179          model_->allSlackBasis();
1180          return;
1181     }
1182     double fixTolerance = IDIOT_FIX_TOLERANCE;
1183#ifdef COIN_DEVELOP
1184     double startTime = CoinCpuTime();
1185#endif
1186     ClpSimplex * saveModel = NULL;
1187     ClpMatrixBase * matrix = model_->clpMatrix();
1188     const int * row = matrix->getIndices();
1189     const CoinBigIndex * columnStart = matrix->getVectorStarts();
1190     const int * columnLength = matrix->getVectorLengths();
1191     const double * element = matrix->getElements();
1192     const double * rowupper = model_->getRowUpper();
1193     int nrows = model_->getNumRows();
1194     int ncols = model_->getNumCols();
1195     double * rowsol, * colsol;
1196     // different for Osi
1197     double * lower = model_->columnLower();
1198     double * upper = model_->columnUpper();
1199     const double * rowlower = model_->getRowLower();
1200     int * whenUsed = whenUsed_;
1201     rowsol = model_->primalRowSolution();
1202     colsol = model_->primalColumnSolution();;
1203     double * cost = model_->objective();
1204
1205     int slackEnd, ordStart, ordEnd;
1206     int slackStart = countCostedSlacks(model_);
1207
1208     int addAll = mode & 7;
1209     int presolve = 0;
1210
1211     double djTolerance = djTolerance_;
1212     if (djTolerance > 0.0 && djTolerance < 1.0)
1213          djTolerance = 1.0;
1214     int iteration;
1215     int i, n = 0;
1216     double ratio = 1.0;
1217     double objValue = 0.0;
1218     if ((strategy_ & 128) != 0) {
1219          fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
1220     }
1221     if ((mode & 16) != 0 && addAll < 3) presolve = 1;
1222     double * saveUpper = NULL;
1223     double * saveLower = NULL;
1224     double * saveRowUpper = NULL;
1225     double * saveRowLower = NULL;
1226     bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000);
1227     if (addAll < 3) {
1228          saveUpper = new double [ncols];
1229          saveLower = new double [ncols];
1230          CoinMemcpyN(upper, ncols, saveUpper);
1231          CoinMemcpyN(lower, ncols, saveLower);
1232          if (allowInfeasible) {
1233               saveRowUpper = new double [nrows];
1234               saveRowLower = new double [nrows];
1235               CoinMemcpyN(rowupper, nrows, saveRowUpper);
1236               CoinMemcpyN(rowlower, nrows, saveRowLower);
1237               double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
1238               fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas);
1239          }
1240     }
1241     if (slackStart >= 0) {
1242          slackEnd = slackStart + nrows;
1243          if (slackStart) {
1244               ordStart = 0;
1245               ordEnd = slackStart;
1246          } else {
1247               ordStart = nrows;
1248               ordEnd = ncols;
1249          }
1250     } else {
1251          slackEnd = slackStart;
1252          ordStart = 0;
1253          ordEnd = ncols;
1254     }
1255     /* get correct rowsol (without known slacks) */
1256     memset(rowsol, 0, nrows * sizeof(double));
1257     for (i = ordStart; i < ordEnd; i++) {
1258          CoinBigIndex j;
1259          double value = colsol[i];
1260          if (value < lower[i] + fixTolerance) {
1261               value = lower[i];
1262               colsol[i] = value;
1263          }
1264          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1265               int irow = row[j];
1266               rowsol[irow] += value * element[j];
1267          }
1268     }
1269     if (slackStart >= 0) {
1270          for (i = 0; i < nrows; i++) {
1271               if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) {
1272                    ratio = rowlower[i] / rowsol[i];
1273               }
1274          }
1275          for (i = 0; i < nrows; i++) {
1276               rowsol[i] *= ratio;
1277          }
1278          for (i = ordStart; i < ordEnd; i++) {
1279               double value = colsol[i] * ratio;
1280               colsol[i] = value;
1281               objValue += value * cost[i];
1282          }
1283          for (i = 0; i < nrows; i++) {
1284               double value = rowlower[i] - rowsol[i];
1285               colsol[i+slackStart] = value;
1286               objValue += value * cost[i+slackStart];
1287          }
1288          COIN_DETAIL_PRINT(printf("New objective after scaling %g\n", objValue));
1289     }
1290#if 0
1291     maybe put back - but just get feasible ?
1292     // If not many fixed then just exit
1293     int numberFixed = 0;
1294     for (i = ordStart; i < ordEnd; i++) {
1295          if (colsol[i] < lower[i] + fixTolerance)
1296               numberFixed++;
1297          else if (colsol[i] > upper[i] - fixTolerance)
1298               numberFixed++;
1299     }
1300     if (numberFixed < ncols / 2) {
1301          addAll = 3;
1302          presolve = 0;
1303     }
1304#endif
1305#ifdef FEB_TRY
1306     int savePerturbation = model_->perturbation();
1307     int saveOptions = model_->specialOptions();
1308     model_->setSpecialOptions(saveOptions | 8192);
1309     if (savePerturbation_ == 50)
1310          model_->setPerturbation(56);
1311#endif
1312     model_->createStatus();
1313     /* addAll
1314        0 - chosen,all used, all
1315        1 - chosen, all
1316        2 - all
1317        3 - do not do anything  - maybe basis
1318     */
1319     for (i = ordStart; i < ordEnd; i++) {
1320          if (addAll < 2) {
1321               if (colsol[i] < lower[i] + fixTolerance) {
1322                    upper[i] = lower[i];
1323                    colsol[i] = lower[i];
1324               } else if (colsol[i] > upper[i] - fixTolerance) {
1325                    lower[i] = upper[i];
1326                    colsol[i] = upper[i];
1327               }
1328          }
1329          model_->setColumnStatus(i, ClpSimplex::superBasic);
1330     }
1331     if ((strategy_ & 16384) != 0) {
1332          // put in basis
1333          int * posSlack = whenUsed_ + ncols;
1334          int * negSlack = posSlack + nrows;
1335          int * nextSlack = negSlack + nrows;
1336          /* Laci - try both ways - to see what works -
1337             you can change second part as much as you want */
1338#ifndef LACI_TRY  // was #if 1
1339          // Array for sorting out slack values
1340          double * ratio = new double [ncols];
1341          int * which = new int [ncols];
1342          for (i = 0; i < nrows; i++) {
1343               if (posSlack[i] >= 0 || negSlack[i] >= 0) {
1344                    int iCol;
1345                    int nPlus = 0;
1346                    int nMinus = 0;
1347                    bool possible = true;
1348                    // Get sum
1349                    double sum = 0.0;
1350                    iCol = posSlack[i];
1351                    while (iCol >= 0) {
1352                         double value = element[columnStart[iCol]];
1353                         sum += value * colsol[iCol];
1354                         if (lower[iCol]) {
1355                              possible = false;
1356                              break;
1357                         } else {
1358                              nPlus++;
1359                         }
1360                         iCol = nextSlack[iCol];
1361                    }
1362                    iCol = negSlack[i];
1363                    while (iCol >= 0) {
1364                         double value = -element[columnStart[iCol]];
1365                         sum -= value * colsol[iCol];
1366                         if (lower[iCol]) {
1367                              possible = false;
1368                              break;
1369                         } else {
1370                              nMinus++;
1371                         }
1372                         iCol = nextSlack[iCol];
1373                    }
1374                    //printf("%d plus, %d minus",nPlus,nMinus);
1375                    //printf("\n");
1376                    if ((rowsol[i] - rowlower[i] < 1.0e-7 ||
1377                              rowupper[i] - rowsol[i] < 1.0e-7) &&
1378                              nPlus + nMinus < 2)
1379                         possible = false;
1380                    if (possible) {
1381                         // Amount contributed by other varaibles
1382                         sum = rowsol[i] - sum;
1383                         double lo = rowlower[i];
1384                         if (lo > -1.0e20)
1385                              lo -= sum;
1386                         double up = rowupper[i];
1387                         if (up < 1.0e20)
1388                              up -= sum;
1389                         //printf("row bounds %g %g\n",lo,up);
1390                         if (0) {
1391                              double sum = 0.0;
1392                              double x = 0.0;
1393                              for (int k = 0; k < ncols; k++) {
1394                                   CoinBigIndex j;
1395                                   double value = colsol[k];
1396                                   x += value * cost[k];
1397                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1398                                        int irow = row[j];
1399                                        if (irow == i)
1400                                             sum += element[j] * value;
1401                                   }
1402                              }
1403                              printf("Before sum %g <= %g <= %g cost %.18g\n",
1404                                     rowlower[i], sum, rowupper[i], x);
1405                         }
1406                         // set all to zero
1407                         iCol = posSlack[i];
1408                         while (iCol >= 0) {
1409                              colsol[iCol] = 0.0;
1410                              iCol = nextSlack[iCol];
1411                         }
1412                         iCol = negSlack[i];
1413                         while (iCol >= 0) {
1414                              colsol[iCol] = 0.0;
1415                              iCol = nextSlack[iCol];
1416                         }
1417                         {
1418                              int iCol;
1419                              iCol = posSlack[i];
1420                              while (iCol >= 0) {
1421                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1422                                   //     iCol,element[columnStart[iCol]],
1423                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1424                                   iCol = nextSlack[iCol];
1425                              }
1426                              iCol = negSlack[i];
1427                              while (iCol >= 0) {
1428                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1429                                   //     iCol,element[columnStart[iCol]],
1430                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1431                                   iCol = nextSlack[iCol];
1432                              }
1433                         }
1434                         //printf("now what?\n");
1435                         int n = 0;
1436                         bool basic = false;
1437                         if (lo > 0.0) {
1438                              // Add in positive
1439                              iCol = posSlack[i];
1440                              while (iCol >= 0) {
1441                                   double value = element[columnStart[iCol]];
1442                                   ratio[n] = cost[iCol] / value;
1443                                   which[n++] = iCol;
1444                                   iCol = nextSlack[iCol];
1445                              }
1446                              CoinSort_2(ratio, ratio + n, which);
1447                              for (int i = 0; i < n; i++) {
1448                                   iCol = which[i];
1449                                   double value = element[columnStart[iCol]];
1450                                   if (lo >= upper[iCol]*value) {
1451                                        value *= upper[iCol];
1452                                        sum += value;
1453                                        lo -= value;
1454                                        colsol[iCol] = upper[iCol];
1455                                   } else {
1456                                        value = lo / value;
1457                                        sum += lo;
1458                                        lo = 0.0;
1459                                        colsol[iCol] = value;
1460                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1461                                        basic = true;
1462                                   }
1463                                   if (lo < 1.0e-7)
1464                                        break;
1465                              }
1466                         } else if (up < 0.0) {
1467                              // Use lo so coding is more similar
1468                              lo = -up;
1469                              // Add in negative
1470                              iCol = negSlack[i];
1471                              while (iCol >= 0) {
1472                                   double value = -element[columnStart[iCol]];
1473                                   ratio[n] = cost[iCol] / value;
1474                                   which[n++] = iCol;
1475                                   iCol = nextSlack[iCol];
1476                              }
1477                              CoinSort_2(ratio, ratio + n, which);
1478                              for (int i = 0; i < n; i++) {
1479                                   iCol = which[i];
1480                                   double value = -element[columnStart[iCol]];
1481                                   if (lo >= upper[iCol]*value) {
1482                                        value *= upper[iCol];
1483                                        sum += value;
1484                                        lo -= value;
1485                                        colsol[iCol] = upper[iCol];
1486                                   } else {
1487                                        value = lo / value;
1488                                        sum += lo;
1489                                        lo = 0.0;
1490                                        colsol[iCol] = value;
1491                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1492                                        basic = true;
1493                                   }
1494                                   if (lo < 1.0e-7)
1495                                        break;
1496                              }
1497                         }
1498                         if (0) {
1499                              double sum2 = 0.0;
1500                              double x = 0.0;
1501                              for (int k = 0; k < ncols; k++) {
1502                                   CoinBigIndex j;
1503                                   double value = colsol[k];
1504                                   x += value * cost[k];
1505                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1506                                        int irow = row[j];
1507                                        if (irow == i)
1508                                             sum2 += element[j] * value;
1509                                   }
1510                              }
1511                              printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n",
1512                                     rowlower[i], sum2, rowupper[i], x, sum);
1513                         }
1514                         rowsol[i] = sum;
1515                         if (basic) {
1516                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1517                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1518                              else
1519                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1520                         }
1521                    } else {
1522                         int n = 0;
1523                         int iCol;
1524                         iCol = posSlack[i];
1525                         while (iCol >= 0) {
1526                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1527                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1528                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1529                                   n++;
1530                              }
1531                              iCol = nextSlack[iCol];
1532                         }
1533                         iCol = negSlack[i];
1534                         while (iCol >= 0) {
1535                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1536                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1537                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1538                                   n++;
1539                              }
1540                              iCol = nextSlack[iCol];
1541                         }
1542                         if (n) {
1543                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1544                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1545                              else
1546                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1547#ifdef CLP_INVESTIGATE
1548                              if (n > 1)
1549                                   printf("%d basic on row %d!\n", n, i);
1550#endif
1551                         }
1552                    }
1553               }
1554          }
1555          delete [] ratio;
1556          delete [] which;
1557#else
1558          for (i = 0; i < nrows; i++) {
1559               int n = 0;
1560               int iCol;
1561               iCol = posSlack[i];
1562               while (iCol >= 0) {
1563                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1564                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1565                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1566                         n++;
1567                    }
1568                    iCol = nextSlack[iCol];
1569               }
1570               iCol = negSlack[i];
1571               while (iCol >= 0) {
1572                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1573                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1574                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1575                         n++;
1576                    }
1577                    iCol = nextSlack[iCol];
1578               }
1579               if (n) {
1580                    if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1581                         model_->setRowStatus(i, ClpSimplex::atLowerBound);
1582                    else
1583                         model_->setRowStatus(i, ClpSimplex::atUpperBound);
1584#ifdef CLP_INVESTIGATE
1585                    if (n > 1)
1586                         printf("%d basic on row %d!\n", n, i);
1587#endif
1588               }
1589          }
1590#endif
1591     }
1592     double maxmin;
1593     if (model_->getObjSense() == -1.0) {
1594          maxmin = -1.0;
1595     } else {
1596          maxmin = 1.0;
1597     }
1598     bool justValuesPass = majorIterations_ > 1000000;
1599     if (slackStart >= 0) {
1600          for (i = 0; i < nrows; i++) {
1601               model_->setRowStatus(i, ClpSimplex::superBasic);
1602          }
1603          for (i = slackStart; i < slackEnd; i++) {
1604               model_->setColumnStatus(i, ClpSimplex::basic);
1605          }
1606     } else {
1607          /* still try and put singletons rather than artificials in basis */
1608          int ninbas = 0;
1609          for (i = 0; i < nrows; i++) {
1610               model_->setRowStatus(i, ClpSimplex::basic);
1611          }
1612          for (i = 0; i < ncols; i++) {
1613               if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
1614                    CoinBigIndex j = columnStart[i];
1615                    double value = element[j];
1616                    int irow = row[j];
1617                    double rlo = rowlower[irow];
1618                    double rup = rowupper[irow];
1619                    double clo = lower[i];
1620                    double cup = upper[i];
1621                    double csol = colsol[i];
1622                    /* adjust towards feasibility */
1623                    double move = 0.0;
1624                    if (rowsol[irow] > rup) {
1625                         move = (rup - rowsol[irow]) / value;
1626                         if (value > 0.0) {
1627                              /* reduce */
1628                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1629                         } else {
1630                              /* increase */
1631                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1632                         }
1633                    } else if (rowsol[irow] < rlo) {
1634                         move = (rlo - rowsol[irow]) / value;
1635                         if (value > 0.0) {
1636                              /* increase */
1637                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1638                         } else {
1639                              /* reduce */
1640                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1641                         }
1642                    } else {
1643                         /* move to improve objective */
1644                         if (cost[i]*maxmin > 0.0) {
1645                              if (value > 0.0) {
1646                                   move = (rlo - rowsol[irow]) / value;
1647                                   /* reduce */
1648                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1649                              } else {
1650                                   move = (rup - rowsol[irow]) / value;
1651                                   /* increase */
1652                                   if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1653                              }
1654                         } else if (cost[i]*maxmin < 0.0) {
1655                              if (value > 0.0) {
1656                                   move = (rup - rowsol[irow]) / value;
1657                                   /* increase */
1658                                   if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1659                              } else {
1660                                   move = (rlo - rowsol[irow]) / value;
1661                                   /* reduce */
1662                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1663                              }
1664                         }
1665                    }
1666                    rowsol[irow] += move * value;
1667                    colsol[i] += move;
1668                    /* put in basis if row was artificial */
1669                    if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) {
1670                         model_->setRowStatus(irow, ClpSimplex::superBasic);
1671                         model_->setColumnStatus(i, ClpSimplex::basic);
1672                         ninbas++;
1673                    }
1674               }
1675          }
1676          /*printf("%d in basis\n",ninbas);*/
1677     }
1678     bool wantVector = false;
1679     if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
1680          // See if original wanted vector
1681          ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
1682          wantVector = clpMatrixO->wantsSpecialColumnCopy();
1683     }
1684     if (addAll < 3) {
1685          ClpPresolve pinfo;
1686          if (presolve) {
1687               if (allowInfeasible) {
1688                    // fix up so will be feasible
1689                    double * rhs = new double[nrows];
1690                    memset(rhs, 0, nrows * sizeof(double));
1691                    model_->clpMatrix()->times(1.0, colsol, rhs);
1692                    double * rowupper = model_->rowUpper();
1693                    double * rowlower = model_->rowLower();
1694                    saveRowUpper = CoinCopyOfArray(rowupper, nrows);
1695                    saveRowLower = CoinCopyOfArray(rowlower, nrows);
1696                    double sum = 0.0;
1697                    for (i = 0; i < nrows; i++) {
1698                         if (rhs[i] > rowupper[i]) {
1699                              sum += rhs[i] - rowupper[i];
1700                              rowupper[i] = rhs[i];
1701                         }
1702                         if (rhs[i] < rowlower[i]) {
1703                              sum += rowlower[i] - rhs[i];
1704                              rowlower[i] = rhs[i];
1705                         }
1706                    }
1707                    COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum));
1708                    delete [] rhs;
1709               }
1710               saveModel = model_;
1711               pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
1712               model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1713          }
1714          if (model_) {
1715               if (!wantVector) {
1716                    //#define TWO_GOES
1717#ifdef ABC_INHERIT
1718#ifndef TWO_GOES
1719                    model_->dealWithAbc(1,justValuesPass ? 2 : 1);
1720#else
1721                    model_->dealWithAbc(1,1 + 11);
1722#endif
1723#else
1724#ifndef TWO_GOES
1725                    model_->primal(justValuesPass ? 2 : 1);
1726#else
1727                    model_->primal(1 + 11);
1728#endif
1729#endif
1730               } else {
1731                    ClpMatrixBase * matrix = model_->clpMatrix();
1732                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1733                    assert (clpMatrix);
1734                    clpMatrix->makeSpecialColumnCopy();
1735#ifdef ABC_INHERIT
1736                    model_->dealWithAbc(1,1);
1737#else
1738                    model_->primal(1);
1739#endif
1740                    clpMatrix->releaseSpecialColumnCopy();
1741               }
1742               if (presolve) {
1743                 model_->primal();
1744                    pinfo.postsolve(true);
1745                    delete model_;
1746                    model_ = saveModel;
1747                    saveModel = NULL;
1748               }
1749          } else {
1750               // not feasible
1751               addAll = 1;
1752               presolve = 0;
1753               model_ = saveModel;
1754               saveModel = NULL;
1755               if (justValuesPass)
1756#ifdef ABC_INHERIT
1757                    model_->dealWithAbc(1,2);
1758#else
1759                    model_->primal(2);
1760#endif
1761          }
1762          if (allowInfeasible) {
1763               CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper());
1764               CoinMemcpyN(saveRowLower, nrows, model_->rowLower());
1765               delete [] saveRowUpper;
1766               delete [] saveRowLower;
1767               saveRowUpper = NULL;
1768               saveRowLower = NULL;
1769          }
1770          if (addAll < 2) {
1771               n = 0;
1772               if (!addAll ) {
1773                    /* could do scans to get a good number */
1774                    iteration = 1;
1775                    for (i = ordStart; i < ordEnd; i++) {
1776                         if (whenUsed[i] >= iteration) {
1777                              if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1778                                   n++;
1779                                   upper[i] = saveUpper[i];
1780                                   lower[i] = saveLower[i];
1781                              }
1782                         }
1783                    }
1784               } else {
1785                    for (i = ordStart; i < ordEnd; i++) {
1786                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1787                              n++;
1788                              upper[i] = saveUpper[i];
1789                              lower[i] = saveLower[i];
1790                         }
1791                    }
1792                    delete [] saveUpper;
1793                    delete [] saveLower;
1794                    saveUpper = NULL;
1795                    saveLower = NULL;
1796               }
1797#ifdef COIN_DEVELOP
1798               printf("Time so far %g, %d now added from previous iterations\n",
1799                      CoinCpuTime() - startTime, n);
1800#endif
1801               if (justValuesPass)
1802                    return;
1803               if (addAll)
1804                    presolve = 0;
1805               if (presolve) {
1806                    saveModel = model_;
1807                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1808               } else {
1809                    presolve = 0;
1810               }
1811               if (!wantVector) {
1812#ifdef ABC_INHERIT
1813                    model_->dealWithAbc(1,1);
1814#else
1815                    model_->primal(1);
1816#endif
1817               } else {
1818                    ClpMatrixBase * matrix = model_->clpMatrix();
1819                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1820                    assert (clpMatrix);
1821                    clpMatrix->makeSpecialColumnCopy();
1822#ifdef ABC_INHERIT
1823                    model_->dealWithAbc(1,1);
1824#else
1825                    model_->primal(1);
1826#endif
1827                    clpMatrix->releaseSpecialColumnCopy();
1828               }
1829               if (presolve) {
1830                    pinfo.postsolve(true);
1831                    delete model_;
1832                    model_ = saveModel;
1833                    saveModel = NULL;
1834               }
1835               if (!addAll) {
1836                    n = 0;
1837                    for (i = ordStart; i < ordEnd; i++) {
1838                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1839                              n++;
1840                              upper[i] = saveUpper[i];
1841                              lower[i] = saveLower[i];
1842                         }
1843                    }
1844                    delete [] saveUpper;
1845                    delete [] saveLower;
1846                    saveUpper = NULL;
1847                    saveLower = NULL;
1848#ifdef COIN_DEVELOP
1849                    printf("Time so far %g, %d now added from previous iterations\n",
1850                           CoinCpuTime() - startTime, n);
1851#endif
1852               }
1853               if (presolve) {
1854                    saveModel = model_;
1855                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1856               } else {
1857                    presolve = 0;
1858               }
1859               if (!wantVector) {
1860#ifdef ABC_INHERIT
1861                    model_->dealWithAbc(1,1);
1862#else
1863                    model_->primal(1);
1864#endif
1865               } else {
1866                    ClpMatrixBase * matrix = model_->clpMatrix();
1867                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1868                    assert (clpMatrix);
1869                    clpMatrix->makeSpecialColumnCopy();
1870#ifdef ABC_INHERIT
1871                    model_->dealWithAbc(1,1);
1872#else
1873                    model_->primal(1);
1874#endif
1875                    clpMatrix->releaseSpecialColumnCopy();
1876               }
1877               if (presolve) {
1878                    pinfo.postsolve(true);
1879                    delete model_;
1880                    model_ = saveModel;
1881                    saveModel = NULL;
1882               }
1883          }
1884#ifdef COIN_DEVELOP
1885          printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
1886#endif
1887          delete [] saveUpper;
1888          delete [] saveLower;
1889     }
1890#ifdef FEB_TRY
1891     model_->setSpecialOptions(saveOptions);
1892     model_->setPerturbation(savePerturbation);
1893#endif
1894     return ;
1895}
1896#endif
1897/*****************************************************************************/
1898
1899// Default contructor
1900Idiot::Idiot()
1901{
1902     model_ = NULL;
1903     maxBigIts_ = 3;
1904     maxIts_ = 5;
1905     logLevel_ = 1;
1906     logFreq_ = 100;
1907     maxIts2_ = 100;
1908     djTolerance_ = 1e-1;
1909     mu_ = 1e-4;
1910     drop_ = 5.0;
1911     exitDrop_ = -1.0e20;
1912     muFactor_ = 0.3333;
1913     stopMu_ = 1e-12;
1914     smallInfeas_ = 1e-1;
1915     reasonableInfeas_ = 1e2;
1916     muAtExit_ = 1.0e31;
1917     strategy_ = 8;
1918     lambdaIterations_ = 0;
1919     checkFrequency_ = 100;
1920     whenUsed_ = NULL;
1921     majorIterations_ = 30;
1922     exitFeasibility_ = -1.0;
1923     dropEnoughFeasibility_ = 0.02;
1924     dropEnoughWeighted_ = 0.01;
1925     // adjust
1926     double nrows = 10000.0;
1927     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
1928     baseIts = baseIts / 10;
1929     baseIts *= 10;
1930     maxIts2_ = 200 + baseIts + 5;
1931     maxIts2_ = 100;
1932     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
1933     lightWeight_ = 0;
1934}
1935// Constructor from model
1936Idiot::Idiot(OsiSolverInterface &model)
1937{
1938     model_ = & model;
1939     maxBigIts_ = 3;
1940     maxIts_ = 5;
1941     logLevel_ = 1;
1942     logFreq_ = 100;
1943     maxIts2_ = 100;
1944     djTolerance_ = 1e-1;
1945     mu_ = 1e-4;
1946     drop_ = 5.0;
1947     exitDrop_ = -1.0e20;
1948     muFactor_ = 0.3333;
1949     stopMu_ = 1e-12;
1950     smallInfeas_ = 1e-1;
1951     reasonableInfeas_ = 1e2;
1952     muAtExit_ = 1.0e31;
1953     strategy_ = 8;
1954     lambdaIterations_ = 0;
1955     checkFrequency_ = 100;
1956     whenUsed_ = NULL;
1957     majorIterations_ = 30;
1958     exitFeasibility_ = -1.0;
1959     dropEnoughFeasibility_ = 0.02;
1960     dropEnoughWeighted_ = 0.01;
1961     // adjust
1962     double nrows;
1963     if (model_)
1964          nrows = model_->getNumRows();
1965     else
1966          nrows = 10000.0;
1967     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
1968     baseIts = baseIts / 10;
1969     baseIts *= 10;
1970     maxIts2_ = 200 + baseIts + 5;
1971     maxIts2_ = 100;
1972     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
1973     lightWeight_ = 0;
1974}
1975// Copy constructor.
1976Idiot::Idiot(const Idiot &rhs)
1977{
1978     model_ = rhs.model_;
1979     if (model_ && rhs.whenUsed_) {
1980          int numberColumns = model_->getNumCols();
1981          whenUsed_ = new int [numberColumns];
1982          CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
1983     } else {
1984          whenUsed_ = NULL;
1985     }
1986     djTolerance_ = rhs.djTolerance_;
1987     mu_ = rhs.mu_;
1988     drop_ = rhs.drop_;
1989     muFactor_ = rhs.muFactor_;
1990     stopMu_ = rhs.stopMu_;
1991     smallInfeas_ = rhs.smallInfeas_;
1992     reasonableInfeas_ = rhs.reasonableInfeas_;
1993     exitDrop_ = rhs.exitDrop_;
1994     muAtExit_ = rhs.muAtExit_;
1995     exitFeasibility_ = rhs.exitFeasibility_;
1996     dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1997     dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1998     maxBigIts_ = rhs.maxBigIts_;
1999     maxIts_ = rhs.maxIts_;
2000     majorIterations_ = rhs.majorIterations_;
2001     logLevel_ = rhs.logLevel_;
2002     logFreq_ = rhs.logFreq_;
2003     checkFrequency_ = rhs.checkFrequency_;
2004     lambdaIterations_ = rhs.lambdaIterations_;
2005     maxIts2_ = rhs.maxIts2_;
2006     strategy_ = rhs.strategy_;
2007     lightWeight_ = rhs.lightWeight_;
2008}
2009// Assignment operator. This copies the data
2010Idiot &
2011Idiot::operator=(const Idiot & rhs)
2012{
2013     if (this != &rhs) {
2014          delete [] whenUsed_;
2015          model_ = rhs.model_;
2016          if (model_ && rhs.whenUsed_) {
2017               int numberColumns = model_->getNumCols();
2018               whenUsed_ = new int [numberColumns];
2019               CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
2020          } else {
2021               whenUsed_ = NULL;
2022          }
2023          djTolerance_ = rhs.djTolerance_;
2024          mu_ = rhs.mu_;
2025          drop_ = rhs.drop_;
2026          muFactor_ = rhs.muFactor_;
2027          stopMu_ = rhs.stopMu_;
2028          smallInfeas_ = rhs.smallInfeas_;
2029          reasonableInfeas_ = rhs.reasonableInfeas_;
2030          exitDrop_ = rhs.exitDrop_;
2031          muAtExit_ = rhs.muAtExit_;
2032          exitFeasibility_ = rhs.exitFeasibility_;
2033          dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
2034          dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
2035          maxBigIts_ = rhs.maxBigIts_;
2036          maxIts_ = rhs.maxIts_;
2037          majorIterations_ = rhs.majorIterations_;
2038          logLevel_ = rhs.logLevel_;
2039          logFreq_ = rhs.logFreq_;
2040          checkFrequency_ = rhs.checkFrequency_;
2041          lambdaIterations_ = rhs.lambdaIterations_;
2042          maxIts2_ = rhs.maxIts2_;
2043          strategy_ = rhs.strategy_;
2044          lightWeight_ = rhs.lightWeight_;
2045     }
2046     return *this;
2047}
2048Idiot::~Idiot()
2049{
2050     delete [] whenUsed_;
2051}
Note: See TracBrowser for help on using the repository browser.