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

Last change on this file since 2324 was 2324, checked in by forrest, 19 months ago

fix for null problem

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