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

Last change on this file since 1879 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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