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

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 80.2 KB
Line 
1/* $Id: Idiot.cpp 2030 2014-04-15 15:54:11Z forrest $ */
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, 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          if (slackStart) {
620               ordStart = 0;
621               ordEnd = slackStart;
622          } else {
623               ordStart = nrows;
624               ordEnd = ncols;
625          }
626     } else {
627          ordStart = 0;
628          ordEnd = ncols;
629     }
630     if (offset && logLevel > 2) {
631          printf("** Objective offset is %g\n", offset);
632     }
633     /* compute reasonable solution cost */
634     for (i = 0; i < nrows; i++) {
635          rowsol[i] = 1.0e31;
636     }
637     for (i = 0; i < ncols; i++) {
638          CoinBigIndex j;
639          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
640               if (element[j] != 1.0) {
641                    allOnes = 0;
642                    break;
643               }
644          }
645     }
646     if (allOnes) {
647          elemXX = NULL;
648     } else {
649          elemXX = element;
650     }
651     // Do scaling if wanted
652     bool scaled = false;
653#ifndef OSI_IDIOT
654     if ((strategy_ & 32) != 0 && !allOnes) {
655          if (model_->scalingFlag() > 0)
656               scaled = model_->clpMatrix()->scale(model_) == 0;
657          if (scaled) {
658               const double * rowScale = model_->rowScale();
659               const double * columnScale = model_->columnScale();
660               double * oldLower = lower;
661               double * oldUpper = upper;
662               double * oldCost = cost;
663               lower = new double[ncols];
664               upper = new double[ncols];
665               cost = new double[ncols];
666               CoinMemcpyN(oldLower, ncols, lower);
667               CoinMemcpyN(oldUpper, ncols, upper);
668               CoinMemcpyN(oldCost, ncols, cost);
669               int icol, irow;
670               for (icol = 0; icol < ncols; icol++) {
671                    double multiplier = 1.0 / columnScale[icol];
672                    if (lower[icol] > -1.0e50)
673                         lower[icol] *= multiplier;
674                    if (upper[icol] < 1.0e50)
675                         upper[icol] *= multiplier;
676                    colsol[icol] *= multiplier;
677                    cost[icol] *= columnScale[icol];
678               }
679               CoinMemcpyN(model_->rowLower(), nrows, rowlower);
680               for (irow = 0; irow < nrows; irow++) {
681                    double multiplier = rowScale[irow];
682                    if (rowlower[irow] > -1.0e50)
683                         rowlower[irow] *= multiplier;
684                    if (rowupper[irow] < 1.0e50)
685                         rowupper[irow] *= multiplier;
686                    rowsol[irow] *= multiplier;
687               }
688               int length = columnStart[ncols-1] + columnLength[ncols-1];
689               double * elemYY = new double[length];
690               for (i = 0; i < ncols; i++) {
691                    CoinBigIndex j;
692                    double scale = columnScale[i];
693                    for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
694                         int irow = row[j];
695                         elemYY[j] = element[j] * scale * rowScale[irow];
696                    }
697               }
698               elemXX = elemYY;
699          }
700     }
701#endif
702     for (i = 0; i < ncols; i++) {
703          CoinBigIndex j;
704          double dd = columnLength[i];
705          dd = cost[i] / dd;
706          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
707               int irow = row[j];
708               if (dd < rowsol[irow]) {
709                    rowsol[irow] = dd;
710               }
711          }
712     }
713     d2 = 0.0;
714     for (i = 0; i < nrows; i++) {
715          d2 += rowsol[i];
716     }
717     d2 *= 2.0; /* for luck */
718
719     d2 = d2 / static_cast<double> (4 * nrows + 8000);
720     d2 *= 0.5; /* halve with more flexible method */
721     if (d2 < 5.0) d2 = 5.0;
722     if (djExit == 0.0) {
723          djExit = d2;
724     }
725     if ((saveStrategy & 4) != 0) {
726          /* go to relative tolerances - first small */
727          djExit = 1.0e-10;
728          djFlag = 1.0e-5;
729          drop = 1.0e-10;
730     }
731     memset(whenUsed, 0, ncols * sizeof(int));
732     strategy = strategies[strategy];
733     if ((saveStrategy & 8) != 0) strategy |= 64; /* don't allow large theta */
734     CoinMemcpyN(colsol, ncols, saveSol);
735
736     lastResult = IdiSolve(nrows, ncols, rowsol , colsol, pi,
737                           dj, cost, rowlower, rowupper,
738                           lower, upper, elemXX, row, columnStart, columnLength, lambda,
739                           0, mu, drop,
740                           maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
741     // update whenUsed_
742     n = cleanIteration(iteration, ordStart, ordEnd,
743                        colsol,  lower,  upper,
744                        rowlower, rowupper,
745                        cost, elemXX, fixTolerance, lastResult.objval, lastResult.infeas);
746     if ((strategy_ & 16384) != 0) {
747          int * posSlack = whenUsed_ + ncols;
748          int * negSlack = posSlack + nrows;
749          int * nextSlack = negSlack + nrows;
750          double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
751          for (i = 0; i < nrows; i++)
752               rowsol[i] += rowsol2[i];
753     }
754     if ((logLevel_ & 1) != 0) {
755#ifndef OSI_IDIOT
756          if (!handler) {
757#endif
758               printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
759                      iteration, lastResult.infeas, lastResult.objval, mu, lastResult.iteration, n);
760#ifndef OSI_IDIOT
761          } else {
762               handler->message(CLP_IDIOT_ITERATION, *messages)
763                         << iteration << lastResult.infeas << lastResult.objval << mu << lastResult.iteration << n
764                         << CoinMessageEol;
765          }
766#endif
767     }
768     int numberBaseTrys = 0; // for first time
769     int numberAway = -1;
770     iterationTotal = lastResult.iteration;
771     firstInfeas = lastResult.infeas;
772     if ((strategy_ & 1024) != 0) reasonableInfeas = 0.5 * firstInfeas;
773     if (lastResult.infeas < reasonableInfeas) lastResult.infeas = reasonableInfeas;
774     double keepinfeas = 1.0e31;
775     double lastInfeas = 1.0e31;
776     double bestInfeas = 1.0e31;
777     while ((mu > stopMu && lastResult.infeas > smallInfeas) ||
778               (lastResult.infeas <= smallInfeas &&
779                dropping(lastResult, exitDrop, smallInfeas, &badIts)) ||
780               checkIteration < 2 || lambdaIteration < lambdaIterations_) {
781          if (lastResult.infeas <= exitFeasibility_)
782               break;
783          iteration++;
784          checkIteration++;
785          if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
786               bestFeasible = lastResult.objval;
787          }
788          if (lastResult.infeas + mu * lastResult.objval < bestWeighted) {
789               bestWeighted = lastResult.objval + mu * lastResult.objval;
790          }
791          if ((saveStrategy & 4096)) strategy |= 256;
792          if ((saveStrategy & 4) != 0 && iteration > 2) {
793               /* go to relative tolerances */
794               double weighted = 10.0 + fabs(lastWeighted);
795               djExit = djTolerance_ * weighted;
796               djFlag = 2.0 * djExit;
797               drop = drop_ * weighted;
798               djTol = 0.01 * djExit;
799          }
800          CoinMemcpyN(colsol, ncols, saveSol);
801          result = IdiSolve(nrows, ncols, rowsol , colsol, pi, dj,
802                            cost, rowlower, rowupper,
803                            lower, upper, elemXX, row, columnStart, columnLength, lambda,
804                            maxIts, mu, drop,
805                            maxmin, offset, strategy, djTol, djExit, djFlag, randomNumberGenerator);
806          n = cleanIteration(iteration, ordStart, ordEnd,
807                             colsol,  lower,  upper,
808                             rowlower, rowupper,
809                             cost, elemXX, fixTolerance, result.objval, result.infeas);
810          if ((strategy_ & 16384) != 0) {
811               int * posSlack = whenUsed_ + ncols;
812               int * negSlack = posSlack + nrows;
813               int * nextSlack = negSlack + nrows;
814               double * rowsol2 = reinterpret_cast<double *> (nextSlack + ncols);
815               for (i = 0; i < nrows; i++)
816                    rowsol[i] += rowsol2[i];
817          }
818          if ((logLevel_ & 1) != 0) {
819#ifndef OSI_IDIOT
820               if (!handler) {
821#endif
822                    printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
823                           iteration, result.infeas, result.objval, mu, result.iteration, n);
824#ifndef OSI_IDIOT
825               } else {
826                    handler->message(CLP_IDIOT_ITERATION, *messages)
827                              << iteration << result.infeas << result.objval << mu << result.iteration << n
828                              << CoinMessageEol;
829               }
830#endif
831          }
832          if (iteration > 50 && n == numberAway ) {
833            if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) {
834#ifdef CLP_INVESTIGATE
835               printf("infeas small %g\n", result.infeas);
836#endif
837               break; // not much happening
838            }
839          }
840          if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
841               if (lastInfeas != bestInfeas && CoinMin(result.infeas, lastInfeas) > 0.95 * bestInfeas)
842                    majorIterations_ = CoinMin(majorIterations_, iteration); // not getting feasible
843          }
844          lastInfeas = result.infeas;
845          numberAway = n;
846          keepinfeas = result.infeas;
847          lastWeighted = result.weighted;
848          iterationTotal += result.iteration;
849          if (iteration == 1) {
850               if ((strategy_ & 1024) != 0 && mu < 1.0e-10)
851                    result.infeas = firstInfeas * 0.8;
852               if (majorIterations_ >= 50 || dropEnoughFeasibility_ <= 0.0)
853                    result.infeas *= 0.8;
854               if (result.infeas > firstInfeas * 0.9
855                         && result.infeas > reasonableInfeas) {
856                    iteration--;
857                    if (majorIterations_ < 50)
858                         mu *= 1.0e-1;
859                    else
860                         mu *= 0.7;
861                    bestFeasible = 1.0e31;
862                    bestWeighted = 1.0e60;
863                    numberBaseTrys++;
864                    if (mu < 1.0e-30 || (numberBaseTrys > 10 && lightWeight_)) {
865                         // back to all slack basis
866                         lightWeight_ = 2;
867                         break;
868                    }
869                    CoinMemcpyN(saveSol, ncols, colsol);
870               } else {
871                    maxIts = maxIts2;
872                    checkIteration = 0;
873                    if ((strategy_ & 1024) != 0) mu *= 1.0e-1;
874               }
875          } else {
876          }
877          bestInfeas = CoinMin(bestInfeas, result.infeas);
878          if (majorIterations_>100&&majorIterations_<200) {
879            if (iteration==majorIterations_-100) {
880              // redo
881              double muX=mu*10.0;
882              bestInfeas=1.0e3;
883              mu=muX;
884              nTry=0;
885            }
886          }
887          if (iteration) {
888               /* this code is in to force it to terminate sometime */
889               double changeMu = factor;
890               if ((saveStrategy & 64) != 0) {
891                    keepinfeas = 0.0; /* switch off ranga's increase */
892                    fakeSmall = smallInfeas;
893               } else {
894                    fakeSmall = -1.0;
895               }
896               saveLambdaScale = 0.0;
897               if (result.infeas > reasonableInfeas ||
898                         (nTry + 1 == maxBigIts && result.infeas > fakeSmall)) {
899                    if (result.infeas > lastResult.infeas*(1.0 - dropEnoughFeasibility_) ||
900                              nTry + 1 == maxBigIts ||
901                              (result.infeas > lastResult.infeas * 0.9
902                               && result.weighted > lastResult.weighted
903                               - dropEnoughWeighted_ * CoinMax(fabs(lastResult.weighted), fabs(result.weighted)))) {
904                         mu *= changeMu;
905                         if ((saveStrategy & 32) != 0 && result.infeas < reasonableInfeas && 0) {
906                              reasonableInfeas = CoinMax(smallInfeas, reasonableInfeas * sqrt(changeMu));
907                              COIN_DETAIL_PRINT(printf("reasonable infeas now %g\n", reasonableInfeas));
908                         }
909                         result.weighted = 1.0e60;
910                         nTry = 0;
911                         bestFeasible = 1.0e31;
912                         bestWeighted = 1.0e60;
913                         checkIteration = 0;
914                         lambdaIteration = 0;
915#define LAMBDA
916#ifdef LAMBDA
917                         if ((saveStrategy & 2048) == 0) {
918                              memset(lambda, 0, nrows * sizeof(double));
919                         }
920#else
921                         memset(lambda, 0, nrows * sizeof(double));
922#endif
923                    } else {
924                         nTry++;
925                    }
926               } else if (lambdaIterations_ >= 0) {
927                    /* update lambda  */
928                    double scale = 1.0 / mu;
929                    int i, nnz = 0;
930                    saveLambdaScale = scale;
931                    lambdaIteration++;
932                    if ((saveStrategy & 4) == 0) drop = drop_ / 50.0;
933                    if (lambdaIteration > 4 &&
934                              (((lambdaIteration % 10) == 0 && smallInfeas < keepinfeas) ||
935                               ((lambdaIteration % 5) == 0 && 1.5 * smallInfeas < keepinfeas))) {
936                         //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
937                         smallInfeas *= 1.5;
938                    }
939                    if ((saveStrategy & 2048) == 0) {
940                         for (i = 0; i < nrows; i++) {
941                              if (lambda[i]) nnz++;
942                              lambda[i] += scale * rowsol[i];
943                         }
944                    } else {
945                         nnz = 1;
946#ifdef LAMBDA
947                         for (i = 0; i < nrows; i++) {
948                              lambda[i] += scale * rowsol[i];
949                         }
950#else
951                         for (i = 0; i < nrows; i++) {
952                              lambda[i] = scale * rowsol[i];
953                         }
954                         for (i = 0; i < ncols; i++) {
955                              CoinBigIndex j;
956                              double value = cost[i] * maxmin;
957                              for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
958                                   int irow = row[j];
959                                   value += element[j] * lambda[irow];
960                              }
961                              cost[i] = value * maxmin;
962                         }
963                         for (i = 0; i < nrows; i++) {
964                              offset += lambda[i] * rowupper[i];
965                              lambda[i] = 0.0;
966                         }
967#ifdef DEBUG
968                         printf("offset %g\n", offset);
969#endif
970                         model_->setDblParam(OsiObjOffset, offset);
971#endif
972                    }
973                    nTry++;
974                    if (!nnz) {
975                         bestFeasible = 1.0e32;
976                         bestWeighted = 1.0e60;
977                         checkIteration = 0;
978                         result.weighted = 1.0e31;
979                    }
980#ifdef DEBUG
981                    double trueCost = 0.0;
982                    for (i = 0; i < ncols; i++) {
983                         int j;
984                         trueCost += cost[i] * colsol[i];
985                    }
986                    printf("True objective %g\n", trueCost - offset);
987#endif
988               } else {
989                    nTry++;
990               }
991               lastResult = result;
992               if (result.infeas < reasonableInfeas && !belowReasonable) {
993                    belowReasonable = 1;
994                    bestFeasible = 1.0e32;
995                    bestWeighted = 1.0e60;
996                    checkIteration = 0;
997                    result.weighted = 1.0e31;
998               }
999          }
1000          if (iteration >= majorIterations_) {
1001               // If not feasible and crash then dive dive dive
1002               if (mu > 1.0e-12 && result.infeas > 1.0 && majorIterations_ < 40) {
1003                    mu = 1.0e-30;
1004                    majorIterations_ = iteration + 1;
1005                    stopMu = 0.0;
1006               } else {
1007                    if (logLevel > 2)
1008                         printf("Exiting due to number of major iterations\n");
1009                    break;
1010               }
1011          }
1012     }
1013     majorIterations_ = saveMajorIterations;
1014#ifndef OSI_IDIOT
1015     if (scaled) {
1016          // Scale solution and free arrays
1017          const double * rowScale = model_->rowScale();
1018          const double * columnScale = model_->columnScale();
1019          int icol, irow;
1020          for (icol = 0; icol < ncols; icol++) {
1021               colsol[icol] *= columnScale[icol];
1022               saveSol[icol] *= columnScale[icol];
1023               dj[icol] /= columnScale[icol];
1024          }
1025          for (irow = 0; irow < nrows; irow++) {
1026               rowsol[irow] /= rowScale[irow];
1027               pi[irow] *= rowScale[irow];
1028          }
1029          // Don't know why getting Microsoft problems
1030#if defined (_MSC_VER)
1031          delete [] ( double *) elemXX;
1032#else
1033          delete [] elemXX;
1034#endif
1035          model_->setRowScale(NULL);
1036          model_->setColumnScale(NULL);
1037          delete [] lower;
1038          delete [] upper;
1039          delete [] cost;
1040          lower = model_->columnLower();
1041          upper = model_->columnUpper();
1042          cost = model_->objective();
1043          //rowlower = model_->rowLower();
1044     }
1045#endif
1046#define TRYTHIS
1047#ifdef TRYTHIS
1048     if ((saveStrategy & 2048) != 0) {
1049          double offset;
1050          model_->getDblParam(OsiObjOffset, offset);
1051          for (i = 0; i < ncols; i++) {
1052               CoinBigIndex j;
1053               double djval = cost[i] * maxmin;
1054               for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1055                    int irow = row[j];
1056                    djval -= element[j] * lambda[irow];
1057               }
1058               cost[i] = djval;
1059          }
1060          for (i = 0; i < nrows; i++) {
1061               offset += lambda[i] * rowupper[i];
1062          }
1063          model_->setDblParam(OsiObjOffset, offset);
1064     }
1065#endif
1066     if (saveLambdaScale) {
1067          /* back off last update */
1068          for (i = 0; i < nrows; i++) {
1069               lambda[i] -= saveLambdaScale * rowsol[i];
1070          }
1071     }
1072     muAtExit_ = mu;
1073     // For last iteration make as feasible as possible
1074     if (oddSlacks && (strategy_&32768)==0)
1075          strategy_ |= 16384;
1076     // not scaled
1077     n = cleanIteration(iteration, ordStart, ordEnd,
1078                        colsol,  lower,  upper,
1079                        model_->rowLower(), model_->rowUpper(),
1080                        cost, element, fixTolerance, lastResult.objval, lastResult.infeas);
1081#if 0
1082     if ((logLevel & 1) == 0 || (strategy_ & 16384) != 0) {
1083          printf(
1084               "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
1085               iteration, mu, lastResult.infeas, lastResult.objval, n);
1086     }
1087#endif
1088#ifndef OSI_IDIOT
1089     model_->setSumPrimalInfeasibilities(lastResult.infeas);
1090#endif
1091     // Put back more feasible solution
1092     double saveInfeas[] = {0.0, 0.0};
1093     for (int iSol = 0; iSol < 3; iSol++) {
1094          const double * solution = iSol ? colsol : saveSol;
1095          if (iSol == 2 && saveInfeas[0] < saveInfeas[1]) {
1096               // put back best solution
1097               CoinMemcpyN(saveSol, ncols, colsol);
1098          }
1099          double large = 0.0;
1100          int i;
1101          memset(rowsol, 0, nrows * sizeof(double));
1102          for (i = 0; i < ncols; i++) {
1103               CoinBigIndex j;
1104               double value = solution[i];
1105               for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1106                    int irow = row[j];
1107                    rowsol[irow] += element[j] * value;
1108               }
1109          }
1110          for (i = 0; i < nrows; i++) {
1111               if (rowsol[i] > rowupper[i]) {
1112                    double diff = rowsol[i] - rowupper[i];
1113                    if (diff > large)
1114                         large = diff;
1115               } else if (rowsol[i] < rowlower[i]) {
1116                    double diff = rowlower[i] - rowsol[i];
1117                    if (diff > large)
1118                         large = diff;
1119               }
1120          }
1121          if (iSol < 2)
1122               saveInfeas[iSol] = large;
1123          if (logLevel > 2)
1124               printf("largest infeasibility is %g\n", large);
1125     }
1126     /* subtract out lambda */
1127     for (i = 0; i < nrows; i++) {
1128          pi[i] -= lambda[i];
1129     }
1130     for (i = 0; i < ncols; i++) {
1131          CoinBigIndex j;
1132          double djval = cost[i] * maxmin;
1133          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1134               int irow = row[j];
1135               djval -= element[j] * pi[irow];
1136          }
1137          dj[i] = djval;
1138     }
1139     if ((strategy_ & 1024) != 0) {
1140          double ratio = static_cast<double> (ncols) / static_cast<double> (nrows);
1141          COIN_DETAIL_PRINT(printf("col/row ratio %g infeas ratio %g\n", ratio, lastResult.infeas / firstInfeas));
1142          if (lastResult.infeas > 0.01 * firstInfeas * ratio) {
1143               strategy_ &= (~1024);
1144               COIN_DETAIL_PRINT(printf(" - layer off\n"));
1145          } else {
1146            COIN_DETAIL_PRINT(printf(" - layer on\n"));
1147          }
1148     }
1149     delete [] saveSol;
1150     delete [] lambda;
1151     // save solution
1152     // duals not much use - but save anyway
1153#ifndef OSI_IDIOT
1154     CoinMemcpyN(rowsol, nrows, model_->primalRowSolution());
1155     CoinMemcpyN(colsol, ncols, model_->primalColumnSolution());
1156     CoinMemcpyN(pi, nrows, model_->dualRowSolution());
1157     CoinMemcpyN(dj, ncols, model_->dualColumnSolution());
1158#else
1159     model_->setColSolution(colsol);
1160     model_->setRowPrice(pi);
1161     delete [] cost;
1162#endif
1163     delete [] rowsol;
1164     delete [] colsol;
1165     delete [] pi;
1166     delete [] dj;
1167     delete [] rowlower;
1168     delete [] rowupper;
1169     return ;
1170}
1171#ifndef OSI_IDIOT
1172void
1173Idiot::crossOver(int mode)
1174{
1175     if (lightWeight_ == 2) {
1176          // total failure
1177          model_->allSlackBasis();
1178          return;
1179     }
1180     double fixTolerance = IDIOT_FIX_TOLERANCE;
1181#ifdef COIN_DEVELOP
1182     double startTime = CoinCpuTime();
1183#endif
1184     ClpSimplex * saveModel = NULL;
1185     ClpMatrixBase * matrix = model_->clpMatrix();
1186     const int * row = matrix->getIndices();
1187     const CoinBigIndex * columnStart = matrix->getVectorStarts();
1188     const int * columnLength = matrix->getVectorLengths();
1189     const double * element = matrix->getElements();
1190     const double * rowupper = model_->getRowUpper();
1191     int nrows = model_->getNumRows();
1192     int ncols = model_->getNumCols();
1193     double * rowsol, * colsol;
1194     // different for Osi
1195     double * lower = model_->columnLower();
1196     double * upper = model_->columnUpper();
1197     const double * rowlower = model_->getRowLower();
1198     int * whenUsed = whenUsed_;
1199     rowsol = model_->primalRowSolution();
1200     colsol = model_->primalColumnSolution();;
1201     double * cost = model_->objective();
1202
1203     int slackEnd, ordStart, ordEnd;
1204     int slackStart = countCostedSlacks(model_);
1205
1206     int addAll = mode & 7;
1207     int presolve = 0;
1208
1209     double djTolerance = djTolerance_;
1210     if (djTolerance > 0.0 && djTolerance < 1.0)
1211          djTolerance = 1.0;
1212     int iteration;
1213     int i, n = 0;
1214     double ratio = 1.0;
1215     double objValue = 0.0;
1216     if ((strategy_ & 128) != 0) {
1217          fixTolerance = SMALL_IDIOT_FIX_TOLERANCE;
1218     }
1219     if ((mode & 16) != 0 && addAll < 3) presolve = 1;
1220     double * saveUpper = NULL;
1221     double * saveLower = NULL;
1222     double * saveRowUpper = NULL;
1223     double * saveRowLower = NULL;
1224     bool allowInfeasible = ((strategy_ & 8192) != 0) || (majorIterations_ > 1000000);
1225     if (addAll < 3) {
1226          saveUpper = new double [ncols];
1227          saveLower = new double [ncols];
1228          CoinMemcpyN(upper, ncols, saveUpper);
1229          CoinMemcpyN(lower, ncols, saveLower);
1230          if (allowInfeasible) {
1231               saveRowUpper = new double [nrows];
1232               saveRowLower = new double [nrows];
1233               CoinMemcpyN(rowupper, nrows, saveRowUpper);
1234               CoinMemcpyN(rowlower, nrows, saveRowLower);
1235               double averageInfeas = model_->sumPrimalInfeasibilities() / static_cast<double> (model_->numberRows());
1236               fixTolerance = CoinMax(fixTolerance, 1.0e-5 * averageInfeas);
1237          }
1238     }
1239     if (slackStart >= 0) {
1240          slackEnd = slackStart + nrows;
1241          if (slackStart) {
1242               ordStart = 0;
1243               ordEnd = slackStart;
1244          } else {
1245               ordStart = nrows;
1246               ordEnd = ncols;
1247          }
1248     } else {
1249          slackEnd = slackStart;
1250          ordStart = 0;
1251          ordEnd = ncols;
1252     }
1253     /* get correct rowsol (without known slacks) */
1254     memset(rowsol, 0, nrows * sizeof(double));
1255     for (i = ordStart; i < ordEnd; i++) {
1256          CoinBigIndex j;
1257          double value = colsol[i];
1258          if (value < lower[i] + fixTolerance) {
1259               value = lower[i];
1260               colsol[i] = value;
1261          }
1262          for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
1263               int irow = row[j];
1264               rowsol[irow] += value * element[j];
1265          }
1266     }
1267     if (slackStart >= 0) {
1268          for (i = 0; i < nrows; i++) {
1269               if (ratio * rowsol[i] > rowlower[i] && rowsol[i] > 1.0e-8) {
1270                    ratio = rowlower[i] / rowsol[i];
1271               }
1272          }
1273          for (i = 0; i < nrows; i++) {
1274               rowsol[i] *= ratio;
1275          }
1276          for (i = ordStart; i < ordEnd; i++) {
1277               double value = colsol[i] * ratio;
1278               colsol[i] = value;
1279               objValue += value * cost[i];
1280          }
1281          for (i = 0; i < nrows; i++) {
1282               double value = rowlower[i] - rowsol[i];
1283               colsol[i+slackStart] = value;
1284               objValue += value * cost[i+slackStart];
1285          }
1286          COIN_DETAIL_PRINT(printf("New objective after scaling %g\n", objValue));
1287     }
1288#if 0
1289     maybe put back - but just get feasible ?
1290     // If not many fixed then just exit
1291     int numberFixed = 0;
1292     for (i = ordStart; i < ordEnd; i++) {
1293          if (colsol[i] < lower[i] + fixTolerance)
1294               numberFixed++;
1295          else if (colsol[i] > upper[i] - fixTolerance)
1296               numberFixed++;
1297     }
1298     if (numberFixed < ncols / 2) {
1299          addAll = 3;
1300          presolve = 0;
1301     }
1302#endif
1303#ifdef FEB_TRY
1304     int savePerturbation = model_->perturbation();
1305     int saveOptions = model_->specialOptions();
1306     model_->setSpecialOptions(saveOptions | 8192);
1307     if (savePerturbation_ == 50)
1308          model_->setPerturbation(56);
1309#endif
1310     model_->createStatus();
1311     /* addAll
1312        0 - chosen,all used, all
1313        1 - chosen, all
1314        2 - all
1315        3 - do not do anything  - maybe basis
1316     */
1317     for (i = ordStart; i < ordEnd; i++) {
1318          if (addAll < 2) {
1319               if (colsol[i] < lower[i] + fixTolerance) {
1320                    upper[i] = lower[i];
1321                    colsol[i] = lower[i];
1322               } else if (colsol[i] > upper[i] - fixTolerance) {
1323                    lower[i] = upper[i];
1324                    colsol[i] = upper[i];
1325               }
1326          }
1327          model_->setColumnStatus(i, ClpSimplex::superBasic);
1328     }
1329     if ((strategy_ & 16384) != 0) {
1330          // put in basis
1331          int * posSlack = whenUsed_ + ncols;
1332          int * negSlack = posSlack + nrows;
1333          int * nextSlack = negSlack + nrows;
1334          /* Laci - try both ways - to see what works -
1335             you can change second part as much as you want */
1336#ifndef LACI_TRY  // was #if 1
1337          // Array for sorting out slack values
1338          double * ratio = new double [ncols];
1339          int * which = new int [ncols];
1340          for (i = 0; i < nrows; i++) {
1341               if (posSlack[i] >= 0 || negSlack[i] >= 0) {
1342                    int iCol;
1343                    int nPlus = 0;
1344                    int nMinus = 0;
1345                    bool possible = true;
1346                    // Get sum
1347                    double sum = 0.0;
1348                    iCol = posSlack[i];
1349                    while (iCol >= 0) {
1350                         double value = element[columnStart[iCol]];
1351                         sum += value * colsol[iCol];
1352                         if (lower[iCol]) {
1353                              possible = false;
1354                              break;
1355                         } else {
1356                              nPlus++;
1357                         }
1358                         iCol = nextSlack[iCol];
1359                    }
1360                    iCol = negSlack[i];
1361                    while (iCol >= 0) {
1362                         double value = -element[columnStart[iCol]];
1363                         sum -= value * colsol[iCol];
1364                         if (lower[iCol]) {
1365                              possible = false;
1366                              break;
1367                         } else {
1368                              nMinus++;
1369                         }
1370                         iCol = nextSlack[iCol];
1371                    }
1372                    //printf("%d plus, %d minus",nPlus,nMinus);
1373                    //printf("\n");
1374                    if ((rowsol[i] - rowlower[i] < 1.0e-7 ||
1375                              rowupper[i] - rowsol[i] < 1.0e-7) &&
1376                              nPlus + nMinus < 2)
1377                         possible = false;
1378                    if (possible) {
1379                         // Amount contributed by other varaibles
1380                         sum = rowsol[i] - sum;
1381                         double lo = rowlower[i];
1382                         if (lo > -1.0e20)
1383                              lo -= sum;
1384                         double up = rowupper[i];
1385                         if (up < 1.0e20)
1386                              up -= sum;
1387                         //printf("row bounds %g %g\n",lo,up);
1388                         if (0) {
1389                              double sum = 0.0;
1390                              double x = 0.0;
1391                              for (int k = 0; k < ncols; k++) {
1392                                   CoinBigIndex j;
1393                                   double value = colsol[k];
1394                                   x += value * cost[k];
1395                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1396                                        int irow = row[j];
1397                                        if (irow == i)
1398                                             sum += element[j] * value;
1399                                   }
1400                              }
1401                              printf("Before sum %g <= %g <= %g cost %.18g\n",
1402                                     rowlower[i], sum, rowupper[i], x);
1403                         }
1404                         // set all to zero
1405                         iCol = posSlack[i];
1406                         while (iCol >= 0) {
1407                              colsol[iCol] = 0.0;
1408                              iCol = nextSlack[iCol];
1409                         }
1410                         iCol = negSlack[i];
1411                         while (iCol >= 0) {
1412                              colsol[iCol] = 0.0;
1413                              iCol = nextSlack[iCol];
1414                         }
1415                         {
1416                              int iCol;
1417                              iCol = posSlack[i];
1418                              while (iCol >= 0) {
1419                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1420                                   //     iCol,element[columnStart[iCol]],
1421                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1422                                   iCol = nextSlack[iCol];
1423                              }
1424                              iCol = negSlack[i];
1425                              while (iCol >= 0) {
1426                                   //printf("col %d el %g sol %g bounds %g %g cost %g\n",
1427                                   //     iCol,element[columnStart[iCol]],
1428                                   //     colsol[iCol],lower[iCol],upper[iCol],cost[iCol]);
1429                                   iCol = nextSlack[iCol];
1430                              }
1431                         }
1432                         //printf("now what?\n");
1433                         int n = 0;
1434                         bool basic = false;
1435                         if (lo > 0.0) {
1436                              // Add in positive
1437                              iCol = posSlack[i];
1438                              while (iCol >= 0) {
1439                                   double value = element[columnStart[iCol]];
1440                                   ratio[n] = cost[iCol] / value;
1441                                   which[n++] = iCol;
1442                                   iCol = nextSlack[iCol];
1443                              }
1444                              CoinSort_2(ratio, ratio + n, which);
1445                              for (int i = 0; i < n; i++) {
1446                                   iCol = which[i];
1447                                   double value = element[columnStart[iCol]];
1448                                   if (lo >= upper[iCol]*value) {
1449                                        value *= upper[iCol];
1450                                        sum += value;
1451                                        lo -= value;
1452                                        colsol[iCol] = upper[iCol];
1453                                   } else {
1454                                        value = lo / value;
1455                                        sum += lo;
1456                                        lo = 0.0;
1457                                        colsol[iCol] = value;
1458                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1459                                        basic = true;
1460                                   }
1461                                   if (lo < 1.0e-7)
1462                                        break;
1463                              }
1464                         } else if (up < 0.0) {
1465                              // Use lo so coding is more similar
1466                              lo = -up;
1467                              // Add in negative
1468                              iCol = negSlack[i];
1469                              while (iCol >= 0) {
1470                                   double value = -element[columnStart[iCol]];
1471                                   ratio[n] = cost[iCol] / value;
1472                                   which[n++] = iCol;
1473                                   iCol = nextSlack[iCol];
1474                              }
1475                              CoinSort_2(ratio, ratio + n, which);
1476                              for (int i = 0; i < n; i++) {
1477                                   iCol = which[i];
1478                                   double value = -element[columnStart[iCol]];
1479                                   if (lo >= upper[iCol]*value) {
1480                                        value *= upper[iCol];
1481                                        sum += value;
1482                                        lo -= value;
1483                                        colsol[iCol] = upper[iCol];
1484                                   } else {
1485                                        value = lo / value;
1486                                        sum += lo;
1487                                        lo = 0.0;
1488                                        colsol[iCol] = value;
1489                                        model_->setColumnStatus(iCol, ClpSimplex::basic);
1490                                        basic = true;
1491                                   }
1492                                   if (lo < 1.0e-7)
1493                                        break;
1494                              }
1495                         }
1496                         if (0) {
1497                              double sum2 = 0.0;
1498                              double x = 0.0;
1499                              for (int k = 0; k < ncols; k++) {
1500                                   CoinBigIndex j;
1501                                   double value = colsol[k];
1502                                   x += value * cost[k];
1503                                   for (j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
1504                                        int irow = row[j];
1505                                        if (irow == i)
1506                                             sum2 += element[j] * value;
1507                                   }
1508                              }
1509                              printf("after sum %g <= %g <= %g cost %.18g (sum = %g)\n",
1510                                     rowlower[i], sum2, rowupper[i], x, sum);
1511                         }
1512                         rowsol[i] = sum;
1513                         if (basic) {
1514                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1515                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1516                              else
1517                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1518                         }
1519                    } else {
1520                         int n = 0;
1521                         int iCol;
1522                         iCol = posSlack[i];
1523                         while (iCol >= 0) {
1524                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1525                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1526                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1527                                   n++;
1528                              }
1529                              iCol = nextSlack[iCol];
1530                         }
1531                         iCol = negSlack[i];
1532                         while (iCol >= 0) {
1533                              if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1534                                        colsol[iCol] < upper[iCol] - 1.0e-8) {
1535                                   model_->setColumnStatus(iCol, ClpSimplex::basic);
1536                                   n++;
1537                              }
1538                              iCol = nextSlack[iCol];
1539                         }
1540                         if (n) {
1541                              if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1542                                   model_->setRowStatus(i, ClpSimplex::atLowerBound);
1543                              else
1544                                   model_->setRowStatus(i, ClpSimplex::atUpperBound);
1545#ifdef CLP_INVESTIGATE
1546                              if (n > 1)
1547                                   printf("%d basic on row %d!\n", n, i);
1548#endif
1549                         }
1550                    }
1551               }
1552          }
1553          delete [] ratio;
1554          delete [] which;
1555#else
1556          for (i = 0; i < nrows; i++) {
1557               int n = 0;
1558               int iCol;
1559               iCol = posSlack[i];
1560               while (iCol >= 0) {
1561                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1562                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1563                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1564                         n++;
1565                    }
1566                    iCol = nextSlack[iCol];
1567               }
1568               iCol = negSlack[i];
1569               while (iCol >= 0) {
1570                    if (colsol[iCol] > lower[iCol] + 1.0e-8 &&
1571                              colsol[iCol] < upper[iCol] - 1.0e-8) {
1572                         model_->setColumnStatus(iCol, ClpSimplex::basic);
1573                         n++;
1574                    }
1575                    iCol = nextSlack[iCol];
1576               }
1577               if (n) {
1578                    if (fabs(rowsol[i] - rowlower[i]) < fabs(rowsol[i] - rowupper[i]))
1579                         model_->setRowStatus(i, ClpSimplex::atLowerBound);
1580                    else
1581                         model_->setRowStatus(i, ClpSimplex::atUpperBound);
1582#ifdef CLP_INVESTIGATE
1583                    if (n > 1)
1584                         printf("%d basic on row %d!\n", n, i);
1585#endif
1586               }
1587          }
1588#endif
1589     }
1590     double maxmin;
1591     if (model_->getObjSense() == -1.0) {
1592          maxmin = -1.0;
1593     } else {
1594          maxmin = 1.0;
1595     }
1596     bool justValuesPass = majorIterations_ > 1000000;
1597     if (slackStart >= 0) {
1598          for (i = 0; i < nrows; i++) {
1599               model_->setRowStatus(i, ClpSimplex::superBasic);
1600          }
1601          for (i = slackStart; i < slackEnd; i++) {
1602               model_->setColumnStatus(i, ClpSimplex::basic);
1603          }
1604     } else {
1605          /* still try and put singletons rather than artificials in basis */
1606          int ninbas = 0;
1607          for (i = 0; i < nrows; i++) {
1608               model_->setRowStatus(i, ClpSimplex::basic);
1609          }
1610          for (i = 0; i < ncols; i++) {
1611               if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
1612                    CoinBigIndex j = columnStart[i];
1613                    double value = element[j];
1614                    int irow = row[j];
1615                    double rlo = rowlower[irow];
1616                    double rup = rowupper[irow];
1617                    double clo = lower[i];
1618                    double cup = upper[i];
1619                    double csol = colsol[i];
1620                    /* adjust towards feasibility */
1621                    double move = 0.0;
1622                    if (rowsol[irow] > rup) {
1623                         move = (rup - rowsol[irow]) / value;
1624                         if (value > 0.0) {
1625                              /* reduce */
1626                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1627                         } else {
1628                              /* increase */
1629                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1630                         }
1631                    } else if (rowsol[irow] < rlo) {
1632                         move = (rlo - rowsol[irow]) / value;
1633                         if (value > 0.0) {
1634                              /* increase */
1635                              if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1636                         } else {
1637                              /* reduce */
1638                              if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1639                         }
1640                    } else {
1641                         /* move to improve objective */
1642                         if (cost[i]*maxmin > 0.0) {
1643                              if (value > 0.0) {
1644                                   move = (rlo - rowsol[irow]) / value;
1645                                   /* reduce */
1646                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1647                              } else {
1648                                   move = (rup - rowsol[irow]) / value;
1649                                   /* increase */
1650                                   if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1651                              }
1652                         } else if (cost[i]*maxmin < 0.0) {
1653                              if (value > 0.0) {
1654                                   move = (rup - rowsol[irow]) / value;
1655                                   /* increase */
1656                                   if (csol + move > cup) move = CoinMax(0.0, cup - csol);
1657                              } else {
1658                                   move = (rlo - rowsol[irow]) / value;
1659                                   /* reduce */
1660                                   if (csol + move < clo) move = CoinMin(0.0, clo - csol);
1661                              }
1662                         }
1663                    }
1664                    rowsol[irow] += move * value;
1665                    colsol[i] += move;
1666                    /* put in basis if row was artificial */
1667                    if (rup - rlo < 1.0e-7 && model_->getRowStatus(irow) == ClpSimplex::basic) {
1668                         model_->setRowStatus(irow, ClpSimplex::superBasic);
1669                         model_->setColumnStatus(i, ClpSimplex::basic);
1670                         ninbas++;
1671                    }
1672               }
1673          }
1674          /*printf("%d in basis\n",ninbas);*/
1675     }
1676     bool wantVector = false;
1677     if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
1678          // See if original wanted vector
1679          ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
1680          wantVector = clpMatrixO->wantsSpecialColumnCopy();
1681     }
1682     if (addAll < 3) {
1683          ClpPresolve pinfo;
1684          if (presolve) {
1685               if (allowInfeasible) {
1686                    // fix up so will be feasible
1687                    double * rhs = new double[nrows];
1688                    memset(rhs, 0, nrows * sizeof(double));
1689                    model_->clpMatrix()->times(1.0, colsol, rhs);
1690                    double * rowupper = model_->rowUpper();
1691                    double * rowlower = model_->rowLower();
1692                    saveRowUpper = CoinCopyOfArray(rowupper, nrows);
1693                    saveRowLower = CoinCopyOfArray(rowlower, nrows);
1694                    double sum = 0.0;
1695                    for (i = 0; i < nrows; i++) {
1696                         if (rhs[i] > rowupper[i]) {
1697                              sum += rhs[i] - rowupper[i];
1698                              rowupper[i] = rhs[i];
1699                         }
1700                         if (rhs[i] < rowlower[i]) {
1701                              sum += rowlower[i] - rhs[i];
1702                              rowlower[i] = rhs[i];
1703                         }
1704                    }
1705                    COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum));
1706                    delete [] rhs;
1707               }
1708               saveModel = model_;
1709               pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
1710               model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1711          }
1712          if (model_) {
1713               if (!wantVector) {
1714                    //#define TWO_GOES
1715#ifdef ABC_INHERIT
1716#ifndef TWO_GOES
1717                    model_->dealWithAbc(1,justValuesPass ? 2 : 1);
1718#else
1719                    model_->dealWithAbc(1,1 + 11);
1720#endif
1721#else
1722#ifndef TWO_GOES
1723                    model_->primal(justValuesPass ? 2 : 1);
1724#else
1725                    model_->primal(1 + 11);
1726#endif
1727#endif
1728               } else {
1729                    ClpMatrixBase * matrix = model_->clpMatrix();
1730                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1731                    assert (clpMatrix);
1732                    clpMatrix->makeSpecialColumnCopy();
1733#ifdef ABC_INHERIT
1734                    model_->dealWithAbc(1,1);
1735#else
1736                    model_->primal(1);
1737#endif
1738                    clpMatrix->releaseSpecialColumnCopy();
1739               }
1740               if (presolve) {
1741                 model_->primal(1);
1742                    pinfo.postsolve(true);
1743                    delete model_;
1744                    model_ = saveModel;
1745                    saveModel = NULL;
1746               }
1747          } else {
1748               // not feasible
1749               addAll = 1;
1750               presolve = 0;
1751               model_ = saveModel;
1752               saveModel = NULL;
1753               if (justValuesPass)
1754#ifdef ABC_INHERIT
1755                    model_->dealWithAbc(1,2);
1756#else
1757                    model_->primal(2);
1758#endif
1759          }
1760          if (allowInfeasible) {
1761               CoinMemcpyN(saveRowUpper, nrows, model_->rowUpper());
1762               CoinMemcpyN(saveRowLower, nrows, model_->rowLower());
1763               delete [] saveRowUpper;
1764               delete [] saveRowLower;
1765               saveRowUpper = NULL;
1766               saveRowLower = NULL;
1767          }
1768          if (addAll < 2) {
1769               n = 0;
1770               if (!addAll ) {
1771                    /* could do scans to get a good number */
1772                    iteration = 1;
1773                    for (i = ordStart; i < ordEnd; i++) {
1774                         if (whenUsed[i] >= iteration) {
1775                              if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1776                                   n++;
1777                                   upper[i] = saveUpper[i];
1778                                   lower[i] = saveLower[i];
1779                              }
1780                         }
1781                    }
1782               } else {
1783                    for (i = ordStart; i < ordEnd; i++) {
1784                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1785                              n++;
1786                              upper[i] = saveUpper[i];
1787                              lower[i] = saveLower[i];
1788                         }
1789                    }
1790                    delete [] saveUpper;
1791                    delete [] saveLower;
1792                    saveUpper = NULL;
1793                    saveLower = NULL;
1794               }
1795#ifdef COIN_DEVELOP
1796               printf("Time so far %g, %d now added from previous iterations\n",
1797                      CoinCpuTime() - startTime, n);
1798#endif
1799               if (justValuesPass)
1800                    return;
1801               if (addAll)
1802                    presolve = 0;
1803               if (presolve) {
1804                    saveModel = model_;
1805                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1806               } else {
1807                    presolve = 0;
1808               }
1809               if (!wantVector) {
1810#ifdef ABC_INHERIT
1811                    model_->dealWithAbc(1,1);
1812#else
1813                    model_->primal(1);
1814#endif
1815               } else {
1816                    ClpMatrixBase * matrix = model_->clpMatrix();
1817                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1818                    assert (clpMatrix);
1819                    clpMatrix->makeSpecialColumnCopy();
1820#ifdef ABC_INHERIT
1821                    model_->dealWithAbc(1,1);
1822#else
1823                    model_->primal(1);
1824#endif
1825                    clpMatrix->releaseSpecialColumnCopy();
1826               }
1827               if (presolve) {
1828                    pinfo.postsolve(true);
1829                    delete model_;
1830                    model_ = saveModel;
1831                    saveModel = NULL;
1832               }
1833               if (!addAll) {
1834                    n = 0;
1835                    for (i = ordStart; i < ordEnd; i++) {
1836                         if (upper[i] - lower[i] < 1.0e-5 && saveUpper[i] - saveLower[i] > 1.0e-5) {
1837                              n++;
1838                              upper[i] = saveUpper[i];
1839                              lower[i] = saveLower[i];
1840                         }
1841                    }
1842                    delete [] saveUpper;
1843                    delete [] saveLower;
1844                    saveUpper = NULL;
1845                    saveLower = NULL;
1846#ifdef COIN_DEVELOP
1847                    printf("Time so far %g, %d now added from previous iterations\n",
1848                           CoinCpuTime() - startTime, n);
1849#endif
1850               }
1851               if (presolve) {
1852                    saveModel = model_;
1853                    model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
1854               } else {
1855                    presolve = 0;
1856               }
1857               if (!wantVector) {
1858#ifdef ABC_INHERIT
1859                    model_->dealWithAbc(1,1);
1860#else
1861                    model_->primal(1);
1862#endif
1863               } else {
1864                    ClpMatrixBase * matrix = model_->clpMatrix();
1865                    ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1866                    assert (clpMatrix);
1867                    clpMatrix->makeSpecialColumnCopy();
1868#ifdef ABC_INHERIT
1869                    model_->dealWithAbc(1,1);
1870#else
1871                    model_->primal(1);
1872#endif
1873                    clpMatrix->releaseSpecialColumnCopy();
1874               }
1875               if (presolve) {
1876                    pinfo.postsolve(true);
1877                    delete model_;
1878                    model_ = saveModel;
1879                    saveModel = NULL;
1880               }
1881          }
1882#ifdef COIN_DEVELOP
1883          printf("Total time in crossover %g\n", CoinCpuTime() - startTime);
1884#endif
1885          delete [] saveUpper;
1886          delete [] saveLower;
1887     }
1888#ifdef FEB_TRY
1889     model_->setSpecialOptions(saveOptions);
1890     model_->setPerturbation(savePerturbation);
1891#endif
1892     return ;
1893}
1894#endif
1895/*****************************************************************************/
1896
1897// Default contructor
1898Idiot::Idiot()
1899{
1900     model_ = NULL;
1901     maxBigIts_ = 3;
1902     maxIts_ = 5;
1903     logLevel_ = 1;
1904     logFreq_ = 100;
1905     maxIts2_ = 100;
1906     djTolerance_ = 1e-1;
1907     mu_ = 1e-4;
1908     drop_ = 5.0;
1909     exitDrop_ = -1.0e20;
1910     muFactor_ = 0.3333;
1911     stopMu_ = 1e-12;
1912     smallInfeas_ = 1e-1;
1913     reasonableInfeas_ = 1e2;
1914     muAtExit_ = 1.0e31;
1915     strategy_ = 8;
1916     lambdaIterations_ = 0;
1917     checkFrequency_ = 100;
1918     whenUsed_ = NULL;
1919     majorIterations_ = 30;
1920     exitFeasibility_ = -1.0;
1921     dropEnoughFeasibility_ = 0.02;
1922     dropEnoughWeighted_ = 0.01;
1923     // adjust
1924     double nrows = 10000.0;
1925     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
1926     baseIts = baseIts / 10;
1927     baseIts *= 10;
1928     maxIts2_ = 200 + baseIts + 5;
1929     maxIts2_ = 100;
1930     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
1931     lightWeight_ = 0;
1932}
1933// Constructor from model
1934Idiot::Idiot(OsiSolverInterface &model)
1935{
1936     model_ = & model;
1937     maxBigIts_ = 3;
1938     maxIts_ = 5;
1939     logLevel_ = 1;
1940     logFreq_ = 100;
1941     maxIts2_ = 100;
1942     djTolerance_ = 1e-1;
1943     mu_ = 1e-4;
1944     drop_ = 5.0;
1945     exitDrop_ = -1.0e20;
1946     muFactor_ = 0.3333;
1947     stopMu_ = 1e-12;
1948     smallInfeas_ = 1e-1;
1949     reasonableInfeas_ = 1e2;
1950     muAtExit_ = 1.0e31;
1951     strategy_ = 8;
1952     lambdaIterations_ = 0;
1953     checkFrequency_ = 100;
1954     whenUsed_ = NULL;
1955     majorIterations_ = 30;
1956     exitFeasibility_ = -1.0;
1957     dropEnoughFeasibility_ = 0.02;
1958     dropEnoughWeighted_ = 0.01;
1959     // adjust
1960     double nrows;
1961     if (model_)
1962          nrows = model_->getNumRows();
1963     else
1964          nrows = 10000.0;
1965     int baseIts = static_cast<int> (sqrt(static_cast<double>(nrows)));
1966     baseIts = baseIts / 10;
1967     baseIts *= 10;
1968     maxIts2_ = 200 + baseIts + 5;
1969     maxIts2_ = 100;
1970     reasonableInfeas_ = static_cast<double> (nrows) * 0.05;
1971     lightWeight_ = 0;
1972}
1973// Copy constructor.
1974Idiot::Idiot(const Idiot &rhs)
1975{
1976     model_ = rhs.model_;
1977     if (model_ && rhs.whenUsed_) {
1978          int numberColumns = model_->getNumCols();
1979          whenUsed_ = new int [numberColumns];
1980          CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
1981     } else {
1982          whenUsed_ = NULL;
1983     }
1984     djTolerance_ = rhs.djTolerance_;
1985     mu_ = rhs.mu_;
1986     drop_ = rhs.drop_;
1987     muFactor_ = rhs.muFactor_;
1988     stopMu_ = rhs.stopMu_;
1989     smallInfeas_ = rhs.smallInfeas_;
1990     reasonableInfeas_ = rhs.reasonableInfeas_;
1991     exitDrop_ = rhs.exitDrop_;
1992     muAtExit_ = rhs.muAtExit_;
1993     exitFeasibility_ = rhs.exitFeasibility_;
1994     dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
1995     dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
1996     maxBigIts_ = rhs.maxBigIts_;
1997     maxIts_ = rhs.maxIts_;
1998     majorIterations_ = rhs.majorIterations_;
1999     logLevel_ = rhs.logLevel_;
2000     logFreq_ = rhs.logFreq_;
2001     checkFrequency_ = rhs.checkFrequency_;
2002     lambdaIterations_ = rhs.lambdaIterations_;
2003     maxIts2_ = rhs.maxIts2_;
2004     strategy_ = rhs.strategy_;
2005     lightWeight_ = rhs.lightWeight_;
2006}
2007// Assignment operator. This copies the data
2008Idiot &
2009Idiot::operator=(const Idiot & rhs)
2010{
2011     if (this != &rhs) {
2012          delete [] whenUsed_;
2013          model_ = rhs.model_;
2014          if (model_ && rhs.whenUsed_) {
2015               int numberColumns = model_->getNumCols();
2016               whenUsed_ = new int [numberColumns];
2017               CoinMemcpyN(rhs.whenUsed_, numberColumns, whenUsed_);
2018          } else {
2019               whenUsed_ = NULL;
2020          }
2021          djTolerance_ = rhs.djTolerance_;
2022          mu_ = rhs.mu_;
2023          drop_ = rhs.drop_;
2024          muFactor_ = rhs.muFactor_;
2025          stopMu_ = rhs.stopMu_;
2026          smallInfeas_ = rhs.smallInfeas_;
2027          reasonableInfeas_ = rhs.reasonableInfeas_;
2028          exitDrop_ = rhs.exitDrop_;
2029          muAtExit_ = rhs.muAtExit_;
2030          exitFeasibility_ = rhs.exitFeasibility_;
2031          dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
2032          dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
2033          maxBigIts_ = rhs.maxBigIts_;
2034          maxIts_ = rhs.maxIts_;
2035          majorIterations_ = rhs.majorIterations_;
2036          logLevel_ = rhs.logLevel_;
2037          logFreq_ = rhs.logFreq_;
2038          checkFrequency_ = rhs.checkFrequency_;
2039          lambdaIterations_ = rhs.lambdaIterations_;
2040          maxIts2_ = rhs.maxIts2_;
2041          strategy_ = rhs.strategy_;
2042          lightWeight_ = rhs.lightWeight_;
2043     }
2044     return *this;
2045}
2046Idiot::~Idiot()
2047{
2048     delete [] whenUsed_;
2049}
Note: See TracBrowser for help on using the repository browser.